From 27a4713bd41ff34c8a5ef8fe0bf2eee0cf5deae6 Mon Sep 17 00:00:00 2001 From: Jonathan Robey Date: Fri, 15 Jul 2016 11:31:41 -0700 Subject: [PATCH] Create VectorTools:interpolate with mask --- doc/news/changes.h | 7 + include/deal.II/numerics/vector_tools.h | 6 +- .../deal.II/numerics/vector_tools.templates.h | 25 ++- .../numerics/vector_tools_interpolate.inst.in | 12 +- tests/dofs/interpolate_q_system_mask_01.cc | 163 +++++++++++++++++ .../dofs/interpolate_q_system_mask_01.output | 93 ++++++++++ tests/dofs/interpolate_q_system_mask_02.cc | 170 ++++++++++++++++++ .../dofs/interpolate_q_system_mask_02.output | 93 ++++++++++ 8 files changed, 557 insertions(+), 12 deletions(-) create mode 100644 tests/dofs/interpolate_q_system_mask_01.cc create mode 100644 tests/dofs/interpolate_q_system_mask_01.output create mode 100644 tests/dofs/interpolate_q_system_mask_02.cc create mode 100644 tests/dofs/interpolate_q_system_mask_02.output diff --git a/doc/news/changes.h b/doc/news/changes.h index b4a5e0d57e..101d667eb7 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -150,6 +150,13 @@ inconvenience this causes.

General

    +
  1. + Improved: VectorTools:interpolate now takes a ComponentMask to select the + components to interpolate. +
    + (Jonathan Robey, 2016/07/21) +
  2. +
  3. Improved: Split out pattern descriptions for LaTeX and Description ParameterHandler OutputStyles, and add better description text.
    diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index a6e8638a25..38cc561476 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -587,7 +587,8 @@ namespace VectorTools void interpolate (const Mapping &mapping, const DoFHandlerType &dof, const Function &function, - VectorType &vec); + VectorType &vec, + const ComponentMask &component_mask = ComponentMask()); /** * Calls the @p interpolate() function above with @@ -596,7 +597,8 @@ namespace VectorTools template void interpolate (const DoFHandlerType &dof, const Function &function, - VectorType &vec); + VectorType &vec, + const ComponentMask &component_mask = ComponentMask()); /** * Interpolate different finite element spaces. The interpolation of vector diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index cfe4a5756c..ee161501e7 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -80,14 +80,22 @@ namespace VectorTools void interpolate (const Mapping &mapping, const DoFHandlerType &dof, const Function &function, - VectorType &vec) + VectorType &vec, + const ComponentMask &component_mask) { typedef typename VectorType::value_type number; + Assert (component_mask.represents_n_components(dof.get_fe().n_components()), + ExcMessage("The number of components in the mask has to be either " + "zero or equal to the number of components in the finite " + "element.") ); + Assert (vec.size() == dof.n_dofs(), ExcDimensionMismatch (vec.size(), dof.n_dofs())); Assert (dof.get_fe().n_components() == function.n_components, ExcDimensionMismatch(dof.get_fe().n_components(), function.n_components)); + Assert (component_mask.n_selected_components(dof.get_fe().n_components()) > 0, + ComponentMask::ExcNoComponentSelected()); const hp::FECollection fe (dof.get_fe()); const unsigned int n_components = fe.n_components(); @@ -247,9 +255,12 @@ namespace VectorTools { const unsigned int component = fe[fe_index].system_to_component_index(i).first; - const unsigned int rep_dof=dof_to_rep_index_table[fe_index][i]; - vec(dofs_on_cell[i]) - = function_values_system[fe_index][rep_dof](component); + if (component_mask[component] == true) + { + const unsigned int rep_dof=dof_to_rep_index_table[fe_index][i]; + vec(dofs_on_cell[i]) + = function_values_system[fe_index][rep_dof](component); + } } } else @@ -277,13 +288,15 @@ namespace VectorTools template void interpolate (const DoFHandlerType &dof, const Function &function, - VectorType &vec) + VectorType &vec, + const ComponentMask &component_mask) { interpolate(StaticMappingQ1::mapping, dof, function, - vec); + vec, + component_mask); } diff --git a/source/numerics/vector_tools_interpolate.inst.in b/source/numerics/vector_tools_interpolate.inst.in index bc99300ad4..da50372b88 100644 --- a/source/numerics/vector_tools_interpolate.inst.in +++ b/source/numerics/vector_tools_interpolate.inst.in @@ -24,13 +24,15 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens (const Mapping &, const DoFHandler&, const Function&, - VEC&); + VEC&, + const ComponentMask&); template void interpolate (const DoFHandler&, const Function&, - VEC&); + VEC&, + const ComponentMask&); template void interpolate @@ -65,12 +67,14 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens (const Mapping&, const hp::DoFHandler&, const Function&, - VEC&); + VEC&, + const ComponentMask&); template void interpolate (const hp::DoFHandler&, const Function&, - VEC&); + VEC&, + const ComponentMask&); template void interpolate_based_on_material_id(const Mapping&, diff --git a/tests/dofs/interpolate_q_system_mask_01.cc b/tests/dofs/interpolate_q_system_mask_01.cc new file mode 100644 index 0000000000..2e2729075b --- /dev/null +++ b/tests/dofs/interpolate_q_system_mask_01.cc @@ -0,0 +1,163 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2006 - 2015 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 with a component mask works for +// FE_System(FE_Q(p)) elements correctly on a uniformly refined mesh for +// functions of degree q + +#include "../tests.h" +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +template +class F : public Function +{ +public: + F (const unsigned int q, const double adj) : Function(3), q(q), adj(adj) {} + + virtual void vector_value (const Point &p, + Vector &v) const + { + for (unsigned int c=0; c +void test () +{ + Triangulation triangulation; + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (3); + + for (unsigned int p=1; p<6-dim; ++p) + { + FE_Q fe_1(p); + FE_Q fe_2(p+1); + FESystem fe(fe_1, 2, + fe_2, 1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs (fe); + + // Use constant offset to distinguish between masks + const double adj1 = 0.3; + ComponentSelectFunction select_mask1(0,3); + ComponentMask mask1(3, false); + mask1.set(0, true); + + const double adj2 = 1.7; + ComponentSelectFunction select_mask2(std::make_pair(1,3),3); + ComponentMask mask2(3, false); + mask2.set(1, true); + mask2.set(2, true); + + Vector interpolant (dof_handler.n_dofs()); + Vector error (triangulation.n_active_cells()); + for (unsigned int q=0; q<=p+2; ++q) + { + // interpolate the function with mask 1 + VectorTools::interpolate (dof_handler, + F (q, adj1), + interpolant, + mask1); + + // interpolate the function with mask 2 + VectorTools::interpolate (dof_handler, + F (q, adj2), + interpolant, + mask2); + + // then compute the interpolation error for mask 1 + VectorTools::integrate_difference (dof_handler, + interpolant, + F (q, adj1), + error, + QGauss(q+2), + VectorTools::L2_norm, + &select_mask1); + if (q<=p) + Assert (error.l2_norm() < 1e-12*interpolant.l2_norm(), + ExcInternalError()); + + deallog << fe.get_name() << ", Mask 1, P_" << q + << ", rel. error=" << error.l2_norm() / interpolant.l2_norm() + << std::endl; + + // then compute the interpolation error for mask 2 + VectorTools::integrate_difference (dof_handler, + interpolant, + F (q, adj2), + error, + QGauss(q+2), + VectorTools::L2_norm, + &select_mask2); + if (q<=p) + Assert (error.l2_norm() < 1e-12*interpolant.l2_norm(), + ExcInternalError()); + + deallog << fe.get_name() << ", Mask 2, P_" << q + << ", rel. error=" << error.l2_norm() / interpolant.l2_norm() + << std::endl; + } + } +} + + + +int main () +{ + std::ofstream logfile("output"); + deallog << std::setprecision (3); + + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + test<1>(); + test<2>(); + test<3>(); +} + diff --git a/tests/dofs/interpolate_q_system_mask_01.output b/tests/dofs/interpolate_q_system_mask_01.output new file mode 100644 index 0000000000..c44928b7a0 --- /dev/null +++ b/tests/dofs/interpolate_q_system_mask_01.output @@ -0,0 +1,93 @@ + +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], Mask 1, P_2, rel. error=0.000124 +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], Mask 2, P_2, rel. error=0.000124 +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], Mask 1, P_3, rel. error=0.000296 +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], Mask 2, P_3, rel. error=0.000296 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 1, P_2, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 2, P_2, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 1, P_3, rel. error=2.31e-06 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 2, P_3, rel. error=2.31e-06 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 1, P_4, rel. error=6.92e-06 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 2, P_4, rel. error=6.92e-06 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 1, P_2, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 2, P_2, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 1, P_3, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 2, P_3, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 1, P_4, rel. error=5.66e-08 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 2, P_4, rel. error=5.66e-08 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 1, P_5, rel. error=2.03e-07 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 2, P_5, rel. error=2.03e-07 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 1, P_2, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 2, P_2, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 1, P_3, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 2, P_3, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 1, P_4, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 2, P_4, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 1, P_5, rel. error=1.50e-09 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 2, P_5, rel. error=1.50e-09 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 1, P_6, rel. error=6.27e-09 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 2, P_6, rel. error=6.27e-09 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], Mask 1, P_2, rel. error=4.17e-05 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], Mask 2, P_2, rel. error=4.17e-05 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], Mask 1, P_3, rel. error=9.68e-05 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], Mask 2, P_3, rel. error=9.68e-05 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 1, P_2, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 2, P_2, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 1, P_3, rel. error=4.87e-07 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 2, P_3, rel. error=4.87e-07 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 1, P_4, rel. error=1.46e-06 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 2, P_4, rel. error=1.46e-06 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 1, P_2, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 2, P_2, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 1, P_3, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 2, P_3, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 1, P_4, rel. error=1.02e-08 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 2, P_4, rel. error=1.02e-08 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 1, P_5, rel. error=3.68e-08 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 2, P_5, rel. error=3.68e-08 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], Mask 1, P_2, rel. error=1.26e-05 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], Mask 2, P_2, rel. error=1.26e-05 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], Mask 1, P_3, rel. error=2.89e-05 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], Mask 2, P_3, rel. error=2.89e-05 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 1, P_2, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 2, P_2, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 1, P_3, rel. error=1.03e-07 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 2, P_3, rel. error=1.03e-07 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 1, P_4, rel. error=3.08e-07 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 2, P_4, rel. error=3.08e-07 diff --git a/tests/dofs/interpolate_q_system_mask_02.cc b/tests/dofs/interpolate_q_system_mask_02.cc new file mode 100644 index 0000000000..b6b7e58ea1 --- /dev/null +++ b/tests/dofs/interpolate_q_system_mask_02.cc @@ -0,0 +1,170 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2006 - 2015 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 with a component mask works for +// FE_System(FE_Q(p)) elements correctly on an adaptively refined mesh for +// functions of degree q + +#include "../tests.h" +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +template +class F : public Function +{ +public: + F (const unsigned int q, const double adj) : Function(3), q(q), adj(adj) {} + + virtual void vector_value (const Point &p, + Vector &v) const + { + for (unsigned int c=0; c +void test () +{ + Triangulation triangulation; + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (1); + triangulation.begin_active()->set_refine_flag (); + triangulation.execute_coarsening_and_refinement (); + triangulation.refine_global (1); + + for (unsigned int p=1; p<6-dim; ++p) + { + FE_Q fe1(p); + FE_Q fe2(p+1); + FESystem fe(fe1, 2, fe2 , 1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs (fe); + + const double adj1 = 0.3; + ComponentSelectFunction select_mask1(0,3); + ComponentMask mask1(3, false); + mask1.set(0, true); + + const double adj2 = 1.7; + ComponentSelectFunction select_mask2(std::make_pair(1,3),3); + ComponentMask mask2(3, false); + mask2.set(1, true); + mask2.set(2, true); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + Vector interpolant (dof_handler.n_dofs()); + Vector error (triangulation.n_active_cells()); + for (unsigned int q=0; q<=p+2; ++q) + { + // interpolate the function + VectorTools::interpolate (dof_handler, + F (q, adj1), + interpolant, + mask1); + VectorTools::interpolate (dof_handler, + F (q, adj2), + interpolant, + mask2); + constraints.distribute (interpolant); + + + // then compute the interpolation error for mask 1 + VectorTools::integrate_difference (dof_handler, + interpolant, + F (q, adj1), + error, + QGauss(q+2), + VectorTools::L2_norm, + &select_mask1); + if (q<=p) + Assert (error.l2_norm() < 1e-12*interpolant.l2_norm(), + ExcInternalError()); + + deallog << fe.get_name() << ", Mask 1, P_" << q + << ", rel. error=" << error.l2_norm() / interpolant.l2_norm() + << std::endl; + + // then compute the interpolation error for mask 2 + VectorTools::integrate_difference (dof_handler, + interpolant, + F (q, adj2), + error, + QGauss(q+2), + VectorTools::L2_norm, + &select_mask2); + if (q<=p) + Assert (error.l2_norm() < 1e-12*interpolant.l2_norm(), + ExcInternalError()); + + deallog << fe.get_name() << ", Mask 2, P_" << q + << ", rel. error=" << error.l2_norm() / interpolant.l2_norm() + << std::endl; + } + } +} + + + +int main () +{ + std::ofstream logfile("output"); + deallog << std::setprecision (3); + + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + test<1>(); + test<2>(); + test<3>(); +} + diff --git a/tests/dofs/interpolate_q_system_mask_02.output b/tests/dofs/interpolate_q_system_mask_02.output new file mode 100644 index 0000000000..d44ce315d7 --- /dev/null +++ b/tests/dofs/interpolate_q_system_mask_02.output @@ -0,0 +1,93 @@ + +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], Mask 1, P_2, rel. error=0.000425 +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], Mask 2, P_2, rel. error=0.000425 +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], Mask 1, P_3, rel. error=0.00125 +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], Mask 2, P_3, rel. error=0.00125 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 1, P_2, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 2, P_2, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 1, P_3, rel. error=1.57e-05 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 2, P_3, rel. error=1.57e-05 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 1, P_4, rel. error=5.96e-05 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], Mask 2, P_4, rel. error=5.96e-05 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 1, P_2, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 2, P_2, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 1, P_3, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 2, P_3, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 1, P_4, rel. error=7.72e-07 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 2, P_4, rel. error=7.72e-07 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 1, P_5, rel. error=3.54e-06 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], Mask 2, P_5, rel. error=3.54e-06 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 1, P_2, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 2, P_2, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 1, P_3, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 2, P_3, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 1, P_4, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 2, P_4, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 1, P_5, rel. error=4.11e-08 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 2, P_5, rel. error=4.11e-08 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 1, P_6, rel. error=2.21e-07 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], Mask 2, P_6, rel. error=2.21e-07 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], Mask 1, P_2, rel. error=0.000222 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], Mask 2, P_2, rel. error=0.000222 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], Mask 1, P_3, rel. error=0.000561 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], Mask 2, P_3, rel. error=0.000561 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 1, P_2, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 2, P_2, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 1, P_3, rel. error=5.28e-06 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 2, P_3, rel. error=5.28e-06 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 1, P_4, rel. error=1.74e-05 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], Mask 2, P_4, rel. error=1.74e-05 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 1, P_2, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 2, P_2, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 1, P_3, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 2, P_3, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 1, P_4, rel. error=2.26e-07 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 2, P_4, rel. error=2.26e-07 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 1, P_5, rel. error=8.96e-07 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], Mask 2, P_5, rel. error=8.96e-07 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], Mask 1, P_2, rel. error=9.66e-05 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], Mask 2, P_2, rel. error=9.66e-05 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], Mask 1, P_3, rel. error=0.000232 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], Mask 2, P_3, rel. error=0.000232 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 1, P_0, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 2, P_0, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 1, P_1, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 2, P_1, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 1, P_2, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 2, P_2, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 1, P_3, rel. error=1.64e-06 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 2, P_3, rel. error=1.64e-06 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 1, P_4, rel. error=5.15e-06 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], Mask 2, P_4, rel. error=5.15e-06 -- 2.39.5