[PL] Pass bc mesh to boundary conditions.

ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp

@@ -39,36 +39,35 @@ std::unique_ptr<BoundaryCondition> createBoundaryCondition(
     if (type == "Dirichlet")
     {
         return ProcessLib::createDirichletBoundaryCondition(
-            config.config, config.boundary_mesh, dof_table, bulk_mesh.getID(),
-            variable_id, *config.component_id, parameters);
+            config.config, config.boundary_mesh, dof_table, variable_id,
+            *config.component_id, parameters);
     }
     if (type == "Neumann")
     {
         return ProcessLib::createNeumannBoundaryCondition(
             config.config, config.boundary_mesh, dof_table, variable_id,
-            *config.component_id, bulk_mesh.isAxiallySymmetric(),
-            integration_order, shapefunction_order, bulk_mesh.getDimension(),
-            parameters);
+            *config.component_id, integration_order, shapefunction_order,
+            bulk_mesh.getDimension(), parameters);
     }
     if (type == "Robin")
     {
         return ProcessLib::createRobinBoundaryCondition(
             config.config, config.boundary_mesh, dof_table, variable_id,
-            *config.component_id, bulk_mesh.isAxiallySymmetric(),
-            integration_order, shapefunction_order, bulk_mesh.getDimension(),
-            parameters);
+            *config.component_id, integration_order, shapefunction_order,
+            bulk_mesh.getDimension(), parameters);
     }
     if (type == "NonuniformDirichlet")
     {
         return ProcessLib::createNonuniformDirichletBoundaryCondition(
-            config.config, dof_table, variable_id, *config.component_id,
-            bulk_mesh);
+            config.config, config.boundary_mesh, dof_table, variable_id,
+            *config.component_id, bulk_mesh);
     }
     if (type == "NonuniformNeumann")
     {
         return ProcessLib::createNonuniformNeumannBoundaryCondition(
-            config.config, dof_table, variable_id, *config.component_id,
-            integration_order, shapefunction_order, bulk_mesh);
+            config.config, config.boundary_mesh, dof_table, variable_id,
+            *config.component_id, integration_order, shapefunction_order,
+            bulk_mesh);
     }
     if (type == "Python")
     {
@@ -96,8 +95,8 @@ std::unique_ptr<BoundaryCondition> createBoundaryCondition(
         return ProcessLib::NormalTractionBoundaryCondition::
             createNormalTractionBoundaryCondition(
                 config.config, config.boundary_mesh, dof_table, variable_id,
-                bulk_mesh.isAxiallySymmetric(), integration_order,
-                shapefunction_order, bulk_mesh.getDimension(), parameters);
+                integration_order, shapefunction_order,
+                bulk_mesh.getDimension(), parameters);
     }
     if (type == "PhaseFieldIrreversibleDamageOracleBoundaryCondition")
     {

ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp

@@ -22,10 +22,6 @@ void DirichletBoundaryCondition::getEssentialBCValues(
 {
     SpatialPosition pos;
 
-    auto const& bulk_node_ids_map =
-        *_bc_mesh.getProperties().getPropertyVector<std::size_t>(
-            "bulk_node_ids");
-
     bc_values.ids.clear();
     bc_values.values.clear();
 
@@ -35,23 +31,24 @@ void DirichletBoundaryCondition::getEssentialBCValues(
                              _bc_mesh.getNumberOfNodes());
     for (auto const* const node : _bc_mesh.getNodes())
     {
-        auto const id = bulk_node_ids_map[node->getID()];
-        pos.setNodeID(id);
-        MeshLib::Location l(_bulk_mesh_id, MeshLib::MeshItemType::Node, id);
+        auto const id = node->getID();
+        pos.setNodeID(node->getID());
         // TODO: that might be slow, but only done once
-        const auto g_idx =
-            _dof_table.getGlobalIndex(l, _variable_id, _component_id);
-        if (g_idx == NumLib::MeshComponentMap::nop)
+        auto const global_index = _dof_table_boundary->getGlobalIndex(
+            {_bc_mesh.getID(), MeshLib::MeshItemType::Node, id}, _variable_id,
+            _component_id);
+        if (global_index == NumLib::MeshComponentMap::nop)
             continue;
         // For the DDC approach (e.g. with PETSc option), the negative
-        // index of g_idx means that the entry by that index is a ghost one,
-        // which should be dropped. Especially for PETSc routines MatZeroRows
-        // and MatZeroRowsColumns, which are called to apply the Dirichlet BC,
-        // the negative index is not accepted like other matrix or vector
-        // PETSc routines. Therefore, the following if-condition is applied.
-        if (g_idx >= 0)
+        // index of global_index means that the entry by that index is a ghost
+        // one, which should be dropped. Especially for PETSc routines
+        // MatZeroRows and MatZeroRowsColumns, which are called to apply the
+        // Dirichlet BC, the negative index is not accepted like other matrix or
+        // vector PETSc routines. Therefore, the following if-condition is
+        // applied.
+        if (global_index >= 0)
         {
-            bc_values.ids.emplace_back(g_idx);
+            bc_values.ids.emplace_back(global_index);
             bc_values.values.emplace_back(_parameter(t, pos).front());
         }
     }
@@ -59,8 +56,7 @@ void DirichletBoundaryCondition::getEssentialBCValues(
 
 std::unique_ptr<DirichletBoundaryCondition> createDirichletBoundaryCondition(
     BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
-    NumLib::LocalToGlobalIndexMap const& dof_table,
-    std::size_t const bulk_mesh_id, int const variable_id,
+    NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
     int const component_id,
     const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>& parameters)
 {
@@ -75,7 +71,7 @@ std::unique_ptr<DirichletBoundaryCondition> createDirichletBoundaryCondition(
     auto& param = findParameter<double>(param_name, parameters, 1);
 
     return std::make_unique<DirichletBoundaryCondition>(
-        param, bc_mesh, dof_table, bulk_mesh_id, variable_id, component_id);
+        param, bc_mesh, dof_table_bulk, variable_id, component_id);
 }
 
 }  // namespace ProcessLib

ProcessLib/BoundaryCondition/DirichletBoundaryCondition.h

@@ -24,28 +24,26 @@ namespace ProcessLib
 class DirichletBoundaryCondition final : public BoundaryCondition
 {
 public:
-    DirichletBoundaryCondition(Parameter<double> const& parameter,
-                               MeshLib::Mesh const& bc_mesh,
-                               NumLib::LocalToGlobalIndexMap const& dof_table,
-                               std::size_t const bulk_mesh_id,
-                               int const variable_id, int const component_id)
+    DirichletBoundaryCondition(
+        Parameter<double> const& parameter, MeshLib::Mesh const& bc_mesh,
+        NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
+        int const variable_id, int const component_id)
         : _parameter(parameter),
           _bc_mesh(bc_mesh),
-          _dof_table(dof_table),
-          _bulk_mesh_id(bulk_mesh_id),
           _variable_id(variable_id),
           _component_id(component_id)
     {
-        if (variable_id >= static_cast<int>(dof_table.getNumberOfVariables()) ||
+        if (variable_id >=
+                static_cast<int>(dof_table_bulk.getNumberOfVariables()) ||
             component_id >=
-                dof_table.getNumberOfVariableComponents(variable_id))
+                dof_table_bulk.getNumberOfVariableComponents(variable_id))
         {
             OGS_FATAL(
                 "Variable id or component id too high. Actual values: (%d, "
-                "%d), "
-                "maximum values: (%d, %d).",
-                variable_id, component_id, dof_table.getNumberOfVariables(),
-                dof_table.getNumberOfVariableComponents(variable_id));
+                "%d), maximum values: (%d, %d).",
+                variable_id, component_id,
+                dof_table_bulk.getNumberOfVariables(),
+                dof_table_bulk.getNumberOfVariableComponents(variable_id));
         }
 
         if (!_bc_mesh.getProperties().existsPropertyVector<std::size_t>(
@@ -55,6 +53,20 @@ public:
                 "The required bulk node ids map does not exist in the boundary "
                 "mesh '%s'.", _bc_mesh.getName().c_str());
         }
+
+        std::vector<MeshLib::Node*> const& bc_nodes = _bc_mesh.getNodes();
+        DBUG(
+            "Found %d nodes for Dirichlet BCs for the variable %d and "
+            "component "
+            "%d",
+            bc_nodes.size(), variable_id, component_id);
+
+        MeshLib::MeshSubset bc_mesh_subset(_bc_mesh, bc_nodes);
+
+        // Create local DOF table from the bc mesh subset for the given variable
+        // and component id.
+        _dof_table_boundary.reset(dof_table_bulk.deriveBoundaryConstrainedMap(
+            variable_id, {component_id}, std::move(bc_mesh_subset)));
     }
 
     void getEssentialBCValues(
@@ -65,16 +77,15 @@ private:
     Parameter<double> const& _parameter;
 
     MeshLib::Mesh const& _bc_mesh;
-    NumLib::LocalToGlobalIndexMap const& _dof_table;
-    std::size_t const _bulk_mesh_id;
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap const> _dof_table_boundary;
     int const _variable_id;
     int const _component_id;
 };
 
 std::unique_ptr<DirichletBoundaryCondition> createDirichletBoundaryCondition(
     BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
-    NumLib::LocalToGlobalIndexMap const& dof_table, std::size_t const mesh_id,
-    int const variable_id, int const component_id,
+    NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
+    int const component_id,
     const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>& parameters);
 
 }  // namespace ProcessLib

ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h

@@ -7,9 +7,8 @@
  *
  */
 
-#include "GenericNaturalBoundaryCondition.h"
-#include "ProcessLib/Utils/CreateLocalAssemblers.h"
 #include "GenericNaturalBoundaryConditionLocalAssembler.h"
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
 
 namespace ProcessLib
 {
@@ -20,18 +19,16 @@ template <typename Data>
 GenericNaturalBoundaryCondition<BoundaryConditionData,
                                 LocalAssemblerImplementation>::
     GenericNaturalBoundaryCondition(
-        typename std::enable_if<
-            std::is_same<typename std::decay<BoundaryConditionData>::type,
-                         typename std::decay<Data>::type>::value,
-            bool>::type is_axially_symmetric,
         unsigned const integration_order, unsigned const shapefunction_order,
         NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
         int const variable_id, int const component_id,
         unsigned const global_dim, MeshLib::Mesh const& bc_mesh, Data&& data)
-    : _data(std::forward<Data>(data)),
-      _bc_mesh(bc_mesh),
-      _integration_order(integration_order)
+    : _data(std::forward<Data>(data)), _bc_mesh(bc_mesh)
 {
+    static_assert(std::is_same<typename std::decay<BoundaryConditionData>::type,
+                               typename std::decay<Data>::type>::value,
+                  "Type mismatch between declared and passed BC data.");
+
     // check basic data consistency
     if (variable_id >=
             static_cast<int>(dof_table_bulk.getNumberOfVariables()) ||
@@ -45,6 +42,23 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
             dof_table_bulk.getNumberOfVariableComponents(variable_id));
     }
 
+    if (_bc_mesh.getDimension() + 1 != global_dim)
+    {
+        OGS_FATAL(
+            "The dimension of the given boundary mesh (%d) is not by one lower "
+            "than the bulk dimension (%d).",
+            _bc_mesh.getDimension(), global_dim);
+    }
+
+    if (!_bc_mesh.getProperties().template existsPropertyVector<std::size_t>(
+            "bulk_node_ids"))
+    {
+        OGS_FATAL(
+            "The required bulk node ids map does not exist in the boundary "
+            "mesh '%s'.",
+            _bc_mesh.getName().c_str());
+    }
+
     std::vector<MeshLib::Node*> const& bc_nodes = _bc_mesh.getNodes();
     DBUG("Found %d nodes for Natural BCs for the variable %d and component %d",
          bc_nodes.size(), variable_id, component_id);
@@ -58,8 +72,8 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
 
     createLocalAssemblers<LocalAssemblerImplementation>(
         global_dim, _bc_mesh.getElements(), *_dof_table_boundary,
-        shapefunction_order, _local_assemblers, is_axially_symmetric,
-        _integration_order, _data);
+        shapefunction_order, _local_assemblers, _bc_mesh.isAxiallySymmetric(),
+        integration_order, _data);
 }
 
 template <typename BoundaryConditionData,

ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h

@@ -27,10 +27,6 @@ public:
     /// A local DOF-table, a subset of the given one, is constructed.
     template <typename Data>
     GenericNaturalBoundaryCondition(
-        typename std::enable_if<
-            std::is_same<typename std::decay<BoundaryConditionData>::type,
-                         typename std::decay<Data>::type>::value,
-            bool>::type is_axially_symmetric,
         unsigned const integration_order, unsigned const shapefunction_order,
         NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
         int const variable_id, int const component_id,
@@ -52,9 +48,6 @@ private:
     /// participating number of _elements of the boundary condition.
     std::unique_ptr<NumLib::LocalToGlobalIndexMap> _dof_table_boundary;
 
-    /// Integration order for integration over the lower-dimensional elements
-    unsigned const _integration_order;
-
     /// Local assemblers for each element of number of _elements.
     std::vector<
         std::unique_ptr<GenericNaturalBoundaryConditionLocalAssemblerInterface>>

ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h

@@ -7,9 +7,7 @@
  *
  */
 
-#include "GenericNonuniformNaturalBoundaryCondition.h"
 #include "GenericNonuniformNaturalBoundaryConditionLocalAssembler.h"
-#include "MeshLib/MeshSearch/NodeSearch.h"
 #include "ProcessLib/Utils/CreateLocalAssemblers.h"
 
 namespace ProcessLib
@@ -22,63 +20,51 @@ GenericNonuniformNaturalBoundaryCondition<BoundaryConditionData,
                                           LocalAssemblerImplementation>::
     GenericNonuniformNaturalBoundaryCondition(
         unsigned const integration_order, unsigned const shapefunction_order,
-        unsigned const global_dim,
-        std::unique_ptr<MeshLib::Mesh>&& boundary_mesh, Data&& data)
-    : _data(std::forward<Data>(data)), _boundary_mesh(std::move(boundary_mesh))
+        NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
+        int const variable_id, int const component_id,
+        unsigned const global_dim, MeshLib::Mesh const& bc_mesh, Data&& data)
+    : _data(std::forward<Data>(data)), _bc_mesh(bc_mesh)
 {
     static_assert(std::is_same<typename std::decay<BoundaryConditionData>::type,
                                typename std::decay<Data>::type>::value,
                   "Type mismatch between declared and passed BC data.");
 
     // check basic data consistency
-    if (_data.variable_id_bulk >=
-            static_cast<int>(_data.dof_table_bulk.getNumberOfVariables()) ||
-        _data.component_id_bulk >=
-            _data.dof_table_bulk.getNumberOfVariableComponents(
-                _data.variable_id_bulk))
+    if (variable_id >=
+            static_cast<int>(dof_table_bulk.getNumberOfVariables()) ||
+        component_id >=
+            dof_table_bulk.getNumberOfVariableComponents(variable_id))
     {
         OGS_FATAL(
             "Variable id or component id too high. Actual values: (%d, %d), "
             "maximum values: (%d, %d).",
-            _data.variable_id_bulk, _data.component_id_bulk,
-            _data.dof_table_bulk.getNumberOfVariables(),
-            _data.dof_table_bulk.getNumberOfVariableComponents(
-                _data.variable_id_bulk));
+            variable_id, component_id, dof_table_bulk.getNumberOfVariables(),
+            dof_table_bulk.getNumberOfVariableComponents(variable_id));
     }
 
-    if (_boundary_mesh->getDimension() + 1 != global_dim)
+    if (_bc_mesh.getDimension() + 1 != global_dim)
     {
         OGS_FATAL(
             "The dimension of the given boundary mesh (%d) is not by one lower "
             "than the bulk dimension (%d).",
-            _boundary_mesh->getDimension(), global_dim);
+            _bc_mesh.getDimension(), global_dim);
     }
 
-    constructDofTable();
+    std::vector<MeshLib::Node*> const& bc_nodes = _bc_mesh.getNodes();
+    DBUG("Found %d nodes for Natural BCs for the variable %d and component %d",
+         bc_nodes.size(), variable_id, component_id);
 
-    createLocalAssemblers<LocalAssemblerImplementation>(
-        global_dim, _boundary_mesh->getElements(), *_dof_table_boundary,
-        shapefunction_order, _local_assemblers,
-        _boundary_mesh->isAxiallySymmetric(), integration_order, _data);
-}
-
-template <typename BoundaryConditionData,
-          template <typename, typename, unsigned>
-          class LocalAssemblerImplementation>
-void GenericNonuniformNaturalBoundaryCondition<
-    BoundaryConditionData, LocalAssemblerImplementation>::constructDofTable()
-{
-    // construct one-component DOF-table for the surface mesh
-    _mesh_subset_all_nodes.reset(
-        new MeshLib::MeshSubset(*_boundary_mesh, _boundary_mesh->getNodes()));
+    MeshLib::MeshSubset bc_mesh_subset(_bc_mesh, bc_nodes);
 
-    std::vector<MeshLib::MeshSubset> all_mesh_subsets{*_mesh_subset_all_nodes};
+    // Create local DOF table from the bc mesh subset for the given variable and
+    // component id.
+    _dof_table_boundary.reset(dof_table_bulk.deriveBoundaryConstrainedMap(
+        variable_id, {component_id}, std::move(bc_mesh_subset)));
 
-    std::vector<int> vec_var_n_components{1};
-
-    _dof_table_boundary = std::make_unique<NumLib::LocalToGlobalIndexMap>(
-        std::move(all_mesh_subsets), vec_var_n_components,
-        NumLib::ComponentOrder::BY_LOCATION);
+    createLocalAssemblers<LocalAssemblerImplementation>(
+        global_dim, _bc_mesh.getElements(), *_dof_table_boundary,
+        shapefunction_order, _local_assemblers, _bc_mesh.isAxiallySymmetric(),
+        integration_order, _data);
 }
 
 template <typename BoundaryConditionData,

ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h

@@ -28,8 +28,9 @@ public:
     template <typename Data>
     GenericNonuniformNaturalBoundaryCondition(
         unsigned const integration_order, unsigned const shapefunction_order,
-        unsigned const global_dim,
-        std::unique_ptr<MeshLib::Mesh>&& boundary_mesh, Data&& data);
+        NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
+        int const variable_id, int const component_id,
+        unsigned const global_dim, MeshLib::Mesh const& bc_mesh, Data&& data);
 
     /// Calls local assemblers which calculate their contributions to the global
     /// matrix and the right-hand-side.
@@ -37,16 +38,14 @@ public:
                         GlobalVector& b, GlobalMatrix* Jac) override;
 
 private:
-    void constructDofTable();
-
     /// Data used in the assembly of the specific boundary condition.
     BoundaryConditionData _data;
 
-    std::unique_ptr<MeshLib::Mesh> _boundary_mesh;
-
-    std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_all_nodes;
+    /// A lower-dimensional mesh on which the boundary condition is defined.
+    MeshLib::Mesh const& _bc_mesh;
 
-    /// DOF-table (one single-component variable) of the boundary mesh.
+    /// Local dof table, a subset of the global one restricted to the
+    /// participating number of _elements of the boundary condition.
     std::unique_ptr<NumLib::LocalToGlobalIndexMap> _dof_table_boundary;
 
     /// Local assemblers for each element of the boundary mesh.

ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryConditionLocalAssembler.h

@@ -79,12 +79,16 @@ public:
     GenericNonuniformNaturalBoundaryConditionLocalAssembler(
         MeshLib::Element const& e, bool is_axially_symmetric,
         unsigned const integration_order)
-        : _ns_and_weights(
+        : _integration_method(integration_order),
+          _element(e),
+          _ns_and_weights(
               initNsAndWeights(e, is_axially_symmetric, integration_order))
     {
     }
 
 protected:
+    IntegrationMethod const _integration_method;
+    MeshLib::Element const& _element;
     std::vector<NAndWeight, Eigen::aligned_allocator<NAndWeight>> const
         _ns_and_weights;
 };

ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp

@@ -15,9 +15,8 @@ namespace ProcessLib
 std::unique_ptr<NeumannBoundaryCondition> createNeumannBoundaryCondition(
     BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
     NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
-    int const component_id, bool is_axially_symmetric,
-    unsigned const integration_order, unsigned const shapefunction_order,
-    unsigned const global_dim,
+    int const component_id, unsigned const integration_order,
+    unsigned const shapefunction_order, unsigned const global_dim,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters)
 {
     DBUG("Constructing Neumann BC from config.");
@@ -31,8 +30,8 @@ std::unique_ptr<NeumannBoundaryCondition> createNeumannBoundaryCondition(
     auto const& param = findParameter<double>(param_name, parameters, 1);
 
     return std::make_unique<NeumannBoundaryCondition>(
-        is_axially_symmetric, integration_order, shapefunction_order, dof_table,
-        variable_id, component_id, global_dim, bc_mesh, param);
+        integration_order, shapefunction_order, dof_table, variable_id,
+        component_id, global_dim, bc_mesh, param);
 }
 
 }  // ProcessLib

ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h

@@ -21,9 +21,8 @@ using NeumannBoundaryCondition = GenericNaturalBoundaryCondition<
 std::unique_ptr<NeumannBoundaryCondition> createNeumannBoundaryCondition(
     BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
     NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
-    int const component_id, bool is_axially_symmetric,
-    unsigned const integration_order, unsigned const shapefunction_order,
-    unsigned const global_dim,
+    int const component_id, unsigned const integration_order,
+    unsigned const shapefunction_order, unsigned const global_dim,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters);
 
 }  // ProcessLib

ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.cpp

@@ -18,7 +18,7 @@ namespace ProcessLib
 {
 std::unique_ptr<NonuniformDirichletBoundaryCondition>
 createNonuniformDirichletBoundaryCondition(
-    BaseLib::ConfigTree const& config,
+    BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh,
     NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
     int const component_id, MeshLib::Mesh const& bulk_mesh)
 {
@@ -26,23 +26,18 @@ createNonuniformDirichletBoundaryCondition(
     //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
     config.checkConfigParameter("type", "NonuniformDirichlet");
 
-    // TODO handle paths correctly
-    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__NonuniformDirichlet__mesh}
-    auto const mesh_file = config.getConfigParameter<std::string>("mesh");
-
-    std::unique_ptr<MeshLib::Mesh> boundary_mesh(
-        MeshLib::IO::readMeshFromFile(mesh_file));
-
-    if (!boundary_mesh)
-    {
-        OGS_FATAL("Error reading mesh `%s'", mesh_file.c_str());
-    }
-
     // The axial symmetry is not used in the Dirichlet BC but kept here for
     // consistency.
     //
     // Surface mesh and bulk mesh must have equal axial symmetry flags!
-    boundary_mesh->setAxiallySymmetric(bulk_mesh.isAxiallySymmetric());
+    if (boundary_mesh.isAxiallySymmetric() != bulk_mesh.isAxiallySymmetric())
+    {
+        OGS_FATAL(
+            "The boundary mesh %s axially symmetric but the bulk mesh %s. Both "
+            "must have an equal axial symmetry property.",
+            boundary_mesh.isAxiallySymmetric() ? "is" : "is not",
+            bulk_mesh.isAxiallySymmetric() ? "is" : "is not");
+    }
 
     // TODO finally use ProcessLib::Parameter here
     auto const field_name =
@@ -50,12 +45,12 @@ createNonuniformDirichletBoundaryCondition(
         config.getConfigParameter<std::string>("field_name");
 
     auto const* const property =
-        boundary_mesh->getProperties().getPropertyVector<double>(field_name);
+        boundary_mesh.getProperties().getPropertyVector<double>(field_name);
 
     if (!property)
     {
         OGS_FATAL("A property with name `%s' does not exist in `%s'.",
-                  field_name.c_str(), mesh_file.c_str());
+                  field_name.c_str(), boundary_mesh.getName().c_str());
     }
 
     if (property->getMeshItemType() != MeshLib::MeshItemType::Node)
@@ -71,24 +66,8 @@ createNonuniformDirichletBoundaryCondition(
         OGS_FATAL("`%s' is not a one-component field.", field_name.c_str());
     }
 
-    std::string const mapping_to_bulk_nodes_property =
-        "OriginalSubsurfaceNodeIDs";
-    auto const* const mapping_to_bulk_nodes =
-        boundary_mesh->getProperties().getPropertyVector<std::size_t>(
-            mapping_to_bulk_nodes_property);
-
-    if (!(mapping_to_bulk_nodes && mapping_to_bulk_nodes->getMeshItemType() ==
-                                       MeshLib::MeshItemType::Node) &&
-        mapping_to_bulk_nodes->getNumberOfComponents() == 1)
-    {
-        OGS_FATAL("Field `%s' is not set up properly.",
-                  mapping_to_bulk_nodes_property.c_str());
-    }
-
     return std::make_unique<NonuniformDirichletBoundaryCondition>(
-        // bulk_mesh.getDimension(),
-        std::move(boundary_mesh), *property, bulk_mesh.getID(),
-        *mapping_to_bulk_nodes, dof_table, variable_id, component_id);
+        boundary_mesh, *property, dof_table, variable_id, component_id);
 }
 
 }  // ProcessLib

ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h

@@ -27,36 +27,51 @@ class NonuniformDirichletBoundaryCondition final : public BoundaryCondition
 public:
     NonuniformDirichletBoundaryCondition(
         // int const bulk_mesh_dimension,
-        std::unique_ptr<MeshLib::Mesh>
-            boundary_mesh,
+        MeshLib::Mesh const& boundary_mesh,
         MeshLib::PropertyVector<double> const& values,
-        std::size_t const bulk_mesh_id,
-        MeshLib::PropertyVector<std::size_t> const& mapping_to_bulk_nodes,
         NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
-        int const variable_id_bulk,
-        int const component_id_bulk)
-        : _bulk_mesh_id(bulk_mesh_id),
-          _values(values),
-          _boundary_mesh(std::move(boundary_mesh)),
-          _mapping_to_bulk_nodes(mapping_to_bulk_nodes),
-          _dof_table_bulk(dof_table_bulk),
-          _variable_id_bulk(variable_id_bulk),
-          _component_id_bulk(component_id_bulk)
+        int const variable_id,
+        int const component_id)
+        : _values(values),
+          _boundary_mesh(boundary_mesh),
+          _variable_id(variable_id),
+          _component_id(component_id)
     {
-        if (_variable_id_bulk >=
-                static_cast<int>(_dof_table_bulk.getNumberOfVariables()) ||
-            _component_id_bulk >= _dof_table_bulk.getNumberOfVariableComponents(
-                                      _variable_id_bulk))
+        if (_variable_id >=
+                static_cast<int>(dof_table_bulk.getNumberOfVariables()) ||
+            _component_id >=
+                dof_table_bulk.getNumberOfVariableComponents(_variable_id))
         {
             OGS_FATAL(
                 "Variable id or component id too high. Actual values: (%d, "
                 "%d), "
                 "maximum values: (%d, %d).",
-                _variable_id_bulk, _component_id_bulk,
-                _dof_table_bulk.getNumberOfVariables(),
-                _dof_table_bulk.getNumberOfVariableComponents(
-                    _variable_id_bulk));
+                _variable_id, _component_id,
+                dof_table_bulk.getNumberOfVariables(),
+                dof_table_bulk.getNumberOfVariableComponents(_variable_id));
         }
+
+        if (!_boundary_mesh.getProperties().existsPropertyVector<std::size_t>(
+                "bulk_node_ids"))
+        {
+            OGS_FATAL(
+                "The required bulk node ids map does not exist in the boundary "
+                "mesh '%s'.",
+                _boundary_mesh.getName().c_str());
+        }
+
+        std::vector<MeshLib::Node*> const& bc_nodes = _boundary_mesh.getNodes();
+        DBUG(
+            "Found %d nodes for Natural BCs for the variable %d and component "
+            "%d",
+            bc_nodes.size(), variable_id, component_id);
+
+        MeshLib::MeshSubset boundary_mesh_subset(_boundary_mesh, bc_nodes);
+
+        // Create local DOF table from the bc mesh subset for the given variable
+        // and component id.
+        _dof_table_boundary.reset(dof_table_bulk.deriveBoundaryConstrainedMap(
+            variable_id, {component_id}, std::move(boundary_mesh_subset)));
     }
 
     void getEssentialBCValues(
@@ -72,35 +87,42 @@ public:
 
         // Map boundary dof indices to bulk dof indices and the corresponding
         // values.
-        for (std::size_t i = 0; i < _boundary_mesh->getNumberOfNodes(); ++i)
+        for (auto const* const node : _boundary_mesh.getNodes())
         {
-            auto const bulk_node_id = _mapping_to_bulk_nodes.getComponent(i, 0);
-
-            MeshLib::Location const l{
-                _bulk_mesh_id, MeshLib::MeshItemType::Node, bulk_node_id};
-
-            auto const global_index = _dof_table_bulk.getGlobalIndex(
-                l, _variable_id_bulk, _component_id_bulk);
-            assert(global_index != NumLib::MeshComponentMap::nop);
-
-            bc_values.ids.push_back(global_index);
-            bc_values.values.push_back(_values.getComponent(i, 0));
+            auto const node_id = node->getID();
+            auto const global_index = _dof_table_boundary->getGlobalIndex(
+                {_boundary_mesh.getID(), MeshLib::MeshItemType::Node, node_id},
+                _variable_id, _component_id);
+            if (global_index == NumLib::MeshComponentMap::nop)
+                continue;
+            // For the DDC approach (e.g. with PETSc option), the negative index
+            // of global_index means that the entry by that index is a ghost
+            // one, which should be dropped. Especially for PETSc routines
+            // MatZeroRows and MatZeroRowsColumns, which are called to apply the
+            // Dirichlet BC, the negative index is not accepted like other
+            // matrix or vector PETSc routines. Therefore, the following
+            // if-condition is applied.
+            if (global_index >= 0)
+            {
+                bc_values.ids.emplace_back(global_index);
+                bc_values.values.push_back(_values[node_id]);
+                DBUG("global index %d, value %g", global_index,
+                     _values[node_id]);
+            }
         }
     }
 
 private:
-    std::size_t _bulk_mesh_id;
     MeshLib::PropertyVector<double> const& _values;
-    std::unique_ptr<MeshLib::Mesh> _boundary_mesh;
-    MeshLib::PropertyVector<std::size_t> const& _mapping_to_bulk_nodes;
-    NumLib::LocalToGlobalIndexMap const& _dof_table_bulk;
-    int const _variable_id_bulk;
-    int const _component_id_bulk;
+    MeshLib::Mesh const& _boundary_mesh;
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap const> _dof_table_boundary;
+    int const _variable_id;
+    int const _component_id;
 };
 
 std::unique_ptr<NonuniformDirichletBoundaryCondition>
 createNonuniformDirichletBoundaryCondition(
-    BaseLib::ConfigTree const& config,
+    BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh,
     NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
     int const component_id, const MeshLib::Mesh& bulk_mesh);
 

ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp

@@ -16,7 +16,7 @@ namespace ProcessLib
 {
 std::unique_ptr<NonuniformNeumannBoundaryCondition>
 createNonuniformNeumannBoundaryCondition(
-    BaseLib::ConfigTree const& config,
+    BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh,
     NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
     int const component_id, unsigned const integration_order,
     unsigned const shapefunction_order, MeshLib::Mesh const& bulk_mesh)
@@ -25,33 +25,28 @@ createNonuniformNeumannBoundaryCondition(
     //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
     config.checkConfigParameter("type", "NonuniformNeumann");
 
-    // TODO handle paths correctly
-    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__NonuniformNeumann__mesh}
-    auto const mesh_file = config.getConfigParameter<std::string>("mesh");
-
-    std::unique_ptr<MeshLib::Mesh> boundary_mesh(
-        MeshLib::IO::readMeshFromFile(mesh_file));
-
-    if (!boundary_mesh)
+    // Surface mesh and bulk mesh must have equal axial symmetry flags!
+    if (boundary_mesh.isAxiallySymmetric() != bulk_mesh.isAxiallySymmetric())
     {
-        OGS_FATAL("Error reading mesh `%s'", mesh_file.c_str());
+        OGS_FATAL(
+            "The boundary mesh %s axially symmetric but the bulk mesh %s. Both "
+            "must have an equal axial symmetry property.",
+            boundary_mesh.isAxiallySymmetric() ? "is" : "is not",
+            bulk_mesh.isAxiallySymmetric() ? "is" : "is not");
     }
 
-    // Surface mesh and bulk mesh must have equal axial symmetry flags!
-    boundary_mesh->setAxiallySymmetric(bulk_mesh.isAxiallySymmetric());
-
     // TODO finally use ProcessLib::Parameter here
     auto const field_name =
         //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__NonuniformNeumann__field_name}
         config.getConfigParameter<std::string>("field_name");
 
     auto const* const property =
-        boundary_mesh->getProperties().getPropertyVector<double>(field_name);
+        boundary_mesh.getProperties().getPropertyVector<double>(field_name);
 
     if (!property)
     {
         OGS_FATAL("A property with name `%s' does not exist in `%s'.",
-                  field_name.c_str(), mesh_file.c_str());
+                  field_name.c_str(), boundary_mesh.getName().c_str());
     }
 
     if (property->getMeshItemType() != MeshLib::MeshItemType::Node)
@@ -67,14 +62,13 @@ createNonuniformNeumannBoundaryCondition(
         OGS_FATAL("`%s' is not a one-component field.", field_name.c_str());
     }
 
-    std::string const mapping_to_bulk_nodes_property = "OriginalSubsurfaceNodeIDs";
+    std::string const mapping_to_bulk_nodes_property = "bulk_node_ids";
     auto const* const mapping_to_bulk_nodes =
-        boundary_mesh->getProperties().getPropertyVector<std::size_t>(
+        boundary_mesh.getProperties().getPropertyVector<std::size_t>(
             mapping_to_bulk_nodes_property);
 
-    if (!(mapping_to_bulk_nodes &&
-          mapping_to_bulk_nodes->getMeshItemType() ==
-              MeshLib::MeshItemType::Node) &&
+    if (!(mapping_to_bulk_nodes && mapping_to_bulk_nodes->getMeshItemType() ==
+                                       MeshLib::MeshItemType::Node) &&
         mapping_to_bulk_nodes->getNumberOfComponents() == 1)
     {
         OGS_FATAL("Field `%s' is not set up properly.",
@@ -82,8 +76,8 @@ createNonuniformNeumannBoundaryCondition(
     }
 
     return std::make_unique<NonuniformNeumannBoundaryCondition>(
-        integration_order, shapefunction_order, bulk_mesh.getDimension(),
-        std::move(boundary_mesh),
+        integration_order, shapefunction_order, dof_table, variable_id,
+        component_id, bulk_mesh.getDimension(), boundary_mesh,
         NonuniformNeumannBoundaryConditionData{
             *property, bulk_mesh.getID(), *mapping_to_bulk_nodes, dof_table,
             variable_id, component_id});

ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h

@@ -22,7 +22,7 @@ using NonuniformNeumannBoundaryCondition =
 
 std::unique_ptr<NonuniformNeumannBoundaryCondition>
 createNonuniformNeumannBoundaryCondition(
-    BaseLib::ConfigTree const& config,
+    BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh,
     NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
     int const component_id, unsigned const integration_order,
     unsigned const shapefunction_order, const MeshLib::Mesh& bulk_mesh);

ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h

@@ -12,6 +12,7 @@
 #include "MeshLib/PropertyVector.h"
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Function/Interpolation.h"
+#include "ProcessLib/Parameter/MeshNodeParameter.h"
 
 #include "GenericNonuniformNaturalBoundaryConditionLocalAssembler.h"
 
@@ -38,6 +39,7 @@ class NonuniformNeumannBoundaryConditionLocalAssembler final
 {
     using Base = GenericNonuniformNaturalBoundaryConditionLocalAssembler<
         ShapeFunction, IntegrationMethod, GlobalDim>;
+    using NodalVectorType = typename Base::NodalVectorType;
 
 public:
     /// The neumann_bc_value factor is directly integrated into the local
@@ -45,7 +47,7 @@ public:
     NonuniformNeumannBoundaryConditionLocalAssembler(
         MeshLib::Element const& e,
         std::size_t const local_matrix_size,
-        bool is_axially_symmetric,
+        bool const is_axially_symmetric,
         unsigned const integration_order,
         NonuniformNeumannBoundaryConditionData const& data)
         : Base(e, is_axially_symmetric, integration_order),
@@ -56,53 +58,35 @@ public:
 
     void assemble(std::size_t const id,
                   NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
-                  double const /*t*/, const GlobalVector& /*x*/,
+                  double const t, const GlobalVector& /*x*/,
                   GlobalMatrix& /*K*/, GlobalVector& b,
                   GlobalMatrix* /*Jac*/) override
     {
         _local_rhs.setZero();
 
-        auto indices = NumLib::getIndices(id, dof_table_boundary);
+        MeshNodeParameter<double> neumann_values{"NeumannValues", _data.values};
+        // Get element nodes for the interpolation from nodes to
+        // integration point.
+        NodalVectorType parameter_node_values =
+            neumann_values.getNodalValuesOnElement(Base::_element, t);
 
-        // TODO lots of heap allocations.
-        std::vector<double> neumann_param_nodal_values_local;
-        neumann_param_nodal_values_local.reserve(indices.size());
-        for (auto i : indices)
-        {
-            neumann_param_nodal_values_local.push_back(
-                _data.values.getComponent(i, 0));
-        }
-
-        auto const n_integration_points = Base::_ns_and_weights.size();
+        unsigned const n_integration_points =
+            Base::_integration_method.getNumberOfPoints();
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
             auto const& n_and_weight = Base::_ns_and_weights[ip];
-            double flux;
-            NumLib::shapeFunctionInterpolate(neumann_param_nodal_values_local,
-                                             n_and_weight.N, flux);
-            _local_rhs.noalias() +=
-                n_and_weight.N * (n_and_weight.weight * flux);
+            auto const& N = n_and_weight.N;
+            auto const& w = n_and_weight.weight;
+            _local_rhs.noalias() += N * parameter_node_values.dot(N) * w;
         }
 
-        // map boundary dof indices to bulk dof indices
-        for (auto& i : indices)
-        {
-            auto const bulk_node_id =
-                _data.mapping_to_bulk_nodes.getComponent(i, 0);
-
-            MeshLib::Location const l{
-                _data.bulk_mesh_id, MeshLib::MeshItemType::Node, bulk_node_id};
-
-            i = _data.dof_table_bulk.getGlobalIndex(l, _data.variable_id_bulk,
-                                                    _data.component_id_bulk);
-            assert(i != NumLib::MeshComponentMap::nop);
-        }
+        auto const indices = NumLib::getIndices(id, dof_table_boundary);
         b.add(indices, _local_rhs);
     }
 
 private:
     NonuniformNeumannBoundaryConditionData const& _data;
-    typename Base::NodalVectorType _local_rhs;
+    NodalVectorType _local_rhs;
 
 public:
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;

ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h

@@ -26,8 +26,7 @@ template <template <typename, typename, unsigned>
           class LocalAssemblerImplementation>
 NormalTractionBoundaryCondition<LocalAssemblerImplementation>::
     NormalTractionBoundaryCondition(
-        bool const is_axially_symmetric, unsigned const integration_order,
-        unsigned const shapefunction_order,
+        unsigned const integration_order, unsigned const shapefunction_order,
         NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
         int const variable_id, unsigned const global_dim,
         MeshLib::Mesh const& bc_mesh, Parameter<double> const& pressure)
@@ -55,7 +54,7 @@ NormalTractionBoundaryCondition<LocalAssemblerImplementation>::
 
     createLocalAssemblers<LocalAssemblerImplementation>(
         global_dim, _bc_mesh.getElements(), *_dof_table_boundary,
-        shapefunction_order, _local_assemblers, is_axially_symmetric,
+        shapefunction_order, _local_assemblers, _bc_mesh.isAxiallySymmetric(),
         _integration_order, _pressure);
 }
 
@@ -78,8 +77,8 @@ std::unique_ptr<NormalTractionBoundaryCondition<
 createNormalTractionBoundaryCondition(
     BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
     NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
-    bool is_axially_symmetric, unsigned const integration_order,
-    unsigned const shapefunction_order, unsigned const global_dim,
+    unsigned const integration_order, unsigned const shapefunction_order,
+    unsigned const global_dim,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters)
 {
     DBUG("Constructing NormalTractionBoundaryCondition from config.");
@@ -94,8 +93,8 @@ createNormalTractionBoundaryCondition(
     auto const& pressure = findParameter<double>(parameter_name, parameters, 1);
     return std::make_unique<NormalTractionBoundaryCondition<
         NormalTractionBoundaryConditionLocalAssembler>>(
-        is_axially_symmetric, integration_order, shapefunction_order, dof_table,
-        variable_id, global_dim, bc_mesh, pressure);
+        integration_order, shapefunction_order, dof_table, variable_id,
+        global_dim, bc_mesh, pressure);
 }
 
 }  // namespace NormalTractionBoundaryCondition

ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition.h

@@ -36,8 +36,7 @@ public:
     /// DOF-table, and a mesh subset for a given variable and its component.
     /// A local DOF-table, a subset of the given one, is constructed.
     NormalTractionBoundaryCondition(
-        bool const is_axially_symmetric, unsigned const integration_order,
-        unsigned const shapefunction_order,
+        unsigned const integration_order, unsigned const shapefunction_order,
         NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
         int const variable_id, unsigned const global_dim,
         MeshLib::Mesh const& bc_mesh, Parameter<double> const& pressure);
@@ -76,8 +75,8 @@ std::unique_ptr<NormalTractionBoundaryCondition<
 createNormalTractionBoundaryCondition(
     BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
     NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
-    bool is_axially_symmetric, unsigned const integration_order,
-    unsigned const shapefunction_order, unsigned const global_dim,
+    unsigned const integration_order, unsigned const shapefunction_order,
+    unsigned const global_dim,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters);
 
 }  // namespace NormalTractionBoundaryCondition

ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h

@@ -65,7 +65,7 @@ public:
     NormalTractionBoundaryConditionLocalAssembler(
         MeshLib::Element const& e,
         std::size_t const local_matrix_size,
-        bool is_axially_symmetric,
+        bool const is_axially_symmetric,
         unsigned const integration_order,
         Parameter<double> const& pressure)
         : _integration_method(integration_order), _pressure(pressure)

ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp

@@ -15,9 +15,8 @@ namespace ProcessLib
 std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
     BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
     NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
-    int const component_id, bool is_axially_symmetric,
-    unsigned const integration_order, unsigned const shapefunction_order,
-    unsigned const global_dim,
+    int const component_id, unsigned const integration_order,
+    unsigned const shapefunction_order, unsigned const global_dim,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters)
 {
     DBUG("Constructing RobinBcConfig from config.");
@@ -33,8 +32,8 @@ std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
     auto const& u_0 = findParameter<double>(u_0_name, parameters, 1);
 
     return std::make_unique<RobinBoundaryCondition>(
-        is_axially_symmetric, integration_order, shapefunction_order, dof_table,
-        variable_id, component_id, global_dim, bc_mesh,
+        integration_order, shapefunction_order, dof_table, variable_id,
+        component_id, global_dim, bc_mesh,
         RobinBoundaryConditionData{alpha, u_0});
 }
 

ProcessLib/BoundaryCondition/RobinBoundaryCondition.h

@@ -32,9 +32,8 @@ using RobinBoundaryCondition = GenericNaturalBoundaryCondition<
 std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
     BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
     NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
-    int const component_id, bool is_axially_symmetric,
-    unsigned const integration_order, unsigned const shapefunction_order,
-    unsigned const global_dim,
+    int const component_id, unsigned const integration_order,
+    unsigned const shapefunction_order, unsigned const global_dim,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters);
 
 }  // ProcessLib