From de774856a0f72ab046be961a7acd78cfb8dcdeca Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Mon, 25 Sep 2017 17:43:38 -0500 Subject: [PATCH] reorganize generalized_interpolation tests, add a variant --- tests/numerics/generalized_interpolation.h | 97 +++++++++++++++++ .../numerics/generalized_interpolation_02.cc | 103 +++++------------- .../generalized_interpolation_02.output | 1 + .../numerics/generalized_interpolation_03.cc | 83 ++------------ .../numerics/generalized_interpolation_04.cc | 45 ++++++++ .../generalized_interpolation_04.output | 3 + 6 files changed, 179 insertions(+), 153 deletions(-) create mode 100644 tests/numerics/generalized_interpolation.h create mode 100644 tests/numerics/generalized_interpolation_04.cc create mode 100644 tests/numerics/generalized_interpolation_04.output diff --git a/tests/numerics/generalized_interpolation.h b/tests/numerics/generalized_interpolation.h new file mode 100644 index 0000000000..d1a37b5099 --- /dev/null +++ b/tests/numerics/generalized_interpolation.h @@ -0,0 +1,97 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +// This header contains a non-trivial testfunction F and a small test +// procedure shared among all generalized_interpolation_* tests. + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +template +class F : public Function +{ +public: + F(const unsigned int n_comp, const unsigned int q) + : Function(n_comp), 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 order_mapping, + bool distort_mesh, + bool print_function_values = false) +{ + 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); + + // Print function value in origin: + if (print_function_values) + { + Functions::FEFieldFunction f2(dof_handler, interpolant, mapping); + deallog << "Function value at (0.0,0.0): "; + for (unsigned int i = 0; i < fe.n_components(); ++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: + + 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; +} diff --git a/tests/numerics/generalized_interpolation_02.cc b/tests/numerics/generalized_interpolation_02.cc index 41d074be10..4c96f17711 100644 --- a/tests/numerics/generalized_interpolation_02.cc +++ b/tests/numerics/generalized_interpolation_02.cc @@ -16,96 +16,43 @@ // 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 #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; -} - +#include "generalized_interpolation.h" int main () { - deallog.depth_console(3); + initlog(); - 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), ConstantFunction<2>(1.0, 1), 1, false, true); + test<2>(FE_Q<2>(2), ConstantFunction<2>(1.0, 1), 2, false, true); + test<3>(FE_Q<3>(3), ConstantFunction<3>(1.0, 1), 3, false, true); - 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_Q<2>(1), ConstantFunction<2>(1.0, 1), 1, true, true); + test<2>(FE_Q<2>(2), ConstantFunction<2>(1.0, 1), 2, true, true); + test<3>(FE_Q<3>(3), ConstantFunction<3>(1.0, 1), 3, true, 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), ConstantFunction<2>(1.0, 2), 1, false, true); + test<2>(FE_RaviartThomas<2>(1), ConstantFunction<2>(1.0, 2), 2, false, true); + test<2>(FE_RaviartThomas<2>(2), ConstantFunction<2>(1.0, 2), 3, false, true); + test<3>(FE_RaviartThomas<3>(0), ConstantFunction<3>(1.0, 3), 1, false, true); + test<3>(FE_RaviartThomas<3>(1), ConstantFunction<3>(1.0, 3), 2, false, true); - 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); + test<2>(FE_RaviartThomas<2>(0), ConstantFunction<2>(1.0, 2), 1, true, true); + test<2>(FE_RaviartThomas<2>(1), ConstantFunction<2>(1.0, 2), 2, true, true); + test<2>(FE_RaviartThomas<2>(2), ConstantFunction<2>(1.0, 2), 3, true, 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); - - 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); + test<3>(FE_RaviartThomas<3>(1), ConstantFunction<3>(1.0, 3), 2, true, true); + + test<2>(FE_Nedelec<2>(0), ConstantFunction<2>(1.0, 2), 1, false, true); + test<2>(FE_Nedelec<2>(1), ConstantFunction<2>(1.0, 2), 2, false, true); + test<2>(FE_Nedelec<2>(2), ConstantFunction<2>(1.0, 2), 3, false, true); + test<2>(FE_Nedelec<2>(0), ConstantFunction<2>(1.0, 2), 1, true, true); + test<2>(FE_Nedelec<2>(1), ConstantFunction<2>(1.0, 2), 2, true, true); + test<2>(FE_Nedelec<2>(2), ConstantFunction<2>(1.0, 2), 3, true, true); } diff --git a/tests/numerics/generalized_interpolation_02.output b/tests/numerics/generalized_interpolation_02.output index 343f3d8675..29294e52f2 100644 --- a/tests/numerics/generalized_interpolation_02.output +++ b/tests/numerics/generalized_interpolation_02.output @@ -1,3 +1,4 @@ + DEAL::dim 2 FE_Q<2>(1) DEAL::Function value at (0.0,0.0): 1.00000 DEAL::Check projection property: 0.00000 diff --git a/tests/numerics/generalized_interpolation_03.cc b/tests/numerics/generalized_interpolation_03.cc index cbc2c24f88..385c67fe4f 100644 --- a/tests/numerics/generalized_interpolation_03.cc +++ b/tests/numerics/generalized_interpolation_03.cc @@ -14,90 +14,23 @@ // --------------------------------------------------------------------- // Check projection property of VectorTools::interpolate for -// Hdiv conforming spaces on something nontrivial. +// Hdiv and Hcurl 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; -} - +#include "generalized_interpolation.h" 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<2>(FE_RaviartThomas<2>(0), F<2>(2, 1), 1, false); + test<2>(FE_RaviartThomas<2>(1), F<2>(2, 0), 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); + test<3>(FE_RaviartThomas<3>(0), F<3>(3, 0), 1, false); + test<3>(FE_RaviartThomas<3>(1), F<3>(3, 0), 2, false); + test<3>(FE_RaviartThomas<3>(1), F<3>(3, 2), 2, false); } diff --git a/tests/numerics/generalized_interpolation_04.cc b/tests/numerics/generalized_interpolation_04.cc new file mode 100644 index 0000000000..cb3601b108 --- /dev/null +++ b/tests/numerics/generalized_interpolation_04.cc @@ -0,0 +1,45 @@ +// --------------------------------------------------------------------- +// +// 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 a complex, +// staggered system of Hdiv / Hcurl / L2 conforming spaces. + +#include +#include +#include +#include + +#include "../tests.h" +#include "generalized_interpolation.h" + +int main () +{ + initlog(); + + static const int dim = 2; + + FESystem fe(FE_RaviartThomas(1), + 1, + FE_Q(1), + 1, + FE_Nedelec(2), + 2, + FESystem(FE_Q(1), dim), + 1); + + const unsigned int n_comp = fe.n_components(); + + test<2>(fe, F<2>(n_comp, 1), 1, false); +} diff --git a/tests/numerics/generalized_interpolation_04.output b/tests/numerics/generalized_interpolation_04.output new file mode 100644 index 0000000000..59027680b3 --- /dev/null +++ b/tests/numerics/generalized_interpolation_04.output @@ -0,0 +1,3 @@ + +DEAL::dim 2 FESystem<2>[FE_RaviartThomas<2>(1)^2-FE_Q<2>(1)-FESystem<2>[FE_Nedelec<2>(2)^2]-FESystem<2>[FE_Q<2>(1)^2]] +DEAL::Check projection property: 20.8500 -- 2.39.5