Refactored map_dofs_to_support_points to allow for the use of a mask and removed...

diff --git a/include/deal.II/dofs/dof_tools.h b/include/deal.II/dofs/dof_tools.h
index 25a07560c2c4f62f965a4f9993791d3f81c5dc4a..fcfe83f1816a6aff36b7a79a721cfb0f82bbfd3c 100644 (file)
--- a/include/deal.II/dofs/dof_tools.h
+++ b/include/deal.II/dofs/dof_tools.h
@@ -2310,7 +2310,8 @@ namespace DoFTools
   void
   map_dofs_to_support_points(const Mapping<dim, spacedim> &   mapping,
                              const DoFHandler<dim, spacedim> &dof_handler,
-                             std::vector<Point<spacedim>> &   support_points);
+                             std::vector<Point<spacedim>> &   support_points,
+                             const ComponentMask &mask = ComponentMask());
 
   /**
    * Same as the previous function but for the hp case.
@@ -2321,7 +2322,8 @@ namespace DoFTools
   map_dofs_to_support_points(
     const dealii::hp::MappingCollection<dim, spacedim> &mapping,
     const hp::DoFHandler<dim, spacedim> &               dof_handler,
-    std::vector<Point<spacedim>> &                      support_points);
+    std::vector<Point<spacedim>> &                      support_points,
+    const ComponentMask &                               mask = ComponentMask());
 
   /**
    * This function is a version of the above map_dofs_to_support_points
@@ -2349,70 +2351,16 @@ namespace DoFTools
    * @param support_points A map that for every locally relevant DoF index
    * contains the corresponding location in real space coordinates. Previous
    * content of this object is deleted in this function.
+   * @param component_mask An optional component mask that restricts the
+   * components from which the support points are extracted
    */
-  template <int dim, int spacedim>
-  void
-  map_dofs_to_support_points(
-    const Mapping<dim, spacedim> &                      mapping,
-    const DoFHandler<dim, spacedim> &                   dof_handler,
-    std::map<types::global_dof_index, Point<spacedim>> &support_points);
-
-  // Same as above but with support for a component mask
-
   template <int dim, int spacedim>
   void
   map_dofs_to_support_points(
     const Mapping<dim, spacedim> &                      mapping,
     const DoFHandler<dim, spacedim> &                   dof_handler,
     std::map<types::global_dof_index, Point<spacedim>> &support_points,
-    const ComponentMask &                               mask)
-  {
-    const auto &tria = dof_handler.get_triangulation();
-    const auto &fe   = dof_handler.get_fe();
-
-    // Check if the FE has support points
-    Assert(fe.has_support_points(),
-           typename FiniteElement<dim>::ExcFEHasNoSupportPoints());
-
-    // Create quadrature at support point
-    Quadrature<dim> quad(fe.get_unit_support_points());
-
-    // Now loop over all cells and enquire the support points on each
-    // of these. we use dummy quadrature formulas where the quadrature
-    // points are located at the unit support points to enquire the
-    // location of the support points in real space.
-    //
-    // The weights of the quadrature rule have been set to invalid
-    // values by the used constructor.
-    //
-    // The mask is used to assess if a component should be kept
-    // or not
-    FEValues<dim, spacedim> fe_values(mapping,
-                                      fe,
-                                      quad,
-                                      update_quadrature_points);
-
-    std::vector<types::global_dof_index> local_dof_indices;
-    for (const auto &cell : dof_handler.active_cell_iterators())
-      // only work on locally relevant cells
-      if (cell->is_artificial() == false)
-        {
-          fe_values.reinit(cell);
-
-          local_dof_indices.resize(cell->get_fe().dofs_per_cell);
-          cell->get_dof_indices(local_dof_indices);
-
-          const std::vector<Point<spacedim>> &points =
-            fe_values.get_quadrature_points();
-          for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell; ++i)
-            {
-              unsigned int dof_comp = fe.system_to_component_index(i).first;
-              // insert the values into the map
-              if (mask[dof_comp])
-                support_points[local_dof_indices[i]] = points[i];
-            }
-        }
-  }
+    const ComponentMask &                               mask = ComponentMask());
 
   /**
    * Same as the previous function but for the hp case.
@@ -2422,7 +2370,8 @@ namespace DoFTools
   map_dofs_to_support_points(
     const dealii::hp::MappingCollection<dim, spacedim> &mapping,
     const hp::DoFHandler<dim, spacedim> &               dof_handler,
-    std::map<types::global_dof_index, Point<spacedim>> &support_points);
+    std::map<types::global_dof_index, Point<spacedim>> &support_points,
+    const ComponentMask &                               mask = ComponentMask());
 
 
   /**

diff --git a/source/dofs/dof_tools.cc b/source/dofs/dof_tools.cc
index be74d9c8f29b96c99fb54644ba5af259ad20fdbe..1ffb78db1e6b2638411631a2735fcdbefda31307 100644 (file)
--- a/source/dofs/dof_tools.cc
+++ b/source/dofs/dof_tools.cc
@@ -2077,7 +2077,8 @@ namespace DoFTools
                                     DoFHandlerType::space_dimension> &mapping,
         const DoFHandlerType &                            dof_handler,
         std::map<types::global_dof_index,
-                 Point<DoFHandlerType::space_dimension>> &support_points)
+                 Point<DoFHandlerType::space_dimension>> &support_points,
+        const ComponentMask &                             in_mask)
       {
         const unsigned int dim      = DoFHandlerType::dimension;
         const unsigned int spacedim = DoFHandlerType::space_dimension;
@@ -2096,6 +2097,12 @@ namespace DoFTools
               fe_collection[fe_index].get_unit_support_points()));
           }
 
+        // Take care of components
+        const ComponentMask mask =
+          (in_mask.size() == 0 ?
+             ComponentMask(fe_collection.n_components(), true) :
+             in_mask);
+
         // Now loop over all cells and enquire the support points on each
         // of these. we use dummy quadrature formulas where the quadrature
         // points are located at the unit support points to enquire the
@@ -2126,8 +2133,14 @@ namespace DoFTools
               const std::vector<Point<spacedim>> &points =
                 fe_values.get_quadrature_points();
               for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell; ++i)
-                // insert the values into the map
-                support_points[local_dof_indices[i]] = points[i];
+                {
+                  unsigned int dof_comp =
+                    cell->get_fe().system_to_component_index(i).first;
+
+                  // insert the values into the map if it is a valid component
+                  if (mask[dof_comp])
+                    support_points[local_dof_indices[i]] = points[i];
+                }
             }
       }
 
@@ -2138,13 +2151,17 @@ namespace DoFTools
         const hp::MappingCollection<DoFHandlerType::dimension,
                                     DoFHandlerType::space_dimension> &mapping,
         const DoFHandlerType &                               dof_handler,
-        std::vector<Point<DoFHandlerType::space_dimension>> &support_points)
+        std::vector<Point<DoFHandlerType::space_dimension>> &support_points,
+        const ComponentMask &                                mask)
       {
         // get the data in the form of the map as above
         std::map<types::global_dof_index,
                  Point<DoFHandlerType::space_dimension>>
           x_support_points;
-        map_dofs_to_support_points(mapping, dof_handler, x_support_points);
+        map_dofs_to_support_points(mapping,
+                                   dof_handler,
+                                   x_support_points,
+                                   mask);
 
         // now convert from the map to the linear vector. make sure every
         // entry really appeared in the map
@@ -2152,6 +2169,7 @@ namespace DoFTools
           {
             Assert(x_support_points.find(i) != x_support_points.end(),
                    ExcInternalError());
+
             support_points[i] = x_support_points[i];
           }
       }
@@ -2162,7 +2180,8 @@ namespace DoFTools
   void
   map_dofs_to_support_points(const Mapping<dim, spacedim> &   mapping,
                              const DoFHandler<dim, spacedim> &dof_handler,
-                             std::vector<Point<spacedim>> &   support_points)
+                             std::vector<Point<spacedim>> &   support_points,
+                             const ComponentMask &            mask)
   {
     AssertDimension(support_points.size(), dof_handler.n_dofs());
     Assert((dynamic_cast<
@@ -2178,7 +2197,8 @@ namespace DoFTools
 
     internal::map_dofs_to_support_points(mapping_collection,
                                          dof_handler,
-                                         support_points);
+                                         support_points,
+                                         mask);
   }
 
 
@@ -2187,7 +2207,8 @@ namespace DoFTools
   map_dofs_to_support_points(
     const hp::MappingCollection<dim, spacedim> &mapping,
     const hp::DoFHandler<dim, spacedim> &       dof_handler,
-    std::vector<Point<spacedim>> &              support_points)
+    std::vector<Point<spacedim>> &              support_points,
+    const ComponentMask &                       mask)
   {
     AssertDimension(support_points.size(), dof_handler.n_dofs());
     Assert((dynamic_cast<
@@ -2199,7 +2220,10 @@ namespace DoFTools
 
     // Let the internal function do all the work, just make sure that it
     // gets a MappingCollection
-    internal::map_dofs_to_support_points(mapping, dof_handler, support_points);
+    internal::map_dofs_to_support_points(mapping,
+                                         dof_handler,
+                                         support_points,
+                                         mask);
   }
 
 
@@ -2208,7 +2232,8 @@ namespace DoFTools
   map_dofs_to_support_points(
     const Mapping<dim, spacedim> &                      mapping,
     const DoFHandler<dim, spacedim> &                   dof_handler,
-    std::map<types::global_dof_index, Point<spacedim>> &support_points)
+    std::map<types::global_dof_index, Point<spacedim>> &support_points,
+    const ComponentMask &                               mask)
   {
     support_points.clear();
 
@@ -2218,7 +2243,8 @@ namespace DoFTools
 
     internal::map_dofs_to_support_points(mapping_collection,
                                          dof_handler,
-                                         support_points);
+                                         support_points,
+                                         mask);
   }
 
 
@@ -2227,13 +2253,17 @@ namespace DoFTools
   map_dofs_to_support_points(
     const hp::MappingCollection<dim, spacedim> &        mapping,
     const hp::DoFHandler<dim, spacedim> &               dof_handler,
-    std::map<types::global_dof_index, Point<spacedim>> &support_points)
+    std::map<types::global_dof_index, Point<spacedim>> &support_points,
+    const ComponentMask &                               mask)
   {
     support_points.clear();
 
     // Let the internal function do all the work, just make sure that it
     // gets a MappingCollection
-    internal::map_dofs_to_support_points(mapping, dof_handler, support_points);
+    internal::map_dofs_to_support_points(mapping,
+                                         dof_handler,
+                                         support_points,
+                                         mask);
   }
 
   template <int spacedim>

diff --git a/source/dofs/dof_tools.inst.in b/source/dofs/dof_tools.inst.in
index bcde5579b4a5688ed78689919898594c54f3fd05..f87fe38a54cf94f42c1a8fa5e2b4045c8558e3a1 100644 (file)
--- a/source/dofs/dof_tools.inst.in
+++ b/source/dofs/dof_tools.inst.in
@@ -489,25 +489,29 @@ for (deal_II_dimension : DIMENSIONS)
     template void DoFTools::map_dofs_to_support_points<deal_II_dimension>(
       const Mapping<deal_II_dimension, deal_II_dimension> &,
       const DoFHandler<deal_II_dimension> &,
-      std::vector<Point<deal_II_dimension>> &);
+      std::vector<Point<deal_II_dimension>> &,
+      const ComponentMask &);
 
 
     template void DoFTools::map_dofs_to_support_points<deal_II_dimension>(
       const hp::MappingCollection<deal_II_dimension, deal_II_dimension> &,
       const hp::DoFHandler<deal_II_dimension> &,
-      std::vector<Point<deal_II_dimension>> &);
+      std::vector<Point<deal_II_dimension>> &,
+      const ComponentMask &);
 
 
     template void DoFTools::map_dofs_to_support_points<deal_II_dimension>(
       const Mapping<deal_II_dimension, deal_II_dimension> &,
       const DoFHandler<deal_II_dimension> &,
-      std::map<types::global_dof_index, Point<deal_II_dimension>> &);
+      std::map<types::global_dof_index, Point<deal_II_dimension>> &,
+      const ComponentMask &);
 
 
     template void DoFTools::map_dofs_to_support_points<deal_II_dimension>(
       const hp::MappingCollection<deal_II_dimension, deal_II_dimension> &,
       const hp::DoFHandler<deal_II_dimension> &,
-      std::map<types::global_dof_index, Point<deal_II_dimension>> &);
+      std::map<types::global_dof_index, Point<deal_II_dimension>> &,
+      const ComponentMask &);
 
 #if deal_II_dimension < 3
 
@@ -515,25 +519,29 @@ for (deal_II_dimension : DIMENSIONS)
                                                        deal_II_dimension + 1>(
       const Mapping<deal_II_dimension, deal_II_dimension + 1> &,
       const DoFHandler<deal_II_dimension, deal_II_dimension + 1> &,
-      std::vector<Point<deal_II_dimension + 1>> &);
+      std::vector<Point<deal_II_dimension + 1>> &,
+      const ComponentMask &);
 
     template void DoFTools::map_dofs_to_support_points<deal_II_dimension,
                                                        deal_II_dimension + 1>(
       const hp::MappingCollection<deal_II_dimension, deal_II_dimension + 1> &,
       const hp::DoFHandler<deal_II_dimension, deal_II_dimension + 1> &,
-      std::vector<Point<deal_II_dimension + 1>> &);
+      std::vector<Point<deal_II_dimension + 1>> &,
+      const ComponentMask &);
 
     template void DoFTools::map_dofs_to_support_points<deal_II_dimension,
                                                        deal_II_dimension + 1>(
       const Mapping<deal_II_dimension, deal_II_dimension + 1> &,
       const DoFHandler<deal_II_dimension, deal_II_dimension + 1> &,
-      std::map<types::global_dof_index, Point<deal_II_dimension + 1>> &);
+      std::map<types::global_dof_index, Point<deal_II_dimension + 1>> &,
+      const ComponentMask &);
 
     template void DoFTools::map_dofs_to_support_points<deal_II_dimension,
                                                        deal_II_dimension + 1>(
       const hp::MappingCollection<deal_II_dimension, deal_II_dimension + 1> &,
       const hp::DoFHandler<deal_II_dimension, deal_II_dimension + 1> &,
-      std::map<types::global_dof_index, Point<deal_II_dimension + 1>> &);
+      std::map<types::global_dof_index, Point<deal_II_dimension + 1>> &,
+      const ComponentMask &);
 
 
     template void DoFTools::count_dofs_per_block<
@@ -554,7 +562,7 @@ for (deal_II_dimension : DIMENSIONS)
 #if deal_II_dimension == 3
 
     template void DoFTools::map_dofs_to_support_points<1, 3>(
-      const Mapping<1, 3> &, const DoFHandler<1, 3> &, std::vector<Point<3>> &);
+      const Mapping<1, 3> &, const DoFHandler<1, 3> &, std::vector<Point<3>> &,const ComponentMask &);
 
     template void DoFTools::count_dofs_per_block<DoFHandler<1, 3>>(
       const DoFHandler<1, 3> &,