From 70d7386e74afa2bddda9d5f60d82ef7e3cdcdcfc Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Wed, 6 Sep 2017 17:57:39 -0500 Subject: [PATCH] add tests --- .../numerics/generalized_interpolation_02.cc | 112 ++++++++++++++++++ .../generalized_interpolation_02.output | 45 +++++++ .../numerics/generalized_interpolation_03.cc | 103 ++++++++++++++++ .../generalized_interpolation_03.output | 12 ++ 4 files changed, 272 insertions(+) create mode 100644 tests/numerics/generalized_interpolation_02.cc create mode 100644 tests/numerics/generalized_interpolation_02.output create mode 100644 tests/numerics/generalized_interpolation_03.cc create mode 100644 tests/numerics/generalized_interpolation_03.output diff --git a/tests/numerics/generalized_interpolation_02.cc b/tests/numerics/generalized_interpolation_02.cc new file mode 100644 index 0000000000..a57205be73 --- /dev/null +++ b/tests/numerics/generalized_interpolation_02.cc @@ -0,0 +1,112 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Check that VectorTools::interpolate correctly recovers a +// constant vector field for H1, Hdiv and Hcurl conforming elements. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../tests.h" + +template +void test(const FiniteElement &fe, + const unsigned int n_comp, + const unsigned int order_mapping, + bool distort_mesh) +{ + deallog << "dim " << dim << " " << fe.get_name() << std::endl; + + Triangulation triangulation; + GridGenerator::hyper_cube(triangulation, -0.3, 0.7); + triangulation.refine_global(dim == 2 ? 2 : 1); + if (distort_mesh) + GridTools::distort_random(0.03, triangulation); + + ConstantFunction f (1.0, n_comp); + + MappingQ mapping(order_mapping); + + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + + Vector interpolant(dof_handler.n_dofs()); + VectorTools::interpolate(mapping, dof_handler, f, interpolant); + + // Check that the interpoland really returns 1.0... + + Functions::FEFieldFunction f2(dof_handler, interpolant, mapping); + deallog << "Function value at (0.0,0.0): "; + for (unsigned int i = 0; i < n_comp; ++i) + deallog << f2.value(Point(), i) << " "; + deallog << std::endl; + + // Check that VectorTools::interpolate is in fact a + // projection, i.e. applying the interpolation twice results in the same + // vector: + + Vector interpolant2(dof_handler.n_dofs()); + VectorTools::interpolate( + mapping, dof_handler, f2, interpolant2); + + interpolant2 -= interpolant; + deallog << "Check projection property: " << interpolant2.linfty_norm() + << std::endl; +} + + +int main () +{ + deallog.depth_console(3); + + test<2>(FE_Q<2>(1), 1, 1, false); + test<2>(FE_Q<2>(2), 1, 2, false); + test<3>(FE_Q<3>(3), 1, 3, false); + + test<2>(FE_Q<2>(1), 1, 1, true); + test<2>(FE_Q<2>(2), 1, 2, true); + test<3>(FE_Q<3>(3), 1, 3, true); + + test<2>(FE_RaviartThomas<2>(0), 2, 1, false); + test<2>(FE_RaviartThomas<2>(1), 2, 2, false); + test<2>(FE_RaviartThomas<2>(2), 2, 3, false); + test<3>(FE_RaviartThomas<3>(0), 3, 1, false); + test<3>(FE_RaviartThomas<3>(1), 3, 2, false); + + test<2>(FE_RaviartThomas<2>(0), 2, 1, true); + test<2>(FE_RaviartThomas<2>(1), 2, 2, true); + test<2>(FE_RaviartThomas<2>(2), 2, 3, true); + // lowest order RT in 3D does not contain constant 1 function on a + // distorted mesh. + test<3>(FE_RaviartThomas<3>(1), 3, 2, true); + + // FIXME: Reenable, when FE_Nedelec is fixed :-/ + // test<2>(FE_Nedelec<2>(0), 2, 1, false); + // test<2>(FE_Nedelec<2>(1), 2, 2, false); + // test<2>(FE_Nedelec<2>(2), 2, 3, false); + // test<2>(FE_Nedelec<2>(0), 2, 1, true); + // test<2>(FE_Nedelec<2>(1), 2, 2, true); + // test<2>(FE_Nedelec<2>(2), 2, 3, true); +} diff --git a/tests/numerics/generalized_interpolation_02.output b/tests/numerics/generalized_interpolation_02.output new file mode 100644 index 0000000000..0e70a57e3f --- /dev/null +++ b/tests/numerics/generalized_interpolation_02.output @@ -0,0 +1,45 @@ +DEAL::dim 2 FE_Q<2>(1) +DEAL::Function value at (0.0,0.0): 1.00000 +DEAL::Check projection property: 0.00000 +DEAL::dim 2 FE_Q<2>(2) +DEAL::Function value at (0.0,0.0): 1.00000 +DEAL::Check projection property: 0.00000 +DEAL::dim 3 FE_Q<3>(3) +DEAL::Function value at (0.0,0.0): 1.00000 +DEAL::Check projection property: 0.00000 +DEAL::dim 2 FE_Q<2>(1) +DEAL::Function value at (0.0,0.0): 1.00000 +DEAL::Check projection property: 0.00000 +DEAL::dim 2 FE_Q<2>(2) +DEAL::Function value at (0.0,0.0): 1.00000 +DEAL::Check projection property: 0.00000 +DEAL::dim 3 FE_Q<3>(3) +DEAL::Function value at (0.0,0.0): 1.00000 +DEAL::Check projection property: 0.00000 +DEAL::dim 2 FE_RaviartThomas<2>(0) +DEAL::Function value at (0.0,0.0): 1.00000 1.00000 +DEAL::Check projection property: 0.00000 +DEAL::dim 2 FE_RaviartThomas<2>(1) +DEAL::Function value at (0.0,0.0): 1.00000 1.00000 +DEAL::Check projection property: 0.00000 +DEAL::dim 2 FE_RaviartThomas<2>(2) +DEAL::Function value at (0.0,0.0): 1.00000 1.00000 +DEAL::Check projection property: 0.00000 +DEAL::dim 3 FE_RaviartThomas<3>(0) +DEAL::Function value at (0.0,0.0): 1.00000 1.00000 1.00000 +DEAL::Check projection property: 0.00000 +DEAL::dim 3 FE_RaviartThomas<3>(1) +DEAL::Function value at (0.0,0.0): 1.00000 1.00000 1.00000 +DEAL::Check projection property: 0.00000 +DEAL::dim 2 FE_RaviartThomas<2>(0) +DEAL::Function value at (0.0,0.0): 1.00000 1.00000 +DEAL::Check projection property: 0.00000 +DEAL::dim 2 FE_RaviartThomas<2>(1) +DEAL::Function value at (0.0,0.0): 1.00000 1.00000 +DEAL::Check projection property: 0.00000 +DEAL::dim 2 FE_RaviartThomas<2>(2) +DEAL::Function value at (0.0,0.0): 1.00000 1.00000 +DEAL::Check projection property: 0.00000 +DEAL::dim 3 FE_RaviartThomas<3>(1) +DEAL::Function value at (0.0,0.0): 1.00000 1.00000 1.00000 +DEAL::Check projection property: 0.00000 diff --git a/tests/numerics/generalized_interpolation_03.cc b/tests/numerics/generalized_interpolation_03.cc new file mode 100644 index 0000000000..cbc2c24f88 --- /dev/null +++ b/tests/numerics/generalized_interpolation_03.cc @@ -0,0 +1,103 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Check projection property of VectorTools::interpolate for +// Hdiv conforming spaces on something nontrivial. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../tests.h" + +template +class F : public Function +{ +public: + F(const unsigned int q) : Function(dim), q(q) {} + + virtual double value(const Point &p, const unsigned int) const + { + double v = 0; + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int i = 0; i <= q; ++i) + v += (d + 1) * (i + 1) * std::pow(p[d], 1. * i); + return v; + } + +private: + const unsigned int q; +}; + +template +void test(const FiniteElement &fe, + const T &f, + const unsigned int n_comp, + const unsigned int order_mapping, + bool distort_mesh) +{ + deallog << "dim " << dim << " " << fe.get_name() << std::endl; + + Triangulation triangulation; + GridGenerator::hyper_cube(triangulation, -0.3, 0.7); + triangulation.refine_global(dim == 2 ? 2 : 1); + if (distort_mesh) + GridTools::distort_random(0.03, triangulation); + + MappingQ mapping(order_mapping); + + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + + Vector interpolant(dof_handler.n_dofs()); + VectorTools::interpolate(mapping, dof_handler, f, interpolant); + + // Check that VectorTools::interpolate is in fact a + // projection, i.e. applying the interpolation twice results in the same + // vector: + + Functions::FEFieldFunction f2(dof_handler, interpolant, mapping); + + Vector interpolant2(dof_handler.n_dofs()); + VectorTools::interpolate( + mapping, dof_handler, f2, interpolant2); + + interpolant2 -= interpolant; + deallog << "Check projection property: " << interpolant2.linfty_norm() + << std::endl; +} + + +int main () +{ + deallog.depth_console(3); + + test<2>(FE_RaviartThomas<2>(0), F<2>(1), 2, 1, false); + test<2>(FE_RaviartThomas<2>(1), F<2>(0), 2, 2, false); + test<2>(FE_RaviartThomas<2>(1), F<2>(2), 2, 2, false); + + test<3>(FE_RaviartThomas<3>(0), F<3>(0), 3, 1, false); + test<3>(FE_RaviartThomas<3>(1), F<3>(0), 3, 2, false); + test<3>(FE_RaviartThomas<3>(1), F<3>(2), 3, 2, false); +} diff --git a/tests/numerics/generalized_interpolation_03.output b/tests/numerics/generalized_interpolation_03.output new file mode 100644 index 0000000000..c9c7ab1284 --- /dev/null +++ b/tests/numerics/generalized_interpolation_03.output @@ -0,0 +1,12 @@ +DEAL::dim 2 FE_RaviartThomas<2>(0) +DEAL::Check projection property: 0.00000 +DEAL::dim 2 FE_RaviartThomas<2>(1) +DEAL::Check projection property: 0.00000 +DEAL::dim 2 FE_RaviartThomas<2>(1) +DEAL::Check projection property: 0.00000 +DEAL::dim 3 FE_RaviartThomas<3>(0) +DEAL::Check projection property: 0.00000 +DEAL::dim 3 FE_RaviartThomas<3>(1) +DEAL::Check projection property: 0.00000 +DEAL::dim 3 FE_RaviartThomas<3>(1) +DEAL::Check projection property: 0.00000 -- 2.39.5