<h3>General</h3>
<ol>
+ <li>
+ Improved: VectorTools:interpolate now takes a ComponentMask to select the
+ components to interpolate.
+ <br>
+ (Jonathan Robey, 2016/07/21)
+ </li>
+
<li> Improved: Split out pattern descriptions for LaTeX and Description
ParameterHandler OutputStyles, and add better description text.
<br>
void interpolate (const Mapping<dim,spacedim> &mapping,
const DoFHandlerType<dim,spacedim> &dof,
const Function<spacedim,typename VectorType::value_type> &function,
- VectorType &vec);
+ VectorType &vec,
+ const ComponentMask &component_mask = ComponentMask());
/**
* Calls the @p interpolate() function above with
template <typename VectorType, typename DoFHandlerType>
void interpolate (const DoFHandlerType &dof,
const Function<DoFHandlerType::space_dimension,typename VectorType::value_type> &function,
- VectorType &vec);
+ VectorType &vec,
+ const ComponentMask &component_mask = ComponentMask());
/**
* Interpolate different finite element spaces. The interpolation of vector
void interpolate (const Mapping<dim,spacedim> &mapping,
const DoFHandlerType<dim,spacedim> &dof,
const Function<spacedim, typename VectorType::value_type> &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<dim,spacedim> fe (dof.get_fe());
const unsigned int n_components = fe.n_components();
{
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
template <typename VectorType, typename DoFHandlerType>
void interpolate (const DoFHandlerType &dof,
const Function<DoFHandlerType::space_dimension,typename VectorType::value_type> &function,
- VectorType &vec)
+ VectorType &vec,
+ const ComponentMask &component_mask)
{
interpolate(StaticMappingQ1<DoFHandlerType::dimension,
DoFHandlerType::space_dimension>::mapping,
dof,
function,
- vec);
+ vec,
+ component_mask);
}
(const Mapping<deal_II_dimension,deal_II_space_dimension> &,
const DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
const Function<deal_II_space_dimension, VEC::value_type>&,
- VEC&);
+ VEC&,
+ const ComponentMask&);
template
void interpolate
(const DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
const Function<deal_II_space_dimension, VEC::value_type>&,
- VEC&);
+ VEC&,
+ const ComponentMask&);
template
void interpolate
(const Mapping<deal_II_dimension>&,
const hp::DoFHandler<deal_II_dimension>&,
const Function<deal_II_dimension, VEC::value_type>&,
- VEC&);
+ VEC&,
+ const ComponentMask&);
template
void interpolate
(const hp::DoFHandler<deal_II_dimension>&,
const Function<deal_II_dimension, VEC::value_type>&,
- VEC&);
+ VEC&,
+ const ComponentMask&);
template
void interpolate_based_on_material_id(const Mapping<deal_II_dimension, deal_II_space_dimension>&,
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <fstream>
+#include <vector>
+
+
+template <int dim>
+class F : public Function<dim>
+{
+public:
+ F (const unsigned int q, const double adj) : Function<dim>(3), q(q), adj(adj) {}
+
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &v) const
+ {
+ for (unsigned int c=0; c<v.size(); ++c)
+ {
+ v(c) = 0;
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int i=0; i<=q; ++i)
+ v(c) += (d+1)*(i+1)*std::pow (p[d], 1.*i)+c+adj;
+ }
+ }
+
+private:
+ const unsigned int q;
+ const double adj;
+};
+
+
+
+template <int dim>
+void test ()
+{
+ Triangulation<dim> triangulation;
+ GridGenerator::hyper_cube (triangulation);
+ triangulation.refine_global (3);
+
+ for (unsigned int p=1; p<6-dim; ++p)
+ {
+ FE_Q<dim> fe_1(p);
+ FE_Q<dim> fe_2(p+1);
+ FESystem<dim> fe(fe_1, 2,
+ fe_2, 1);
+ DoFHandler<dim> dof_handler(triangulation);
+ dof_handler.distribute_dofs (fe);
+
+ // Use constant offset to distinguish between masks
+ const double adj1 = 0.3;
+ ComponentSelectFunction<dim> select_mask1(0,3);
+ ComponentMask mask1(3, false);
+ mask1.set(0, true);
+
+ const double adj2 = 1.7;
+ ComponentSelectFunction<dim> select_mask2(std::make_pair(1,3),3);
+ ComponentMask mask2(3, false);
+ mask2.set(1, true);
+ mask2.set(2, true);
+
+ Vector<double> interpolant (dof_handler.n_dofs());
+ Vector<float> 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<dim> (q, adj1),
+ interpolant,
+ mask1);
+
+ // interpolate the function with mask 2
+ VectorTools::interpolate (dof_handler,
+ F<dim> (q, adj2),
+ interpolant,
+ mask2);
+
+ // then compute the interpolation error for mask 1
+ VectorTools::integrate_difference (dof_handler,
+ interpolant,
+ F<dim> (q, adj1),
+ error,
+ QGauss<dim>(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<dim> (q, adj2),
+ error,
+ QGauss<dim>(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>();
+}
+
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <fstream>
+#include <vector>
+
+
+template <int dim>
+class F : public Function<dim>
+{
+public:
+ F (const unsigned int q, const double adj) : Function<dim>(3), q(q), adj(adj) {}
+
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &v) const
+ {
+ for (unsigned int c=0; c<v.size(); ++c)
+ {
+ v(c) = 0;
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int i=0; i<=q; ++i)
+ v(c) += (d+1)*(i+1)*std::pow (p[d], 1.*i)+c+adj;
+ }
+ }
+
+private:
+ const unsigned int q;
+ const double adj;
+};
+
+
+
+template <int dim>
+void test ()
+{
+ Triangulation<dim> 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<dim> fe1(p);
+ FE_Q<dim> fe2(p+1);
+ FESystem<dim> fe(fe1, 2, fe2 , 1);
+ DoFHandler<dim> dof_handler(triangulation);
+ dof_handler.distribute_dofs (fe);
+
+ const double adj1 = 0.3;
+ ComponentSelectFunction<dim> select_mask1(0,3);
+ ComponentMask mask1(3, false);
+ mask1.set(0, true);
+
+ const double adj2 = 1.7;
+ ComponentSelectFunction<dim> 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<double> interpolant (dof_handler.n_dofs());
+ Vector<float> error (triangulation.n_active_cells());
+ for (unsigned int q=0; q<=p+2; ++q)
+ {
+ // interpolate the function
+ VectorTools::interpolate (dof_handler,
+ F<dim> (q, adj1),
+ interpolant,
+ mask1);
+ VectorTools::interpolate (dof_handler,
+ F<dim> (q, adj2),
+ interpolant,
+ mask2);
+ constraints.distribute (interpolant);
+
+
+ // then compute the interpolation error for mask 1
+ VectorTools::integrate_difference (dof_handler,
+ interpolant,
+ F<dim> (q, adj1),
+ error,
+ QGauss<dim>(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<dim> (q, adj2),
+ error,
+ QGauss<dim>(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>();
+}
+
--- /dev/null
+
+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