Hide individual component masks

include/aspect/boundary_traction/interface.h

@@ -136,13 +136,17 @@ namespace aspect
         get_active_plugin_boundary_indicators() const;
 
         /**
-         * Return a list of component masks that indicate for
-         * each active plugin which components it is responsible for.
+         * Return a component masks that indicates for
+         * each boundary which components are prescribed by
+         * this manager class. All plugins that are responsible
+         * for this boundary use the same component mask.
          * The list of plugin objects can be
-         * requested by calling get_active_plugins().
+         * requested by calling get_active_plugins() and the
+         * list of boundaries they are responsible for is
+         * returnd by get_active_plugin_boundary_indicators().
          */
-        const std::vector<ComponentMask> &
-        get_active_plugin_component_masks() const;
+        ComponentMask
+        get_component_mask(const types::boundary_id boundary_id) const;
 
         /**
          * Declare the parameters of all known boundary traction plugins, as

include/aspect/boundary_velocity/interface.h

@@ -161,13 +161,17 @@ namespace aspect
         get_active_plugin_boundary_indicators() const;
 
         /**
-         * Return a list of component masks that indicate for
-         * each active plugin which components it is responsible for.
+         * Return a component masks that indicates for
+         * each boundary which components are prescribed by
+         * this manager class. All plugins that are responsible
+         * for this boundary use the same component mask.
          * The list of plugin objects can be
-         * requested by calling get_active_plugins().
+         * requested by calling get_active_plugins() and the
+         * list of boundaries they are responsible for is
+         * returnd by get_active_plugin_boundary_indicators().
          */
-        const std::vector<ComponentMask> &
-        get_active_plugin_component_masks() const;
+        ComponentMask
+        get_component_mask(const types::boundary_id boundary_id) const;
 
         /**
          * Return a set of boundary indicators for which boundary

source/boundary_traction/interface.cc

@@ -132,10 +132,30 @@ namespace aspect
 
 
     template <int dim>
-    const std::vector<ComponentMask> &
-    Manager<dim>::get_active_plugin_component_masks() const
+    ComponentMask
+    Manager<dim>::get_component_mask(const types::boundary_id boundary_id) const
     {
-      return component_masks;
+      Assert(prescribed_traction_boundary_indicators.find(boundary_id) != prescribed_traction_boundary_indicators.end(),
+             ExcMessage("The boundary traction manager class was asked for the "
+                        "component mask of boundary indicator <"
+                        +
+                        Utilities::int_to_string(boundary_id)
+                        +
+                        "> with symbolic name <"
+                        +
+                        this->get_geometry_model().translate_id_to_symbol_name(boundary_id)
+                        +
+                        ">, but this boundary is not part of the active boundary traction plugins."));
+
+      // Since all component masks of plugins at the same boundary are identical, we can use
+      // the component mask of the first plugin we find that is responsible for this boundary.
+      for (unsigned int i=0; i<boundary_indicators.size(); ++i)
+        if (boundary_indicators[i] == boundary_id)
+          return component_masks[i];
+
+      // We should never get here if plugins and boundary indicators were set up correctly.
+      AssertThrow(false, ExcInternalError());
+      return ComponentMask();
     }
 
 

source/boundary_velocity/interface.cc

@@ -129,10 +129,30 @@ namespace aspect
 
 
     template <int dim>
-    const std::vector<ComponentMask> &
-    Manager<dim>::get_active_plugin_component_masks() const
+    ComponentMask
+    Manager<dim>::get_component_mask(const types::boundary_id boundary_id) const
     {
-      return component_masks;
+      Assert(prescribed_velocity_boundary_indicators.find(boundary_id) != prescribed_velocity_boundary_indicators.end(),
+             ExcMessage("The boundary velocity manager class was asked for the "
+                        "component mask of boundary indicator <"
+                        +
+                        Utilities::int_to_string(boundary_id)
+                        +
+                        "> with symbolic name <"
+                        +
+                        this->get_geometry_model().translate_id_to_symbol_name(boundary_id)
+                        +
+                        ">, but this boundary is not part of the active boundary velocity plugins."));
+
+      // Since all component masks of plugins at the same boundary are identical, we can use
+      // the component mask of the first plugin we find that is responsible for this boundary.
+      for (unsigned int i=0; i<boundary_indicators.size(); ++i)
+        if (boundary_indicators[i] == boundary_id)
+          return component_masks[i];
+
+      // We should never get here if plugins and boundary indicators were set up correctly.
+      AssertThrow(false, ExcInternalError());
+      return ComponentMask();
     }
 
 

source/simulator/assemblers/stokes.cc

@@ -939,10 +939,9 @@ namespace aspect
 
       const typename DoFHandler<dim>::face_iterator face = scratch.cell->face(scratch.face_number);
 
-      const auto &traction_bis = this->get_boundary_traction_manager().get_active_plugin_boundary_indicators();
+      const auto &traction_bis = this->get_boundary_traction_manager().get_prescribed_boundary_traction_indicators();
 
-      if (std::find(traction_bis.begin(), traction_bis.end(), face->boundary_id())
-          != traction_bis.end())
+      if (traction_bis.find(face->boundary_id()) != traction_bis.end())
         {
           for (unsigned int q=0; q<scratch.face_finite_element_values.n_quadrature_points; ++q)
             {
@@ -952,13 +951,15 @@ namespace aspect
                                      scratch.face_finite_element_values.quadrature_point(q),
                                      scratch.face_finite_element_values.normal_vector(q));
 
+              const double JxW = scratch.face_finite_element_values.JxW(q);
+
               for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
                 {
                   if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
                     {
                       data.local_rhs(i_stokes) += scratch.face_finite_element_values[introspection.extractors.velocities].value(i,q) *
                                                   traction *
-                                                  scratch.face_finite_element_values.JxW(q);
+                                                  JxW;
                       ++i_stokes;
                     }
                   ++i;

source/simulator/assembly.cc

@@ -127,7 +127,7 @@ namespace aspect
                              " defined that handles this formulation."));
 
     // add the terms for traction boundary conditions
-    if (!boundary_traction_manager.get_active_plugins().empty())
+    if (!boundary_traction_manager.get_prescribed_boundary_traction_indicators().empty())
       {
         assemblers->stokes_system_on_boundary_face.push_back(
           std::make_unique<aspect::Assemblers::StokesBoundaryTraction<dim>>());
@@ -767,7 +767,7 @@ namespace aspect
       = (
           // see if we need to assemble traction boundary conditions.
           // only if so do we actually need to have an FEFaceValues object
-          !boundary_traction_manager.get_active_plugins().empty()
+          !boundary_traction_manager.get_prescribed_boundary_traction_indicators().empty()
           ?
           update_values |
           update_quadrature_points |

source/simulator/core.cc

@@ -1376,32 +1376,6 @@ namespace aspect
     // for the prescribed velocity fields
     boundary_velocity_manager.update();
 
-    // translate from component mask per plugin to component mask per boundary id
-    std::map<types::boundary_id, ComponentMask> boundary_component_masks;
-    for (unsigned int i=0; i<boundary_velocity_manager.get_active_plugins().size(); ++i)
-      {
-        const types::boundary_id boundary_id = boundary_velocity_manager.get_active_plugin_boundary_indicators()[i];
-        const auto &component_mask = boundary_velocity_manager.get_active_plugin_component_masks()[i];
-
-        const auto it = boundary_component_masks.find(boundary_id);
-        if (it != boundary_component_masks.end())
-          {
-            AssertThrow (it->second == component_mask,
-                         ExcMessage("Boundary indicator <"
-                                    +
-                                    Utilities::int_to_string(boundary_id)
-                                    +
-                                    "> with symbolic name <"
-                                    +
-                                    geometry_model->translate_id_to_symbol_name (boundary_id)
-                                    +
-                                    "> is listed as having different component masks "
-                                    "for different boundary velocity plugins. This is not allowed."));
-          }
-        else
-          boundary_component_masks[boundary_id] = component_mask;
-      }
-
     // put boundary conditions into constraints object for each boundary
     for (const auto boundary_id: boundary_velocity_manager.get_prescribed_boundary_velocity_indicators())
       {
@@ -1422,7 +1396,7 @@ namespace aspect
                                                   boundary_id,
                                                   vel,
                                                   constraints,
-                                                  boundary_component_masks.at(boundary_id));
+                                                  boundary_velocity_manager.get_component_mask(boundary_id));
       }
   }
 

source/simulator/helper_functions.cc

@@ -2387,38 +2387,42 @@ namespace aspect
     boundary_indicator_lists.emplace_back(boundary_velocity_manager.get_prescribed_boundary_velocity_indicators());
     boundary_indicator_lists.emplace_back(boundary_traction_manager.get_prescribed_boundary_traction_indicators());
 
-    // Make sure that each combination of boundary velocity and boundary traction plugin
+    // Make sure that each combination of boundary velocity and boundary traction condition
     // either refers to different boundary indicators or to different components
-    const auto &velocity_boundary_ids = boundary_velocity_manager.get_active_plugin_boundary_indicators();
-    const auto &traction_boundary_ids = boundary_traction_manager.get_active_plugin_boundary_indicators();
-
-    for (unsigned int i=0; i<velocity_boundary_ids.size(); ++i)
-      for (unsigned int j=0; j<traction_boundary_ids.size(); ++j)
-        {
-          if (velocity_boundary_ids[i] == traction_boundary_ids[j])
-            {
-              // if boundary ids are identical, make sure that the components are different
-              AssertThrow((boundary_velocity_manager.get_active_plugin_component_masks()[i] &
-                           boundary_traction_manager.get_active_plugin_component_masks()[j]) ==
-                          ComponentMask(introspection.n_components, false),
-                          ExcMessage("Boundary indicator <"
-                                     +
-                                     Utilities::int_to_string(velocity_boundary_ids[i])
-                                     +
-                                     "> with symbolic name <"
-                                     +
-                                     geometry_model->translate_id_to_symbol_name (velocity_boundary_ids[i])
-                                     +
-                                     "> is listed as having both "
-                                     "velocity and traction boundary conditions in the input file."));
-
-              // we have ensured the prescribed velocity and prescribed traction boundary conditions
-              // are compatible. In order to check them against the other boundary conditions, we
-              // need to remove the boundary indicator from one of the lists to make sure it only
-              // appears in one of them. We choose to remove the boundary indicator from the traction list.
-              boundary_indicator_lists[3].erase(traction_boundary_ids[j]);
-            }
-        }
+    for (const auto velocity_boundary_id: boundary_indicator_lists[2])
+      {
+        bool found_compatible_duplicate_boundary_id = false;
+        for (const auto traction_boundary_id: boundary_indicator_lists[3])
+          {
+            if (velocity_boundary_id == traction_boundary_id)
+              {
+                // if boundary ids are identical, make sure that the components are different
+                AssertThrow((boundary_velocity_manager.get_component_mask(velocity_boundary_id) &
+                             boundary_traction_manager.get_component_mask(traction_boundary_id)) ==
+                            ComponentMask(introspection.n_components, false),
+                            ExcMessage("Boundary indicator <"
+                                       +
+                                       Utilities::int_to_string(velocity_boundary_id)
+                                       +
+                                       "> with symbolic name <"
+                                       +
+                                       geometry_model->translate_id_to_symbol_name (velocity_boundary_id)
+                                       +
+                                       "> is listed as having both "
+                                       "velocity and traction boundary conditions in the input file."));
+
+                found_compatible_duplicate_boundary_id = true;
+              }
+          }
+        // we have ensured the prescribed velocity and prescribed traction boundary conditions
+        // for the current boundary id are compatible. In order to check them against the other
+        // boundary conditions, we need to remove the boundary indicator from one of the lists
+        // to make sure it only appears in one of them. We choose to remove the boundary
+        // indicator from the traction list. We cannot do that in the loop above, because it
+        // invalidates the range of the loop.
+        if (found_compatible_duplicate_boundary_id)
+          boundary_indicator_lists[3].erase(velocity_boundary_id);
+      }
 
     // for each combination of velocity boundary indicator lists, make sure that the
     // intersection is empty

source/simulator/melt.cc

@@ -1034,10 +1034,9 @@ namespace aspect
       Assert(face_no != numbers::invalid_unsigned_int,ExcInternalError());
 
       const typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
-      const auto &traction_bis = this->get_boundary_traction_manager().get_active_plugin_boundary_indicators();
+      const auto &traction_bis = this->get_boundary_traction_manager().get_prescribed_boundary_traction_indicators();
 
-      if (std::find(traction_bis.begin(), traction_bis.end(), face->boundary_id())
-          != traction_bis.end())
+      if (traction_bis.find(face->boundary_id()) != traction_bis.end())
         {
           scratch.face_finite_element_values.reinit (cell, face_no);
 
@@ -1689,7 +1688,7 @@ namespace aspect
       std::make_unique<aspect::Assemblers::MeltStokesSystemBoundary<dim>>());
 
     // add the terms for traction boundary conditions
-    if (!this->get_boundary_traction_manager().get_active_plugins().empty())
+    if (!this->get_boundary_traction_manager().get_prescribed_boundary_traction_indicators().empty())
       {
         assemblers.stokes_system_on_boundary_face.push_back(
           std::make_unique<Assemblers::MeltBoundaryTraction<dim>> ());

source/simulator/newton.cc

@@ -97,7 +97,7 @@ namespace aspect
                              " defined that handles this formulation."));
 
     // add the terms for traction boundary conditions
-    if (!this->get_boundary_traction_manager().get_active_plugins().empty())
+    if (!this->get_boundary_traction_manager().get_prescribed_boundary_traction_indicators().empty())
       {
         assemblers.stokes_system_on_boundary_face.push_back(
           std::make_unique<aspect::Assemblers::StokesBoundaryTraction<dim>>());

source/simulator/stokes_matrix_free.cc

@@ -2487,24 +2487,23 @@ namespace aspect
       mg_constrained_dofs_A_block.initialize(dof_handler_v);
 
       std::set<types::boundary_id> dirichlet_boundary = sim.boundary_velocity_manager.get_zero_boundary_velocity_indicators();
-      for (unsigned int i=0; i<sim.boundary_velocity_manager.get_active_plugins().size(); ++i)
+      for (const auto boundary_id: sim.boundary_velocity_manager.get_prescribed_boundary_velocity_indicators())
         {
-          const types::boundary_id bdryid = sim.boundary_velocity_manager.get_active_plugin_boundary_indicators()[i];
-          const ComponentMask component_mask = sim.boundary_velocity_manager.get_active_plugin_component_masks()[i];
+          const ComponentMask component_mask = sim.boundary_velocity_manager.get_component_mask(boundary_id);
 
           if (component_mask != ComponentMask(sim.introspection.n_components, false))
             {
-              ComponentMask mask(fe_v.n_components(), false);
+              ComponentMask velocity_mask(fe_v.n_components(), false);
 
               for (unsigned int i=0; i<dim; ++i)
-                mask.set(i, component_mask[sim.introspection.component_indices.velocities[i]]);
+                velocity_mask.set(i, component_mask[sim.introspection.component_indices.velocities[i]]);
 
-              mg_constrained_dofs_A_block.make_zero_boundary_constraints(dof_handler_v, {bdryid}, mask);
+              mg_constrained_dofs_A_block.make_zero_boundary_constraints(dof_handler_v, {boundary_id}, velocity_mask);
             }
           else
             {
               // no mask given: add at the end
-              dirichlet_boundary.insert(bdryid);
+              dirichlet_boundary.insert(boundary_id);
             }
         }