From eafb31352d15f280239bfad7693c44fb4ea870cd Mon Sep 17 00:00:00 2001
From: Timo Heister <timo.heister@gmail.com>
Date: Tue, 16 May 2017 10:48:01 -0400
Subject: [PATCH] move and rename extract_component_subset

---
 include/aspect/utilities.h | 13 ++++++
 source/simulator/core.cc   | 93 ++------------------------------------
 source/utilities.cc        | 76 +++++++++++++++++++++++++++++++
 3 files changed, 92 insertions(+), 90 deletions(-)

diff --git a/include/aspect/utilities.h b/include/aspect/utilities.h
index 6a2c0879163..23fc5380324 100644
--- a/include/aspect/utilities.h
+++ b/include/aspect/utilities.h
@@ -78,6 +78,19 @@ namespace aspect
                          const IndexSet &whole_set,
                          std::vector<IndexSet> &partitioned);
 
+
+    /**
+     * Returns an IndexSet that contains all locally active DoFs that belong to
+     * the given component_mask.
+     *
+     * This function should be moved into deal.II at some point.
+     */
+    template <int dim>
+    IndexSet extract_locally_active_dofs_with_component(const DoFHandler<dim> &dof_handler,
+                                                        const ComponentMask &component_mask);
+
+
+
     namespace Coordinates
     {
 
diff --git a/source/simulator/core.cc b/source/simulator/core.cc
index 0b665c71a9f..f057956f4f9 100644
--- a/source/simulator/core.cc
+++ b/source/simulator/core.cc
@@ -1427,94 +1427,7 @@ namespace aspect
   }
 
 
-  /**
-   * This is an internal deal.II function stolen from dof_tools.cc
-   */
-  template <int dim, int spacedim>
-  std::vector<unsigned char>
-  get_local_component_association (const FiniteElement<dim,spacedim>  &fe,
-                                   const ComponentMask        &component_mask)
-  {
-    std::vector<unsigned char> local_component_association (fe.dofs_per_cell,
-                                                            (unsigned char)(-1));
-
-    // compute the component each local dof belongs to.
-    // if the shape function is primitive, then this
-    // is simple and we can just associate it with
-    // what system_to_component_index gives us
-    for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-      if (fe.is_primitive(i))
-        local_component_association[i] =
-          fe.system_to_component_index(i).first;
-      else
-        // if the shape function is not primitive, then either use the
-        // component of the first nonzero component corresponding
-        // to this shape function (if the component is not specified
-        // in the component_mask), or the first component of this block
-        // that is listed in the component_mask (if the block this
-        // component corresponds to is indeed specified in the component
-        // mask)
-        {
-          const unsigned int first_comp =
-            fe.get_nonzero_components(i).first_selected_component();
-
-          if ((fe.get_nonzero_components(i)
-               &
-               component_mask).n_selected_components(fe.n_components()) == 0)
-            local_component_association[i] = first_comp;
-          else
-            // pick the component selected. we know from the previous 'if'
-            // that within the components that are nonzero for this
-            // shape function there must be at least one for which the
-            // mask is true, so we will for sure run into the break()
-            // at one point
-            for (unsigned int c=first_comp; c<fe.n_components(); ++c)
-              if (component_mask[c] == true)
-                {
-                  local_component_association[i] = c;
-                  break;
-                }
-        }
-
-    Assert (std::find (local_component_association.begin(),
-                       local_component_association.end(),
-                       (unsigned char)(-1))
-            ==
-            local_component_association.end(),
-            ExcInternalError());
 
-    return local_component_association;
-  }
-
-  /**
-   * Returns an IndexSet that contains all locally active DoFs that belong to
-   * the given component_mask.
-   *
-   * This function should be moved into deal.II at some point.
-   */
-  template <int dim>
-  IndexSet extract_component_subset(DoFHandler<dim> &dof_handler, const ComponentMask &component_mask)
-  {
-    std::vector<unsigned char> local_asoc =
-      get_local_component_association (dof_handler.get_fe(),
-                                       ComponentMask(dof_handler.get_fe().n_components(), true));
-
-    IndexSet ret(dof_handler.n_dofs());
-
-    unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
-    std::vector<types::global_dof_index> indices(dofs_per_cell);
-    for (typename DoFHandler<dim>::active_cell_iterator cell=dof_handler.begin_active();
-         cell!=dof_handler.end(); ++cell)
-      if (cell->is_locally_owned())
-        {
-          cell->get_dof_indices(indices);
-          for (unsigned int i=0; i<dofs_per_cell; ++i)
-            if (component_mask[local_asoc[i]])
-              ret.add_index(indices[i]);
-        }
-
-    return ret;
-  }
 
 
   template <int dim>
@@ -1541,17 +1454,17 @@ namespace aspect
         {
           // locally_owned_melt_pressure_dofs is only used if not using melt
           if (parameters.use_direct_stokes_solver)
-            introspection.index_sets.locally_owned_pressure_dofs = system_index_set & extract_component_subset(dof_handler, introspection.component_masks.pressure);
+            introspection.index_sets.locally_owned_pressure_dofs = system_index_set & Utilities::extract_locally_active_dofs_with_component(dof_handler, introspection.component_masks.pressure);
           else
             introspection.index_sets.locally_owned_pressure_dofs = introspection.index_sets.system_partitioning[introspection.block_indices.pressure];
         }
 
       if (parameters.include_melt_transport)
         {
-          introspection.index_sets.locally_owned_melt_pressure_dofs = system_index_set & extract_component_subset(dof_handler,
+          introspection.index_sets.locally_owned_melt_pressure_dofs = system_index_set & Utilities::extract_locally_active_dofs_with_component(dof_handler,
                                                                       introspection.variable("fluid pressure").component_mask|
                                                                       introspection.variable("compaction pressure").component_mask);
-          introspection.index_sets.locally_owned_fluid_pressure_dofs = system_index_set & extract_component_subset(dof_handler,
+          introspection.index_sets.locally_owned_fluid_pressure_dofs = system_index_set & Utilities::extract_locally_active_dofs_with_component(dof_handler,
                                                                        introspection.variable("fluid pressure").component_mask);
         }
 
diff --git a/source/utilities.cc b/source/utilities.cc
index 6c81d31acc8..406ccc7c7f0 100644
--- a/source/utilities.cc
+++ b/source/utilities.cc
@@ -74,6 +74,72 @@ namespace aspect
         }
     }
 
+
+    /**
+    * This is an internal deal.II function stolen from dof_tools.cc
+    *
+    * Return an array that for each dof on the reference cell lists the
+    * corresponding vector component.
+    */
+    template <int dim, int spacedim>
+    std::vector<unsigned char>
+    get_local_component_association (const FiniteElement<dim,spacedim>  &fe,
+                                     const ComponentMask        &component_mask)
+    {
+      std::vector<unsigned char> local_component_association (fe.dofs_per_cell,
+                                                              (unsigned char)(-1));
+
+      // compute the component each local dof belongs to.
+      // if the shape function is primitive, then this
+      // is simple and we can just associate it with
+      // what system_to_component_index gives us
+      for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+        {
+          // see the deal.II version if we ever need non-primitive FEs.
+          Assert (fe.is_primitive(i), ExcNotImplemented());
+          local_component_association[i] =
+            fe.system_to_component_index(i).first;
+
+          Assert (std::find (local_component_association.begin(),
+                             local_component_association.end(),
+                             (unsigned char)(-1))
+                  ==
+                  local_component_association.end(),
+                  ExcInternalError());
+        }
+
+
+      return local_component_association;
+    }
+
+
+    template <int dim>
+    IndexSet extract_locally_active_dofs_with_component(const DoFHandler<dim> &dof_handler,
+                                                        const ComponentMask &component_mask)
+    {
+      std::vector<unsigned char> local_asoc =
+        get_local_component_association (dof_handler.get_fe(),
+                                         ComponentMask(dof_handler.get_fe().n_components(), true));
+
+      IndexSet ret(dof_handler.n_dofs());
+
+      unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+      std::vector<types::global_dof_index> indices(dofs_per_cell);
+      for (typename DoFHandler<dim>::active_cell_iterator cell=dof_handler.begin_active();
+           cell!=dof_handler.end(); ++cell)
+        if (cell->is_locally_owned())
+          {
+            cell->get_dof_indices(indices);
+            for (unsigned int i=0; i<dofs_per_cell; ++i)
+              if (component_mask[local_asoc[i]])
+                ret.add_index(indices[i]);
+          }
+
+      return ret;
+    }
+
+
+
     namespace Coordinates
     {
 
@@ -2246,6 +2312,16 @@ namespace aspect
 
 
 // Explicit instantiations
+
+#define INSTANTIATE(dim) \
+  template \
+  IndexSet extract_locally_active_dofs_with_component(const DoFHandler<dim> &, \
+                                                      const ComponentMask &);
+
+    ASPECT_INSTANTIATE(INSTANTIATE)
+
+
+
     template class AsciiDataLookup<1>;
     template class AsciiDataLookup<2>;
     template class AsciiDataLookup<3>;