]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Create VectorTools:interpolate with mask
authorJonathan Robey <class4kayaker@gmail.com>
Fri, 15 Jul 2016 18:31:41 +0000 (11:31 -0700)
committerJonathan Robey <class4kayaker@gmail.com>
Sat, 23 Jul 2016 23:54:46 +0000 (16:54 -0700)
doc/news/changes.h
include/deal.II/numerics/vector_tools.h
include/deal.II/numerics/vector_tools.templates.h
source/numerics/vector_tools_interpolate.inst.in
tests/dofs/interpolate_q_system_mask_01.cc [new file with mode: 0644]
tests/dofs/interpolate_q_system_mask_01.output [new file with mode: 0644]
tests/dofs/interpolate_q_system_mask_02.cc [new file with mode: 0644]
tests/dofs/interpolate_q_system_mask_02.output [new file with mode: 0644]

index b4a5e0d57e0e80dfea6f8e5b1b1e7c40d53ecb0b..101d667eb71e28b812b6b68fa8e603b6eb2cb416 100644 (file)
@@ -150,6 +150,13 @@ inconvenience this causes.
 <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>
index a6e8638a25f65249f3c0234c1a1602f03636345a..38cc56147618a8a13f69a41a03418a541a0f5f9f 100644 (file)
@@ -587,7 +587,8 @@ namespace VectorTools
   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
@@ -596,7 +597,8 @@ namespace VectorTools
   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
index cfe4a5756c8cb637ffc4326ac5489162938c52bc..ee161501e772f4489048a9c70359327571e0684a 100644 (file)
@@ -80,14 +80,22 @@ namespace VectorTools
   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();
@@ -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 <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);
   }
 
 
index bc99300ad40ba68584faaf8532533b7ce5a17da8..da50372b88f76fd53da7e69a995c88f2090d6f37 100644 (file)
@@ -24,13 +24,15 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens
         (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
@@ -65,12 +67,14 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens
         (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>&,
diff --git a/tests/dofs/interpolate_q_system_mask_01.cc b/tests/dofs/interpolate_q_system_mask_01.cc
new file mode 100644 (file)
index 0000000..2e27290
--- /dev/null
@@ -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 <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>();
+}
+
diff --git a/tests/dofs/interpolate_q_system_mask_01.output b/tests/dofs/interpolate_q_system_mask_01.output
new file mode 100644 (file)
index 0000000..c44928b
--- /dev/null
@@ -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 (file)
index 0000000..b6b7e58
--- /dev/null
@@ -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 <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>();
+}
+
diff --git a/tests/dofs/interpolate_q_system_mask_02.output b/tests/dofs/interpolate_q_system_mask_02.output
new file mode 100644 (file)
index 0000000..d44ce31
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.