From 96b009e16a5c6e1ea971b737af258b9cf9ee891e Mon Sep 17 00:00:00 2001
From: Rene Gassmoeller <rene.gassmoeller@mailbox.org>
Date: Tue, 22 Oct 2024 13:19:55 +0200
Subject: [PATCH] tmp velocity

---
 include/aspect/boundary_traction/interface.h |  12 +
 include/aspect/boundary_velocity/interface.h | 101 ++++++--
 include/aspect/parameters.h                  |   8 -
 source/boundary_traction/interface.cc        |  26 +-
 source/boundary_velocity/ascii_data.cc       |  15 +-
 source/boundary_velocity/interface.cc        | 116 ++++++---
 source/mesh_deformation/free_surface.cc      |  17 +-
 source/simulator/assembly.cc                 |   2 +-
 source/simulator/core.cc                     |  99 ++++---
 source/simulator/helper_functions.cc         | 255 ++++++++-----------
 source/simulator/solver_schemes.cc           |   7 +-
 source/simulator/stokes_matrix_free.cc       |  36 +--
 tests/ascii_boundary_member.cc               |  21 +-
 13 files changed, 375 insertions(+), 340 deletions(-)

diff --git a/include/aspect/boundary_traction/interface.h b/include/aspect/boundary_traction/interface.h
index 85a8cdd13d6..a1235e0b073 100644
--- a/include/aspect/boundary_traction/interface.h
+++ b/include/aspect/boundary_traction/interface.h
@@ -119,6 +119,13 @@ namespace aspect
         const std::map<types::boundary_id,std::vector<std::unique_ptr<BoundaryTraction::Interface<dim>>>> &
         get_active_boundary_traction_conditions () const;
 
+        /**
+         * Return a set of boundary indicators for which boundary
+         * tractions are prescribed.
+         */
+        const std::set<types::boundary_id> &
+        get_prescribed_boundary_traction_indicators () const;
+
         /**
          * Return a list of boundary indicators that indicate for
          * each active plugin which boundary id
@@ -220,6 +227,11 @@ namespace aspect
          */
         std::vector<ComponentMask> component_masks;
 
+        /**
+         * A set of boundary indicators, on which tractions are prescribed.
+         */
+        std::set<types::boundary_id> prescribed_traction_boundary_indicators;
+
         /**
          * A list of boundary traction objects that have been requested in the
          * parameter file.
diff --git a/include/aspect/boundary_velocity/interface.h b/include/aspect/boundary_velocity/interface.h
index 9527c9bf453..8f1ccf69fb5 100644
--- a/include/aspect/boundary_velocity/interface.h
+++ b/include/aspect/boundary_velocity/interface.h
@@ -79,28 +79,9 @@ namespace aspect
      * @ingroup BoundaryVelocities
      */
     template <int dim>
-    class Manager : public SimulatorAccess<dim>
+    class Manager : public Plugins::ManagerBase<Interface<dim>>, public SimulatorAccess<dim>
     {
       public:
-        /**
-         * Destructor. Made virtual since this class has virtual member
-         * functions.
-         */
-        ~Manager () override;
-
-        /**
-         * A function that is called at the beginning of each time step and
-         * calls the corresponding functions of all created plugins.
-         *
-         * The point of this function is to allow complex boundary velocity
-         * models to do an initialization step once at the beginning of each
-         * time step. An example would be a model that needs to call an
-         * external program to compute the velocity change at a boundary.
-         */
-        virtual
-        void
-        update ();
-
         /**
          * A function that calls the boundary_velocity functions of all the
          * individual boundary velocity objects and uses the stored operators
@@ -146,7 +127,11 @@ namespace aspect
          * If there are no prescribed boundary velocity plugins
          * for a particular boundary, this boundary identifier will not appear
          * in the map.
+         *
+         * @deprecated: This function will be removed. Use the function
+         * get_active_plugin_names() of the base class ManagerBase instead.
          */
+        DEAL_II_DEPRECATED
         const std::map<types::boundary_id, std::pair<std::string,std::vector<std::string>>> &
         get_active_boundary_velocity_names () const;
 
@@ -158,10 +143,39 @@ namespace aspect
          * boundary models for this boundary. If there are no prescribed
          * boundary velocity plugins for a particular boundary this boundary
          * identifier will not appear in the map.
+         *
+         * @deprecated: This function has been removed. Use the function
+         * get_active_plugins() of the base class ManagerBase instead.
          */
+        DEAL_II_DEPRECATED
         const std::map<types::boundary_id,std::vector<std::unique_ptr<BoundaryVelocity::Interface<dim>>>> &
         get_active_boundary_velocity_conditions () const;
 
+        /**
+         * Return a list of boundary indicators that indicate for
+         * each active plugin which boundary id
+         * it is responsible for. The list of active plugins can be
+         * requested by calling get_active_plugins().
+         */
+        const std::vector<types::boundary_id> &
+        get_active_plugin_boundary_indicators() const;
+
+        /**
+         * Return a list of component masks that indicate for
+         * each active plugin which components it is responsible for.
+         * The list of plugin objects can be
+         * requested by calling get_active_plugins().
+         */
+        const std::vector<ComponentMask> &
+        get_active_plugin_component_masks() const;
+
+        /**
+         * Return a set of boundary indicators for which boundary
+         * velocities are prescribed.
+         */
+        const std::set<types::boundary_id> &
+        get_prescribed_boundary_velocity_indicators () const;
+
         /**
          * Return a list of boundary ids on which the velocity is prescribed
          * to be zero (no-slip).
@@ -190,7 +204,7 @@ namespace aspect
          * then let these objects read their parameters as well.
          */
         void
-        parse_parameters (ParameterHandler &prm);
+        parse_parameters (ParameterHandler &prm) override;
 
         /**
          * Go through the list of all boundary velocity models that have been selected
@@ -200,9 +214,15 @@ namespace aspect
          *
          * This function can only be called if the given template type (the first template
          * argument) is a class derived from the Interface class in this namespace.
+         *
+         * @deprecated Instead of this function, use the
+         *   Plugins::ManagerBase::has_matching_active_plugin() and
+         *   Plugins::ManagerBase::get_matching_active_plugin() functions of the base
+         *   class of the current class.
          */
         template <typename BoundaryVelocityType,
                   typename = typename std::enable_if_t<std::is_base_of<Interface<dim>,BoundaryVelocityType>::value>>
+        DEAL_II_DEPRECATED
         bool
         has_matching_boundary_velocity_model () const;
 
@@ -216,9 +236,15 @@ namespace aspect
          *
          * This function can only be called if the given template type (the first template
          * argument) is a class derived from the Interface class in this namespace.
+         *
+         * @deprecated Instead of this function, use the
+         *   Plugins::ManagerBase::has_matching_active_plugin() and
+         *   Plugins::ManagerBase::get_matching_active_plugin() functions of the base
+         *   class of the current class.
          */
         template <typename BoundaryVelocityType,
                   typename = typename std::enable_if_t<std::is_base_of<Interface<dim>,BoundaryVelocityType>::value>>
+        DEAL_II_DEPRECATED
         const BoundaryVelocityType &
         get_matching_boundary_velocity_model () const;
 
@@ -245,9 +271,32 @@ namespace aspect
                         << arg1
                         << "> among the names of registered boundary velocity objects.");
       private:
+        /**
+         * A list of boundary indicators that indicate for
+         * each plugin in the list of plugin_objects which boundary id
+         * it is responsible for. By default each plugin
+         * is active for all boundaries, but this list
+         * can be modified by derived classes to limit the application
+         * of plugins to specific boundaries.
+         */
+        std::vector<types::boundary_id> boundary_indicators;
+
+        /**
+         * A list of boundary indicators that indicate for
+         * each plugin in the list of plugin_objects which components
+         * it is responsible for. By default each plugin
+         * is active for all components, but this list
+         * can be modified by derived classes to limit the application
+         * of plugins to specific boundaries.
+         */
+        std::vector<ComponentMask> component_masks;
+
         /**
          * A list of boundary velocity objects that have been requested in the
          * parameter file.
+         *
+         * @deprecated: This variable is no longer used, but needed to issue a proper
+         * error message in the function get_active_boundary_velocity_conditions().
          */
         std::map<types::boundary_id,std::vector<std::unique_ptr<BoundaryVelocity::Interface<dim>>>> boundary_velocity_objects;
 
@@ -258,9 +307,19 @@ namespace aspect
          * mapped to one of the plugins of velocity boundary conditions (e.g.
          * "function"). If the components string is empty, it is assumed the
          * plugins are used for all components.
+         *
+         * @deprecated: Remove this variable when the deprecated functions
+         * get_active_boundary_velocity_names and
+         * get_active_boundary_velocity_conditions are removed. Use the base class
+         * variable plugin_names instead.
          */
         std::map<types::boundary_id, std::pair<std::string,std::vector<std::string>>> boundary_velocity_indicators;
 
+        /**
+         * A set of boundary indicators, on which velocities are prescribed.
+         */
+        std::set<types::boundary_id> prescribed_velocity_boundary_indicators;
+
         /**
          * A set of boundary indicators, on which velocities are prescribed to
          * zero (no-slip).
diff --git a/include/aspect/parameters.h b/include/aspect/parameters.h
index 7fdfb7aa265..4c29d9e692c 100644
--- a/include/aspect/parameters.h
+++ b/include/aspect/parameters.h
@@ -637,14 +637,6 @@ namespace aspect
     bool                           enable_additional_stokes_rhs;
     bool                           enable_prescribed_dilation;
 
-    /**
-     * Map from boundary id to a pair "components", "traction boundary type",
-     * where components is of the format "[x][y][z]" and the traction type is
-     * mapped to one of the plugins of traction boundary conditions (e.g.
-     * "function")
-     */
-    std::map<types::boundary_id, std::pair<std::string,std::string>> prescribed_traction_boundary_indicators;
-
     /**
      * A set of boundary ids on which the boundary_heat_flux objects
      * will be applied.
diff --git a/source/boundary_traction/interface.cc b/source/boundary_traction/interface.cc
index 1917009295d..c234437feb5 100644
--- a/source/boundary_traction/interface.cc
+++ b/source/boundary_traction/interface.cc
@@ -113,6 +113,15 @@ namespace aspect
 
 
 
+    template <int dim>
+    const std::set<types::boundary_id> &
+    Manager<dim>::get_prescribed_boundary_traction_indicators () const
+    {
+      return prescribed_traction_boundary_indicators;
+    }
+
+
+
     template <int dim>
     const std::vector<types::boundary_id> &
     Manager<dim>::get_active_plugin_boundary_indicators() const
@@ -273,16 +282,17 @@ namespace aspect
 
             this->plugin_names.push_back(value);
             boundary_indicators.push_back(boundary_id);
+            prescribed_traction_boundary_indicators.insert(boundary_id);
 
-            const bool default_component_mask = comp.empty();
-            ComponentMask component_mask(dim, default_component_mask);
+            ComponentMask component_mask(this->introspection().n_components,
+                                         false);
 
-            if (comp.find('x') != std::string::npos)
-              component_mask.set(0,true);
-            if (comp.find('y') != std::string::npos)
-              component_mask.set(1,true);
-            if (dim == 3 && comp.find('z') != std::string::npos)
-              component_mask.set(2,true);
+            if (comp.empty() || comp.find('x') != std::string::npos)
+              component_mask.set(this->introspection().component_indices.velocities[0],true);
+            if (comp.empty() || comp.find('y') != std::string::npos)
+              component_mask.set(this->introspection().component_indices.velocities[1],true);
+            if (dim == 3 && (comp.empty() || comp.find('z') != std::string::npos))
+              component_mask.set(this->introspection().component_indices.velocities[2],true);
 
             component_masks.push_back(component_mask);
           }
diff --git a/source/boundary_velocity/ascii_data.cc b/source/boundary_velocity/ascii_data.cc
index 5e0843075e1..faab1e46066 100644
--- a/source/boundary_velocity/ascii_data.cc
+++ b/source/boundary_velocity/ascii_data.cc
@@ -38,14 +38,17 @@ namespace aspect
     void
     AsciiData<dim>::initialize ()
     {
-      for (const auto &p : this->get_boundary_velocity_manager().get_active_boundary_velocity_conditions())
+      unsigned int i=0;
+      for (const auto &plugin : this->get_boundary_velocity_manager().get_active_plugins())
         {
-          for (const auto &plugin : p.second)
-            if (plugin.get() == this)
-              boundary_ids.insert(p.first);
+          if (plugin.get() == this)
+            boundary_ids.insert(this->get_boundary_velocity_manager().get_active_plugin_boundary_indicators()[i]);
+
+          ++i;
         }
-      AssertThrow(*(boundary_ids.begin()) != numbers::invalid_boundary_id,
-                  ExcMessage("Did not find the boundary indicator for the prescribed data plugin."));
+
+      AssertThrow(boundary_ids.empty() == false,
+                  ExcMessage("Did not find the boundary indicator for the velocity ascii data plugin."));
 
       Utilities::AsciiDataBoundary<dim>::initialize(boundary_ids,
                                                     dim);
diff --git a/source/boundary_velocity/interface.cc b/source/boundary_velocity/interface.cc
index a80a9245ad4..12b9ad48ba8 100644
--- a/source/boundary_velocity/interface.cc
+++ b/source/boundary_velocity/interface.cc
@@ -36,24 +36,6 @@ namespace aspect
     // ------------------------------ Manager -----------------------------
     // -------------------------------- Deal with registering boundary_velocity models and automating
     // -------------------------------- their setup and selection at run time
-
-    template <int dim>
-    Manager<dim>::~Manager()
-      = default;
-
-
-
-    template <int dim>
-    void
-    Manager<dim>::update ()
-    {
-      for (const auto &boundary : boundary_velocity_objects)
-        for (const auto &p : boundary.second)
-          p->update();
-    }
-
-
-
     namespace
     {
       std::tuple
@@ -84,26 +66,38 @@ namespace aspect
     Manager<dim>::boundary_velocity (const types::boundary_id boundary_indicator,
                                      const Point<dim> &position) const
     {
-      typename std::map<types::boundary_id,std::vector<std::unique_ptr<BoundaryVelocity::Interface<dim>>>>::const_iterator boundary_plugins =
-        boundary_velocity_objects.find(boundary_indicator);
+      Tensor<1,dim> velocity;
 
-      Assert(boundary_plugins != boundary_velocity_objects.end(),
+      bool found_plugin = false;
+      unsigned int i=0;
+      for (const auto &plugin: this->plugin_objects)
+        {
+          if (boundary_indicators[i] == boundary_indicator)
+            {
+              found_plugin = true;
+              const Tensor<1,dim> plugin_velocity = plugin->boundary_velocity(boundary_indicator,
+                                                                              position);
+              for (unsigned int d=0; d<dim; ++d)
+                if (component_masks[i][d] == true)
+                  velocity[d] += plugin_velocity[d];
+            }
+
+          ++i;
+        }
+
+      (void) found_plugin;
+      Assert(found_plugin == true,
              ExcMessage("The boundary velocity manager class was asked for the "
                         "boundary velocity at a boundary that contains no active "
                         "boundary velocity plugin."));
 
-      Tensor<1,dim> velocity = Tensor<1,dim>();
-
-      for (const auto &plugin : boundary_plugins->second)
-        velocity += plugin->boundary_velocity(boundary_indicator,
-                                              position);
-
       return velocity;
     }
 
 
 
     template <int dim>
+    DEAL_II_DEPRECATED
     const std::map<types::boundary_id, std::pair<std::string,std::vector<std::string>>> &
     Manager<dim>::get_active_boundary_velocity_names () const
     {
@@ -113,14 +107,45 @@ namespace aspect
 
 
     template <int dim>
+    DEAL_II_DEPRECATED
     const std::map<types::boundary_id,std::vector<std::unique_ptr<BoundaryVelocity::Interface<dim>>>> &
     Manager<dim>::get_active_boundary_velocity_conditions () const
     {
+      AssertThrow(false, ExcMessage("This function has been removed. Use the function "
+                                    "get_active_plugins() of the base class ManagerBase "
+                                    "instead."));
       return boundary_velocity_objects;
     }
 
 
 
+    template <int dim>
+    const std::vector<types::boundary_id> &
+    Manager<dim>::get_active_plugin_boundary_indicators() const
+    {
+      return boundary_indicators;
+    }
+
+
+
+    template <int dim>
+    const std::vector<ComponentMask> &
+    Manager<dim>::get_active_plugin_component_masks() const
+    {
+      return component_masks;
+    }
+
+
+
+    template <int dim>
+    const std::set<types::boundary_id> &
+    Manager<dim>::get_prescribed_boundary_velocity_indicators () const
+    {
+      return prescribed_velocity_boundary_indicators;
+    }
+
+
+
     template <int dim>
     const std::set<types::boundary_id> &
     Manager<dim>::get_zero_boundary_velocity_indicators () const
@@ -349,27 +374,40 @@ namespace aspect
               {
                 boundary_velocity_indicators[boundary_id] = std::make_pair(comp,std::vector<std::string>(1,value));
               }
+
+            this->plugin_names.push_back(value);
+            boundary_indicators.push_back(boundary_id);
+            prescribed_velocity_boundary_indicators.insert(boundary_id);
+
+            ComponentMask component_mask(this->introspection().n_components,
+                                         false);
+
+            if (comp.empty() || comp.find('x') != std::string::npos)
+              component_mask.set(this->introspection().component_indices.velocities[0],true);
+            if (comp.empty() || comp.find('y') != std::string::npos)
+              component_mask.set(this->introspection().component_indices.velocities[1],true);
+            if (dim == 3 && (comp.empty() || comp.find('z') != std::string::npos))
+              component_mask.set(this->introspection().component_indices.velocities[2],true);
+
+            component_masks.push_back(component_mask);
           }
       }
       prm.leave_subsection();
 
       // go through the list, create objects and let them parse
       // their own parameters
-      for (const auto &boundary_id : boundary_velocity_indicators)
+      for (const auto &plugin_name: this->plugin_names)
         {
-          for (const auto &name : boundary_id.second.second)
-            {
-              boundary_velocity_objects[boundary_id.first].push_back(
-                std::unique_ptr<Interface<dim>> (std::get<dim>(registered_plugins)
-                                                  .create_plugin (name,
-                                                                  "Boundary velocity::Model names")));
+          // create boundary traction objects
+          this->plugin_objects.push_back(std::get<dim>(registered_plugins)
+                                         .create_plugin (plugin_name,
+                                                         "Boundary velocity::Model names"));
 
-              if (SimulatorAccess<dim> *sim = dynamic_cast<SimulatorAccess<dim>*>(boundary_velocity_objects[boundary_id.first].back().get()))
-                sim->initialize_simulator (this->get_simulator());
+          if (SimulatorAccess<dim> *sim = dynamic_cast<SimulatorAccess<dim>*>(this->plugin_objects.back().get()))
+            sim->initialize_simulator (this->get_simulator());
 
-              boundary_velocity_objects[boundary_id.first].back()->parse_parameters (prm);
-              boundary_velocity_objects[boundary_id.first].back()->initialize ();
-            }
+          this->plugin_objects.back()->parse_parameters (prm);
+          this->plugin_objects.back()->initialize ();
         }
     }
 
diff --git a/source/mesh_deformation/free_surface.cc b/source/mesh_deformation/free_surface.cc
index 8d67a0f0557..b1ae199ac5e 100644
--- a/source/mesh_deformation/free_surface.cc
+++ b/source/mesh_deformation/free_surface.cc
@@ -50,20 +50,17 @@ namespace aspect
       std::set<types::boundary_id> velocity_boundary_indicators = this->get_boundary_velocity_manager().get_zero_boundary_velocity_indicators();
 
       // Get the tangential velocity boundary indicators
-      const std::set<types::boundary_id> tmp_tangential_vel_boundary_indicators = this->get_boundary_velocity_manager().get_tangential_boundary_velocity_indicators();
-      velocity_boundary_indicators.insert(tmp_tangential_vel_boundary_indicators.begin(),
-                                          tmp_tangential_vel_boundary_indicators.end());
+      const auto &tangential_boundary_indicators = this->get_boundary_velocity_manager().get_tangential_boundary_velocity_indicators();
+      velocity_boundary_indicators.insert(tangential_boundary_indicators.begin(),
+                                          tangential_boundary_indicators.end());
 
       // Get the active velocity boundary indicators
-      const std::map<types::boundary_id, std::pair<std::string,std::vector<std::string>>>
-      tmp_active_vel_boundary_indicators = this->get_boundary_velocity_manager().get_active_boundary_velocity_names();
-
-      for (const auto &p : tmp_active_vel_boundary_indicators)
-        velocity_boundary_indicators.insert(p.first);
+      const auto &prescribed_boundary_indicators = this->get_boundary_velocity_manager().get_prescribed_boundary_velocity_indicators();
+      velocity_boundary_indicators.insert(prescribed_boundary_indicators.begin(),
+                                          prescribed_boundary_indicators.end());
 
       // Get the mesh deformation boundary indicators
-      const std::set<types::boundary_id> tmp_mesh_deformation_boundary_indicators = this->get_mesh_deformation_boundary_indicators();
-      for (const auto &p : tmp_mesh_deformation_boundary_indicators)
+      for (const auto &p : this->get_mesh_deformation_boundary_indicators())
         AssertThrow(velocity_boundary_indicators.find(p) == velocity_boundary_indicators.end(),
                     ExcMessage("The free surface mesh deformation plugin cannot be used with the current velocity boundary conditions"));
     }
diff --git a/source/simulator/assembly.cc b/source/simulator/assembly.cc
index bb54b960cb6..289c6f1b5c2 100644
--- a/source/simulator/assembly.cc
+++ b/source/simulator/assembly.cc
@@ -743,7 +743,7 @@ namespace aspect
     // we will update the right-hand side with boundary information in
     // StokesMatrixFreeHandler::correct_stokes_rhs().
     if (!stokes_matrix_free)
-      Assert(rebuild_stokes_matrix || boundary_velocity_manager.get_active_boundary_velocity_conditions().size()==0,
+      Assert(rebuild_stokes_matrix || boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().empty(),
              ExcInternalError("If we have inhomogeneous constraints, we must re-assemble the system matrix."));
 
     system_rhs = 0;
diff --git a/source/simulator/core.cc b/source/simulator/core.cc
index e8173dae562..515478a5f88 100644
--- a/source/simulator/core.cc
+++ b/source/simulator/core.cc
@@ -505,8 +505,8 @@ namespace aspect
 
     std::set<types::boundary_id> open_velocity_boundary_indicators
       = geometry_model->get_used_boundary_indicators();
-    for (const auto &p : boundary_velocity_manager.get_active_boundary_velocity_names())
-      open_velocity_boundary_indicators.erase (p.first);
+    for (const auto p : boundary_velocity_manager.get_prescribed_boundary_velocity_indicators())
+      open_velocity_boundary_indicators.erase (p);
     for (const auto p : boundary_velocity_manager.get_zero_boundary_velocity_indicators())
       open_velocity_boundary_indicators.erase (p);
     for (const auto p : boundary_velocity_manager.get_tangential_boundary_velocity_indicators())
@@ -1375,67 +1375,54 @@ namespace aspect
     // set the current time and do the interpolation
     // for the prescribed velocity fields
     boundary_velocity_manager.update();
-    for (const auto &p : boundary_velocity_manager.get_active_boundary_velocity_names())
+
+    // 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())
       {
         Utilities::VectorFunctionFromVelocityFunctionObject<dim> vel
         (introspection.n_components,
          [&] (const dealii::Point<dim> &x) -> Tensor<1,dim>
         {
-          return boundary_velocity_manager.boundary_velocity(p.first, x);
-        });
-
-        // here we create a mask for interpolate_boundary_values out of the 'selector'
-        std::vector<bool> mask(introspection.component_masks.velocities.size(), false);
-        const std::string &comp = p.second.first;
+          if (!assemble_newton_stokes_system || (assemble_newton_stokes_system && nonlinear_iteration == 0))
+            return boundary_velocity_manager.boundary_velocity(boundary_id, x);
+          else
+            return Tensor<1,dim>();
 
-        if (comp.length()>0)
-          {
-            for (const char direction : comp)
-              {
-                switch (direction)
-                  {
-                    case 'x':
-                      mask[introspection.component_indices.velocities[0]] = true;
-                      break;
-                    case 'y':
-                      mask[introspection.component_indices.velocities[1]] = true;
-                      break;
-                    case 'z':
-                      // we must be in 3d, or 'z' should never have gotten through
-                      Assert (dim==3, ExcInternalError());
-                      if (dim==3)
-                        mask[introspection.component_indices.velocities[dim-1]] = true;
-                      break;
-                    default:
-                      Assert (false, ExcInternalError());
-                  }
-              }
-          }
-        else
-          {
-            // no mask given -- take all velocities
-            for (unsigned int i=0; i<introspection.component_masks.velocities.size(); ++i)
-              mask[i]=introspection.component_masks.velocities[i];
-          }
+          return Tensor<1,dim>();
+        });
 
-        if (!assemble_newton_stokes_system || (assemble_newton_stokes_system && nonlinear_iteration == 0))
-          {
-            VectorTools::interpolate_boundary_values (*mapping,
-                                                      dof_handler,
-                                                      p.first,
-                                                      vel,
-                                                      constraints,
-                                                      ComponentMask(mask));
-          }
-        else
-          {
-            VectorTools::interpolate_boundary_values (*mapping,
-                                                      dof_handler,
-                                                      p.first,
-                                                      Functions::ZeroFunction<dim>(introspection.n_components),
-                                                      constraints,
-                                                      ComponentMask(mask));
-          }
+        VectorTools::interpolate_boundary_values (*mapping,
+                                                  dof_handler,
+                                                  boundary_id,
+                                                  vel,
+                                                  constraints,
+                                                  boundary_component_masks.at(boundary_id));
       }
   }
 
diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc
index 8a5a3742909..2215f11a8e3 100644
--- a/source/simulator/helper_functions.cc
+++ b/source/simulator/helper_functions.cc
@@ -2290,7 +2290,7 @@ namespace aspect
       boundary_velocity_manager.get_zero_boundary_velocity_indicators();
 
     const auto &prescribed_velocity_boundaries =
-      boundary_velocity_manager.get_active_boundary_velocity_conditions();
+      boundary_velocity_manager.get_prescribed_boundary_velocity_indicators();
 
     // Loop over all of the boundary faces, ...
     for (const auto &cell : dof_handler.active_cell_iterators())
@@ -2380,104 +2380,50 @@ namespace aspect
   void
   Simulator<dim>::check_consistency_of_boundary_conditions() const
   {
-    // make sure velocity and traction boundary indicators don't appear in multiple lists
-    std::set<types::boundary_id> boundary_indicator_lists[6]
-      = { boundary_velocity_manager.get_zero_boundary_velocity_indicators(),
-          boundary_velocity_manager.get_tangential_boundary_velocity_indicators(),
-          std::set<types::boundary_id>()   // to be prescribed velocity and traction boundary indicators
-        };
-
-    // sets of the boundary indicators only (no selectors and values)
-    std::set<types::boundary_id> velocity_bi;
-    std::set<types::boundary_id> traction_bi;
-
-    for (const auto &p : boundary_velocity_manager.get_active_boundary_velocity_names())
-      velocity_bi.insert(p.first);
-
-    for (const auto &r : parameters.prescribed_traction_boundary_indicators)
-      traction_bi.insert(r.first);
-
-    // are there any indicators that occur in both the prescribed velocity and traction list?
-    std::set<types::boundary_id> intersection;
-    std::set_intersection (velocity_bi.begin(),
-                           velocity_bi.end(),
-                           traction_bi.begin(),
-                           traction_bi.end(),
-                           std::inserter(intersection, intersection.end()));
-
-    // if so, do they have different selectors?
-    if (!intersection.empty())
-      {
-        for (const auto it : intersection)
-          {
-            const std::map<types::boundary_id, std::pair<std::string,std::vector<std::string>>>::const_iterator
-            boundary_velocity_names = boundary_velocity_manager.get_active_boundary_velocity_names().find(it);
-            Assert(boundary_velocity_names != boundary_velocity_manager.get_active_boundary_velocity_names().end(),
-                   ExcInternalError());
-
-            std::set<char> velocity_selector;
-            std::set<char> traction_selector;
-
-            for (const auto it_selector : boundary_velocity_names->second.first)
-              velocity_selector.insert(it_selector);
-
-            for (std::string::const_iterator
-                 it_selector  = parameters.prescribed_traction_boundary_indicators.find(it)->second.first.begin();
-                 it_selector != parameters.prescribed_traction_boundary_indicators.find(it)->second.first.end();
-                 ++it_selector)
-              traction_selector.insert(*it_selector);
-
-            // if there are no selectors specified, throw exception
-            AssertThrow(!velocity_selector.empty() || !traction_selector.empty(),
-                        ExcMessage ("Boundary indicator <"
-                                    +
-                                    Utilities::int_to_string(it)
-                                    +
-                                    "> with symbolic name <"
-                                    +
-                                    geometry_model->translate_id_to_symbol_name (it)
-                                    +
-                                    "> is listed as having both "
-                                    "velocity and traction boundary conditions in the input file."));
-
-            std::set<char> intersection_selector;
-            std::set_intersection (velocity_selector.begin(),
-                                   velocity_selector.end(),
-                                   traction_selector.begin(),
-                                   traction_selector.end(),
-                                   std::inserter(intersection_selector, intersection_selector.end()));
-
-            // if the same selectors are specified, throw exception
-            AssertThrow(intersection_selector.empty(),
-                        ExcMessage ("Selectors of boundary indicator <"
-                                    +
-                                    Utilities::int_to_string(it)
-                                    +
-                                    "> with symbolic name <"
-                                    +
-                                    geometry_model->translate_id_to_symbol_name (it)
-                                    +
-                                    "> are listed as having both "
-                                    "velocity and traction boundary conditions in the input file."));
-          }
-      }
-
-    // remove correct boundary indicators that occur in both the velocity and the traction set
-    // but have different selectors
-    std::set<types::boundary_id> union_set;
-    std::set_union (velocity_bi.begin(),
-                    velocity_bi.end(),
-                    traction_bi.begin(),
-                    traction_bi.end(),
-                    std::inserter(union_set, union_set.end()));
-
-    // assign the prescribed boundary indicator list to the boundary_indicator_lists
-    boundary_indicator_lists[3] = union_set;
+    // a container for the indicators of all boundary conditions
+    std::vector<std::set<types::boundary_id>> boundary_indicator_lists;
+    boundary_indicator_lists.emplace_back(boundary_velocity_manager.get_zero_boundary_velocity_indicators());
+    boundary_indicator_lists.emplace_back(boundary_velocity_manager.get_tangential_boundary_velocity_indicators());
+    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
+    // 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[i]);
+            }
+        }
 
-    // for each combination of boundary indicator lists, make sure that the
+    // for each combination of velocity boundary indicator lists, make sure that the
     // intersection is empty
-    for (unsigned int i=0; i<sizeof(boundary_indicator_lists)/sizeof(boundary_indicator_lists[0]); ++i)
-      for (unsigned int j=i+1; j<sizeof(boundary_indicator_lists)/sizeof(boundary_indicator_lists[0]); ++j)
+    for (unsigned int i=0; i<boundary_indicator_lists.size(); ++i)
+      for (unsigned int j=i+1; j<boundary_indicator_lists.size(); ++j)
         {
           std::set<types::boundary_id> intersection;
           std::set_intersection (boundary_indicator_lists[i].begin(),
@@ -2502,20 +2448,28 @@ namespace aspect
 
     // make sure temperature and heat flux boundary indicators don't appear in multiple lists
     // this is easier than for the velocity/traction, as there are no selectors
-    std::set<types::boundary_id> temperature_bi = boundary_temperature_manager.get_fixed_temperature_boundary_indicators();
-    std::set<types::boundary_id> heat_flux_bi = parameters.fixed_heat_flux_boundary_indicators;
+    boundary_indicator_lists.emplace_back(boundary_temperature_manager.get_fixed_temperature_boundary_indicators());
+    boundary_indicator_lists.emplace_back(parameters.fixed_heat_flux_boundary_indicators);
 
     // are there any indicators that occur in both the prescribed temperature and heat flux list?
     std::set<types::boundary_id> T_intersection;
-    std::set_intersection (temperature_bi.begin(),
-                           temperature_bi.end(),
-                           heat_flux_bi.begin(),
-                           heat_flux_bi.end(),
+    std::set_intersection (boundary_temperature_manager.get_fixed_temperature_boundary_indicators().begin(),
+                           boundary_temperature_manager.get_fixed_temperature_boundary_indicators().end(),
+                           parameters.fixed_heat_flux_boundary_indicators.begin(),
+                           parameters.fixed_heat_flux_boundary_indicators.end(),
                            std::inserter(T_intersection, T_intersection.end()));
 
-    AssertThrow(T_intersection.empty(),
-                ExcMessage ("There is a boundary indicator listed as having both "
-                            "temperature and heat flux boundary conditions in the input file."));
+    AssertThrow (T_intersection.empty(),
+                 ExcMessage ("Boundary indicator <"
+                             +
+                             Utilities::int_to_string(*T_intersection.begin())
+                             +
+                             "> with symbolic name <"
+                             +
+                             geometry_model->translate_id_to_symbol_name (*T_intersection.begin())
+                             +
+                             "> is listed as having more "
+                             "than one type of temperature or heat flux boundary condition in the input file."));
 
     // Check that the periodic boundaries do not have other boundary conditions set
     using periodic_boundary_set
@@ -2524,56 +2478,55 @@ namespace aspect
     periodic_boundary_set pbs = geometry_model->get_periodic_boundary_pairs();
 
     for (const auto &pb : pbs)
-      {
-        // Throw error if we are trying to use the same boundary for more than one boundary condition
-        AssertThrow( is_element( pb.first.first, boundary_temperature_manager.get_fixed_temperature_boundary_indicators() ) == false &&
-                     is_element( pb.first.second, boundary_temperature_manager.get_fixed_temperature_boundary_indicators() ) == false &&
-                     is_element( pb.first.first, boundary_composition_manager.get_fixed_composition_boundary_indicators() ) == false &&
-                     is_element( pb.first.second, boundary_composition_manager.get_fixed_composition_boundary_indicators() ) == false &&
-                     is_element( pb.first.first, boundary_indicator_lists[0] ) == false && // zero velocity
-                     is_element( pb.first.second, boundary_indicator_lists[0] ) == false && // zero velocity
-                     is_element( pb.first.first, boundary_indicator_lists[1] ) == false && // tangential velocity
-                     is_element( pb.first.second, boundary_indicator_lists[1] ) == false && // tangential velocity
-                     is_element( pb.first.first, boundary_indicator_lists[3] ) == false && // prescribed traction or velocity
-                     is_element( pb.first.second, boundary_indicator_lists[3] ) == false,  // prescribed traction or velocity
-                     ExcMessage("Periodic boundaries must not have boundary conditions set."));
-      }
+      for (const auto &boundary_indicators: boundary_indicator_lists)
+        {
+          AssertThrow(is_element(pb.first.first, boundary_indicators) == false,
+                      ExcMessage ("Boundary indicator <"
+                                  +
+                                  Utilities::int_to_string(pb.first.first)
+                                  +
+                                  "> with symbolic name <"
+                                  +
+                                  geometry_model->translate_id_to_symbol_name (pb.first.first)
+                                  +
+                                  "> is listed as having a periodic boundary condition "
+                                  "in the input file, but also has another type of boundary condition. "
+                                  "Periodic boundaries cannot have other boundary conditions."));
+
+          AssertThrow(is_element(pb.first.second, boundary_indicators) == false,
+                      ExcMessage ("Boundary indicator <"
+                                  +
+                                  Utilities::int_to_string(pb.first.second)
+                                  +
+                                  "> with symbolic name <"
+                                  +
+                                  geometry_model->translate_id_to_symbol_name (pb.first.second)
+                                  +
+                                  "> is listed as having a periodic boundary condition "
+                                  "in the input file, but also has another type of boundary condition. "
+                                  "Periodic boundaries cannot have other boundary conditions."));
+        }
 
     const std::set<types::boundary_id> all_boundary_indicators
       = geometry_model->get_used_boundary_indicators();
-    if (parameters.nonlinear_solver != NonlinearSolver::single_Advection_no_Stokes)
-      {
-        // next make sure that all listed indicators are actually used by
-        // this geometry
-        for (const auto &list : boundary_indicator_lists)
-          for (const auto &p : list)
-            AssertThrow (all_boundary_indicators.find (p)
-                         != all_boundary_indicators.end(),
-                         ExcMessage ("One of the boundary indicators listed in the input file "
-                                     "is not used by the geometry model."));
-      }
-    else
+
+    // next make sure that all listed indicators are actually used by
+    // this geometry
+    for (const auto &list : boundary_indicator_lists)
+      for (const auto &p : list)
+        AssertThrow (all_boundary_indicators.find (p)
+                     != all_boundary_indicators.end(),
+                     ExcMessage ("One of the boundary indicators listed in the input file "
+                                 "is not used by the geometry model."));
+
+    if (parameters.nonlinear_solver == NonlinearSolver::single_Advection_no_Stokes)
       {
-        // next make sure that there are no listed indicators
-        for (const auto &list : boundary_indicator_lists)
-          AssertThrow (list.empty(),
+        // make sure that there are no listed velocity boundary conditions
+        for (unsigned int i=0; i<4; ++i)
+          AssertThrow (boundary_indicator_lists[i].empty(),
                        ExcMessage ("With the solver scheme `single Advection, no Stokes', "
                                    "one cannot set boundary conditions for velocity."));
       }
-
-
-    // now do the same for the fixed temperature indicators and the
-    // compositional indicators
-    for (const auto p : boundary_temperature_manager.get_fixed_temperature_boundary_indicators())
-      AssertThrow (all_boundary_indicators.find (p)
-                   != all_boundary_indicators.end(),
-                   ExcMessage ("One of the fixed boundary temperature indicators listed in the input file "
-                               "is not used by the geometry model."));
-    for (const auto p : boundary_composition_manager.get_fixed_composition_boundary_indicators())
-      AssertThrow (all_boundary_indicators.find (p)
-                   != all_boundary_indicators.end(),
-                   ExcMessage ("One of the fixed boundary composition indicators listed in the input file "
-                               "is not used by the geometry model."));
   }
 
 
@@ -2598,8 +2551,8 @@ namespace aspect
     // rebuild the whole system to compute the rhs.
     assemble_newton_stokes_system = true;
     rebuild_stokes_preconditioner = false;
-    rebuild_stokes_matrix = boundary_velocity_manager.get_active_boundary_velocity_conditions().size()!=0;
-    assemble_newton_stokes_matrix = boundary_velocity_manager.get_active_boundary_velocity_conditions().size()!=0;
+    rebuild_stokes_matrix = !boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().empty();
+    assemble_newton_stokes_matrix = !boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().empty();
 
     compute_current_constraints ();
 
diff --git a/source/simulator/solver_schemes.cc b/source/simulator/solver_schemes.cc
index 4c4f71ba4a9..e21aaa61016 100644
--- a/source/simulator/solver_schemes.cc
+++ b/source/simulator/solver_schemes.cc
@@ -370,7 +370,7 @@ namespace aspect
     // don't need to force assembly of the matrix.
     if (stokes_matrix_depends_on_solution()
         ||
-        (boundary_velocity_manager.get_active_boundary_velocity_conditions().size() > 0)
+        (boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().size() > 0)
         || parameters.mesh_deformation_enabled
        )
       rebuild_stokes_matrix = rebuild_stokes_preconditioner = true;
@@ -579,7 +579,7 @@ namespace aspect
             // don't need to force assembly of the matrix.
             if (stokes_matrix_depends_on_solution()
                 ||
-                (nonlinear_iteration == 0 && boundary_velocity_manager.get_active_boundary_velocity_conditions().size() > 0))
+                (nonlinear_iteration == 0 && boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().size() > 0))
               rebuild_stokes_matrix = rebuild_stokes_preconditioner = assemble_newton_stokes_matrix = true;
             else if (parameters.enable_prescribed_dilation)
               // The dilation requires the Stokes matrix (which is on the rhs
@@ -680,8 +680,7 @@ namespace aspect
 
             // Rebuild the rhs to determine the new residual.
             assemble_newton_stokes_matrix = rebuild_stokes_preconditioner = false;
-            rebuild_stokes_matrix = (boundary_velocity_manager.get_active_boundary_velocity_conditions().empty()
-                                     == false);
+            rebuild_stokes_matrix = !boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().empty();
 
             assemble_stokes_system();
 
diff --git a/source/simulator/stokes_matrix_free.cc b/source/simulator/stokes_matrix_free.cc
index f867d462cb4..5a716de77ce 100644
--- a/source/simulator/stokes_matrix_free.cc
+++ b/source/simulator/stokes_matrix_free.cc
@@ -2487,35 +2487,19 @@ 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 (const auto &it: sim.boundary_velocity_manager.get_active_boundary_velocity_names())
+      for (unsigned int i=0; i<sim.boundary_velocity_manager.get_active_plugins().size(); ++i)
         {
-          const types::boundary_id bdryid = it.first;
-          const std::string component=it.second.first;
+          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];
 
-          if (component.length()>0)
+          if (component_mask != ComponentMask(sim.introspection.n_components, false))
             {
-              std::vector<bool> mask(fe_v.n_components(), false);
-              for (const auto &direction : component)
-                {
-                  switch (direction)
-                    {
-                      case 'x':
-                        mask[0] = true;
-                        break;
-                      case 'y':
-                        mask[1] = true;
-                        break;
-                      case 'z':
-                        // we must be in 3d, or 'z' should never have gotten through
-                        Assert (dim==3, ExcInternalError());
-                        if (dim==3)
-                          mask[2] = true;
-                        break;
-                      default:
-                        Assert (false, ExcInternalError());
-                    }
-                }
-              mg_constrained_dofs_A_block.make_zero_boundary_constraints(dof_handler_v, {bdryid}, ComponentMask(mask));
+              ComponentMask 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]]);
+
+              mg_constrained_dofs_A_block.make_zero_boundary_constraints(dof_handler_v, {bdryid}, mask);
             }
           else
             {
diff --git a/tests/ascii_boundary_member.cc b/tests/ascii_boundary_member.cc
index 16b580f44c3..447ccd745c3 100644
--- a/tests/ascii_boundary_member.cc
+++ b/tests/ascii_boundary_member.cc
@@ -98,16 +98,17 @@ namespace aspect
   void
   AsciiBoundaryMember<dim>::initialize ()
   {
-    const std::map<types::boundary_id,std::vector<std::unique_ptr<BoundaryVelocity::Interface<dim>>>> &
-    bvs = this->get_boundary_velocity_manager().get_active_boundary_velocity_conditions();
-    for (const auto &boundary : bvs)
-      for (const auto &p : boundary.second)
-        {
-          if (p.get() == this)
-            boundary_ids.insert(boundary.first);
-        }
-    AssertThrow(*(boundary_ids.begin()) != numbers::invalid_boundary_id,
-                ExcMessage("Did not find the boundary indicator for the prescribed data plugin."));
+    unsigned int i=0;
+    for (const auto &plugin : this->get_boundary_velocity_manager().get_active_plugins())
+      {
+        if (plugin.get() == this)
+          boundary_ids.insert(this->get_boundary_velocity_manager().get_active_plugin_boundary_indicators()[i]);
+
+        ++i;
+      }
+
+    AssertThrow(boundary_ids.empty() == false,
+                ExcMessage("Did not find the boundary indicator for the velocity ascii data plugin."));
 
     member->initialize(boundary_ids,
                        dim);