VT::interpolate() for multigrid levels

diff --git a/doc/news/changes/minor/20240726Munch b/doc/news/changes/minor/20240726Munch
new file mode 100644 (file)
index 0000000..0540303
--- /dev/null
+++ b/doc/news/changes/minor/20240726Munch
@@ -0,0 +1,3 @@
+New: VectorTools::interpolate() can now also used on multigrid levels.
+<br>
+(Peter Munch, Richard Schussnig, 2024/07/26)

diff --git a/include/deal.II/numerics/vector_tools_interpolate.h b/include/deal.II/numerics/vector_tools_interpolate.h
index b4006dfbfa74f40f44695b65e43460b7811db2ec..684fd2b2d03df00976fc432e384d12c5fc61a305 100644 (file)
--- a/include/deal.II/numerics/vector_tools_interpolate.h
+++ b/include/deal.II/numerics/vector_tools_interpolate.h
@@ -77,7 +77,8 @@ namespace VectorTools
     const DoFHandler<dim, spacedim>                           &dof,
     const Function<spacedim, typename VectorType::value_type> &function,
     VectorType                                                &vec,
-    const ComponentMask &component_mask = {});
+    const ComponentMask &component_mask = {},
+    const unsigned int   level          = numbers::invalid_unsigned_int);
 
   /**
    * Same as above but in an hp-context.
@@ -91,7 +92,8 @@ namespace VectorTools
     const DoFHandler<dim, spacedim>                           &dof,
     const Function<spacedim, typename VectorType::value_type> &function,
     VectorType                                                &vec,
-    const ComponentMask &component_mask = {});
+    const ComponentMask &component_mask = {},
+    const unsigned int   level          = numbers::invalid_unsigned_int);
 
 
   /**
@@ -106,7 +108,8 @@ namespace VectorTools
     const DoFHandler<dim, spacedim>                           &dof,
     const Function<spacedim, typename VectorType::value_type> &function,
     VectorType                                                &vec,
-    const ComponentMask &component_mask = {});
+    const ComponentMask &component_mask = {},
+    const unsigned int   level          = numbers::invalid_unsigned_int);
 
   /**
    * Interpolate different finite element spaces. The interpolation of vector

diff --git a/include/deal.II/numerics/vector_tools_interpolate.templates.h b/include/deal.II/numerics/vector_tools_interpolate.templates.h
index 96a617522e723c0b7a15db8f5a8b78f4429f66db..f5fafb1488a592a442a442aef404f9fe134a4a65 100644 (file)
--- a/include/deal.II/numerics/vector_tools_interpolate.templates.h
+++ b/include/deal.II/numerics/vector_tools_interpolate.templates.h
@@ -231,7 +231,8 @@ namespace VectorTools
       const DoFHandler<dim, spacedim>            &dof_handler,
       T                                          &function,
       VectorType                                 &vec,
-      const ComponentMask                        &component_mask)
+      const ComponentMask                        &component_mask,
+      const unsigned int level = numbers::invalid_unsigned_int)
     {
       Assert(component_mask.represents_n_components(
                dof_handler.get_fe_collection().n_components()),
@@ -240,7 +241,10 @@ namespace VectorTools
                "zero or equal to the number of components in the finite "
                "element."));
 
-      AssertDimension(vec.size(), dof_handler.n_dofs());
+      if (level == numbers::invalid_unsigned_int)
+        AssertDimension(vec.size(), dof_handler.n_dofs());
+      else
+        AssertDimension(vec.size(), dof_handler.n_dofs(level));
 
       Assert(component_mask.n_selected_components(
                dof_handler.get_fe_collection().n_components()) > 0,
@@ -354,143 +358,152 @@ namespace VectorTools
       //
       // Now loop over all locally owned, active cells.
       //
+      const auto runner = [&](const auto &cell) {
+        const unsigned int fe_index = cell->active_fe_index();
+
+        // Do nothing if there are no local degrees of freedom.
+        if (fe[fe_index].n_dofs_per_cell() == 0)
+          return;
+
+        // Skip processing of the current cell if the function object is
+        // invalid. This is used by interpolate_by_material_id to skip
+        // interpolating over cells with unknown material id.
+        if (!function(cell))
+          return;
+
+        // Get transformed, generalized support points
+        fe_values.reinit(cell);
+        const std::vector<Point<spacedim>> &generalized_support_points =
+          fe_values.get_present_fe_values().get_quadrature_points();
+
+        // Get indices of the dofs on this cell
+        const auto n_dofs = fe[fe_index].n_dofs_per_cell();
+        dofs_on_cell.resize(n_dofs);
+        cell->get_active_or_mg_dof_indices(dofs_on_cell);
+
+        // Prepare temporary storage
+        auto &function_values = fe_function_values[fe_index];
+        auto &dof_values      = fe_dof_values[fe_index];
+
+        const auto n_components = fe[fe_index].n_components();
+        // Only resize (and create sample entry) if sizes do not match
+        if (function_values.size() != generalized_support_points.size())
+          function_values.resize(generalized_support_points.size(),
+                                 Vector<number>(n_components));
+
+        // Get all function values:
+        AssertDimension(n_components, function(cell)->n_components);
+        function(cell)->vector_value_list(generalized_support_points,
+                                          function_values);
+
+        // For the simple case with elements with support points, we will
+        // simply use the interpolated DoF values in the access loop further
+        // down. Otherwise, we have to transform all function values from
+        // the real cell back to the unit cell. We query the finite element
+        // for the correct transformation. Matters get a bit more
+        // complicated because we have to apply said transformation for
+        // every base element.
+        if (needs_expensive_algorithm[fe_index])
+          {
+            dof_values.resize(n_dofs);
+            const unsigned int offset =
+              apply_transform(fe[fe_index],
+                              /* starting_offset = */ 0,
+                              fe_values,
+                              function_values);
+            (void)offset;
+            Assert(offset == n_components, ExcInternalError());
+
+            FETools::convert_generalized_support_point_values_to_dof_values(
+              fe[fe_index], function_values, dof_values);
+          }
 
-      for (const auto &cell : dof_handler.active_cell_iterators())
-        {
-          // If this cell is not locally owned, do nothing.
-          if (!cell->is_locally_owned())
-            continue;
-
-          const unsigned int fe_index = cell->active_fe_index();
-
-          // Do nothing if there are no local degrees of freedom.
-          if (fe[fe_index].n_dofs_per_cell() == 0)
-            continue;
-
-          // Skip processing of the current cell if the function object is
-          // invalid. This is used by interpolate_by_material_id to skip
-          // interpolating over cells with unknown material id.
-          if (!function(cell))
-            continue;
-
-          // Get transformed, generalized support points
-          fe_values.reinit(cell);
-          const std::vector<Point<spacedim>> &generalized_support_points =
-            fe_values.get_present_fe_values().get_quadrature_points();
-
-          // Get indices of the dofs on this cell
-          const auto n_dofs = fe[fe_index].n_dofs_per_cell();
-          dofs_on_cell.resize(n_dofs);
-          cell->get_dof_indices(dofs_on_cell);
-
-          // Prepare temporary storage
-          auto &function_values = fe_function_values[fe_index];
-          auto &dof_values      = fe_dof_values[fe_index];
-
-          const auto n_components = fe[fe_index].n_components();
-          // Only resize (and create sample entry) if sizes do not match
-          if (function_values.size() != generalized_support_points.size())
-            function_values.resize(generalized_support_points.size(),
-                                   Vector<number>(n_components));
-
-          // Get all function values:
-          AssertDimension(n_components, function(cell)->n_components);
-          function(cell)->vector_value_list(generalized_support_points,
-                                            function_values);
-
-          // For the simple case with elements with support points, we will
-          // simply use the interpolated DoF values in the access loop further
-          // down. Otherwise, we have to transform all function values from
-          // the real cell back to the unit cell. We query the finite element
-          // for the correct transformation. Matters get a bit more
-          // complicated because we have to apply said transformation for
-          // every base element.
-          if (needs_expensive_algorithm[fe_index])
-            {
-              dof_values.resize(n_dofs);
-              const unsigned int offset =
-                apply_transform(fe[fe_index],
-                                /* starting_offset = */ 0,
-                                fe_values,
-                                function_values);
-              (void)offset;
-              Assert(offset == n_components, ExcInternalError());
-
-              FETools::convert_generalized_support_point_values_to_dof_values(
-                fe[fe_index], function_values, dof_values);
-            }
-
-          for (unsigned int i = 0; i < n_dofs; ++i)
-            {
-              const auto &nonzero_components =
-                fe[fe_index].get_nonzero_components(i);
-
-              // Figure out whether the component mask applies. We assume
-              // that we are allowed to set degrees of freedom if at least
-              // one of the components (of the dof) is selected.
-              bool selected = false;
-              for (unsigned int c = 0; c < nonzero_components.size(); ++c)
-                selected =
-                  selected || (nonzero_components[c] && component_mask[c]);
-
-              if (selected)
-                {
+        for (unsigned int i = 0; i < n_dofs; ++i)
+          {
+            const auto &nonzero_components =
+              fe[fe_index].get_nonzero_components(i);
+
+            // Figure out whether the component mask applies. We assume
+            // that we are allowed to set degrees of freedom if at least
+            // one of the components (of the dof) is selected.
+            bool selected = false;
+            for (unsigned int c = 0; c < nonzero_components.size(); ++c)
+              selected =
+                selected || (nonzero_components[c] && component_mask[c]);
+
+            if (selected)
+              {
 #ifdef DEBUG
-                  // make sure that all selected base elements are indeed
-                  // interpolatory
-
-                  if (const auto fe_system =
-                        dynamic_cast<const FESystem<dim> *>(&fe[fe_index]))
-                    {
-                      const auto index =
-                        fe_system->system_to_base_index(i).first.first;
-                      Assert(fe_system->base_element(index)
-                               .has_generalized_support_points(),
-                             ExcMessage("The component mask supplied to "
-                                        "VectorTools::interpolate selects a "
-                                        "non-interpolatory element."));
-                    }
+                // make sure that all selected base elements are indeed
+                // interpolatory
+
+                if (const auto fe_system =
+                      dynamic_cast<const FESystem<dim> *>(&fe[fe_index]))
+                  {
+                    const auto index =
+                      fe_system->system_to_base_index(i).first.first;
+                    Assert(fe_system->base_element(index)
+                             .has_generalized_support_points(),
+                           ExcMessage("The component mask supplied to "
+                                      "VectorTools::interpolate selects a "
+                                      "non-interpolatory element."));
+                  }
 #endif
 
-                  // Add local values to the global vectors
-                  if (needs_expensive_algorithm[fe_index])
-                    ::dealii::internal::ElementAccess<VectorType>::add(
-                      dof_values[i], dofs_on_cell[i], interpolation);
-                  else
-                    {
-                      const auto base_index =
-                        fe[fe_index].system_to_base_index(i);
-                      ::dealii::internal::ElementAccess<VectorType>::add(
-                        function_values[base_index.second]
-                                       [base_index.first.second],
-                        dofs_on_cell[i],
-                        interpolation);
-                    }
+                // Add local values to the global vectors
+                if (needs_expensive_algorithm[fe_index])
                   ::dealii::internal::ElementAccess<VectorType>::add(
-                    typename VectorType::value_type(1.0),
-                    dofs_on_cell[i],
-                    weights);
-                }
-              else
-                {
-                  // If a component is ignored, copy the dof values
-                  // from the vector "vec", but only if they are locally
-                  // available
-                  if (locally_owned_dofs.is_element(dofs_on_cell[i]))
-                    {
-                      const auto value =
-                        ::dealii::internal::ElementAccess<VectorType>::get(
-                          vec, dofs_on_cell[i]);
-                      ::dealii::internal::ElementAccess<VectorType>::add(
-                        value, dofs_on_cell[i], interpolation);
-                      ::dealii::internal::ElementAccess<VectorType>::add(
-                        typename VectorType::value_type(1.0),
-                        dofs_on_cell[i],
-                        weights);
-                    }
-                }
+                    dof_values[i], dofs_on_cell[i], interpolation);
+                else
+                  {
+                    const auto base_index =
+                      fe[fe_index].system_to_base_index(i);
+                    ::dealii::internal::ElementAccess<VectorType>::add(
+                      function_values[base_index.second]
+                                     [base_index.first.second],
+                      dofs_on_cell[i],
+                      interpolation);
+                  }
+                ::dealii::internal::ElementAccess<VectorType>::add(
+                  typename VectorType::value_type(1.0),
+                  dofs_on_cell[i],
+                  weights);
+              }
+            else
+              {
+                // If a component is ignored, copy the dof values
+                // from the vector "vec", but only if they are locally
+                // available
+                if (locally_owned_dofs.is_element(dofs_on_cell[i]))
+                  {
+                    const auto value =
+                      ::dealii::internal::ElementAccess<VectorType>::get(
+                        vec, dofs_on_cell[i]);
+                    ::dealii::internal::ElementAccess<VectorType>::add(
+                      value, dofs_on_cell[i], interpolation);
+                    ::dealii::internal::ElementAccess<VectorType>::add(
+                      typename VectorType::value_type(1.0),
+                      dofs_on_cell[i],
+                      weights);
+                  }
+              }
+          }
+      };
+
+      if (level == numbers::invalid_unsigned_int)
+        {
+          for (const auto &cell : dof_handler.active_cell_iterators())
+            {
+              if (cell->is_locally_owned())
+                runner(cell);
             }
-        } /* loop over dof_handler.active_cell_iterators() */
+        }
+      else
+        {
+          for (const auto &cell : dof_handler.mg_cell_iterators_on_level(level))
+            if (cell->is_locally_owned_on_level())
+              runner(cell);
+        }
 
       interpolation.compress(VectorOperation::add);
       weights.compress(VectorOperation::add);
@@ -527,22 +540,21 @@ namespace VectorTools
     const DoFHandler<dim, spacedim>                           &dof_handler,
     const Function<spacedim, typename VectorType::value_type> &function,
     VectorType                                                &vec,
-    const ComponentMask                                       &component_mask)
+    const ComponentMask                                       &component_mask,
+    const unsigned int                                         level)
   {
     AssertDimension(dof_handler.get_fe_collection().n_components(),
                     function.n_components);
 
     // Create a small lambda capture wrapping function and call the
     // internal implementation
-    const auto function_map =
-      [&function](
-        const typename DoFHandler<dim, spacedim>::active_cell_iterator &)
+    const auto function_map = [&function](const auto &)
       -> const Function<spacedim, typename VectorType::value_type> * {
       return &function;
     };
 
     internal::interpolate(
-      mapping, dof_handler, function_map, vec, component_mask);
+      mapping, dof_handler, function_map, vec, component_mask, level);
   }
 
 
@@ -554,13 +566,15 @@ namespace VectorTools
     const DoFHandler<dim, spacedim>                           &dof_handler,
     const Function<spacedim, typename VectorType::value_type> &function,
     VectorType                                                &vec,
-    const ComponentMask                                       &component_mask)
+    const ComponentMask                                       &component_mask,
+    const unsigned int                                         level)
   {
     interpolate(hp::MappingCollection<dim, spacedim>(mapping),
                 dof_handler,
                 function,
                 vec,
-                component_mask);
+                component_mask,
+                level);
   }
 
 
@@ -571,7 +585,8 @@ namespace VectorTools
     const DoFHandler<dim, spacedim>                           &dof,
     const Function<spacedim, typename VectorType::value_type> &function,
     VectorType                                                &vec,
-    const ComponentMask                                       &component_mask)
+    const ComponentMask                                       &component_mask,
+    const unsigned int                                         level)
   {
     AssertDimension(dof.get_fe_collection().n_components(),
                     function.n_components);
@@ -579,7 +594,8 @@ namespace VectorTools
                 dof,
                 function,
                 vec,
-                component_mask);
+                component_mask,
+                level);
   }
 
 

diff --git a/source/numerics/vector_tools_interpolate.inst.in b/source/numerics/vector_tools_interpolate.inst.in
index d0da45e0e836e7119d8b0a8babd0c1554846ad1d..b79665303fae49983d859ab3fc3f2a0bfdd328f2 100644 (file)
--- a/source/numerics/vector_tools_interpolate.inst.in
+++ b/source/numerics/vector_tools_interpolate.inst.in
@@ -28,7 +28,8 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS;
         const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
         const Function<deal_II_space_dimension, VEC::value_type> &,
         VEC &,
-        const ComponentMask &);
+        const ComponentMask &,
+        const unsigned int);
 
       template void
       interpolate(
@@ -36,14 +37,16 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS;
         const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
         const Function<deal_II_space_dimension, VEC::value_type> &,
         VEC &,
-        const ComponentMask &);
+        const ComponentMask &,
+        const unsigned int);
 
       template void
       interpolate(
         const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
         const Function<deal_II_space_dimension, VEC::value_type> &,
         VEC &,
-        const ComponentMask &);
+        const ComponentMask &,
+        const unsigned int);
 
       template void
       interpolate(

diff --git a/tests/multigrid/mg_data_out_04.cc b/tests/multigrid/mg_data_out_04.cc
index 3b31dc3ba7e736ff9e9e0e275a1014bcc43ae877..a73eff676691fad851dc64d5544bdc4dbfdfcd87 100644 (file)
--- a/tests/multigrid/mg_data_out_04.cc
+++ b/tests/multigrid/mg_data_out_04.cc
@@ -36,7 +36,7 @@
 
 template <int dim>
 void
-do_test()
+do_test(const bool interpolate_with_vector_tools)
 {
   parallel::distributed::Triangulation<dim> triangulation(
     MPI_COMM_WORLD,
@@ -87,6 +87,19 @@ do_test()
     transfer.build(dof_handler);
     transfer.interpolate_to_mg(dof_handler, dof_vector, global_dof_vector);
 
+    if (interpolate_with_vector_tools)
+      for (unsigned int level = 0; level < triangulation.n_global_levels();
+           ++level)
+        {
+          dof_vector[level] = 0.0;
+
+          VectorTools::interpolate(dof_handler,
+                                   Functions::SquareFunction<dim>(),
+                                   dof_vector[level],
+                                   {},
+                                   level);
+        }
+
     for (unsigned int level = 0; level < triangulation.n_global_levels();
          ++level)
       {
@@ -115,7 +128,9 @@ main(int argc, char **argv)
   Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
   MPILogInitAll                    log;
 
-  do_test<2>();
+  do_test<2>(true);
+  do_test<2>(true);
+
   //  do_test<3>();
   return 0;
 }

diff --git a/tests/multigrid/mg_data_out_04.with_p4est=true.mpirun=2.output b/tests/multigrid/mg_data_out_04.with_p4est=true.mpirun=2.output
index 568ac410708ca1e2aa40123eb3fe7bc415562e49..c7bf970230316d5e402a8f46a64b1a62a7f6874e 100644 (file)
--- a/tests/multigrid/mg_data_out_04.with_p4est=true.mpirun=2.output
+++ b/tests/multigrid/mg_data_out_04.with_p4est=true.mpirun=2.output
@@ -100,7 +100,167 @@ DEAL:0::* level 2
 0.500000 0.500000 0.500000 
 
 
+# This file was generated by the deal.II library.
+
+
+#
+# For a description of the GNUPLOT format see the GNUPLOT manual.
+#
+# <x> <y> <data> 
+0.00000 0.00000 0.00000 
+0.250000 0.00000 0.0625000 
+
+0.00000 0.250000 0.0625000 
+0.250000 0.250000 0.125000 
+
+
+0.250000 0.00000 0.0625000 
+0.500000 0.00000 0.250000 
+
+0.250000 0.250000 0.125000 
+0.500000 0.250000 0.312500 
+
+
+0.00000 0.250000 0.0625000 
+0.250000 0.250000 0.125000 
+
+0.00000 0.500000 0.250000 
+0.250000 0.500000 0.312500 
+
+
+0.250000 0.250000 0.125000 
+0.500000 0.250000 0.312500 
+
+0.250000 0.500000 0.312500 
+0.500000 0.500000 0.500000 
+
+
+DEAL:0::* level 0
+# This file was generated by the deal.II library.
+
+
+#
+# For a description of the GNUPLOT format see the GNUPLOT manual.
+#
+# <x> <y> <data> 
+0.00000 0.00000 0.00000 
+1.00000 0.00000 1.00000 
+
+0.00000 1.00000 1.00000 
+1.00000 1.00000 2.00000 
+
+
+DEAL:0::* level 1
+# This file was generated by the deal.II library.
+
+
+#
+# For a description of the GNUPLOT format see the GNUPLOT manual.
+#
+# <x> <y> <data> 
+0.00000 0.00000 0.00000 
+0.500000 0.00000 0.250000 
+
+0.00000 0.500000 0.250000 
+0.500000 0.500000 0.500000 
+
+
+DEAL:0::* level 2
+# This file was generated by the deal.II library.
+
+
+#
+# For a description of the GNUPLOT format see the GNUPLOT manual.
+#
+# <x> <y> <data> 
+0.00000 0.00000 0.00000 
+0.250000 0.00000 0.0625000 
+
+0.00000 0.250000 0.0625000 
+0.250000 0.250000 0.125000 
+
+
+0.250000 0.00000 0.0625000 
+0.500000 0.00000 0.250000 
+
+0.250000 0.250000 0.125000 
+0.500000 0.250000 0.312500 
+
+
+0.00000 0.250000 0.0625000 
+0.250000 0.250000 0.125000 
+
+0.00000 0.500000 0.250000 
+0.250000 0.500000 0.312500 
+
+
+0.250000 0.250000 0.125000 
+0.500000 0.250000 0.312500 
+
+0.250000 0.500000 0.312500 
+0.500000 0.500000 0.500000 
+
+
+
+# This file was generated by the deal.II library.
+
+
+#
+# For a description of the GNUPLOT format see the GNUPLOT manual.
+#
+# <x> <y> <data> 
+0.500000 0.00000 0.250000 
+1.00000 0.00000 1.00000 
+
+0.500000 0.500000 0.500000 
+1.00000 0.500000 1.25000 
+
+
+0.00000 0.500000 0.250000 
+0.500000 0.500000 0.500000 
+
+0.00000 1.00000 1.00000 
+0.500000 1.00000 1.25000 
+
+
+0.500000 0.500000 0.500000 
+1.00000 0.500000 1.25000 
+
+0.500000 1.00000 1.25000 
+1.00000 1.00000 2.00000 
+
+
+DEAL:1::* level 0
+DEAL:1::* level 1
+# This file was generated by the deal.II library.
+
+
+#
+# For a description of the GNUPLOT format see the GNUPLOT manual.
+#
+# <x> <y> <data> 
+0.500000 0.00000 0.250000 
+1.00000 0.00000 1.00000 
+
+0.500000 0.500000 0.500000 
+1.00000 0.500000 1.25000 
+
+
+0.00000 0.500000 0.250000 
+0.500000 0.500000 0.500000 
+
+0.00000 1.00000 1.00000 
+0.500000 1.00000 1.25000 
+
+
+0.500000 0.500000 0.500000 
+1.00000 0.500000 1.25000 
 
+0.500000 1.00000 1.25000 
+1.00000 1.00000 2.00000 
+
+
+DEAL:1::* level 2
 # This file was generated by the deal.II library.