[NL] Fixing the BCs DOF tables for BC-meshes.

NumLib/DOF/LocalToGlobalIndexMap.cpp

@@ -219,14 +219,17 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
 LocalToGlobalIndexMap* LocalToGlobalIndexMap::deriveBoundaryConstrainedMap(
     int const variable_id,
     std::vector<int> const& component_ids,
-    MeshLib::MeshSubset&& mesh_subset,
-    std::vector<MeshLib::Element*> const& elements) const
+    MeshLib::MeshSubset&& new_mesh_subset) const
 {
     DBUG("Construct reduced local to global index map.");
 
     if (component_ids.empty())
         OGS_FATAL("Expected non-empty vector of component ids.");
 
+    // Elements of the new_mesh_subset's mesh.
+    std::vector<MeshLib::Element*> const& elements =
+        new_mesh_subset.getMesh().getElements();
+
     // Create a subset of the current mesh component map.
     std::vector<int> global_component_ids;
 
@@ -234,15 +237,15 @@ LocalToGlobalIndexMap* LocalToGlobalIndexMap::deriveBoundaryConstrainedMap(
         global_component_ids.push_back(
             getGlobalComponent(variable_id, component_id));
 
-    auto mesh_component_map =
-        _mesh_component_map.getSubset(global_component_ids, mesh_subset);
+    auto mesh_component_map = _mesh_component_map.getSubset(
+        _mesh_subsets, new_mesh_subset, global_component_ids);
 
-    // Create copies of the mesh_subset for each of the global components.
+    // Create copies of the new_mesh_subset for each of the global components.
     // The last component is moved after the for-loop.
     std::vector<MeshLib::MeshSubset> all_mesh_subsets;
     for (int i = 0; i < static_cast<int>(global_component_ids.size()) - 1; ++i)
-        all_mesh_subsets.emplace_back(mesh_subset);
-    all_mesh_subsets.emplace_back(std::move(mesh_subset));
+        all_mesh_subsets.emplace_back(new_mesh_subset);
+    all_mesh_subsets.emplace_back(std::move(new_mesh_subset));
 
     return new LocalToGlobalIndexMap(std::move(all_mesh_subsets),
                                      global_component_ids, elements,

NumLib/DOF/LocalToGlobalIndexMap.h

@@ -82,16 +82,13 @@ public:
             vec_var_elements,
         NumLib::ComponentOrder const order);
 
-    /// Derive a LocalToGlobalIndexMap constrained to a set of mesh subsets and
-    /// elements. A new mesh component map will be constructed using the passed
-    /// mesh_subsets for the given variable and component ids.
-    ///
-    /// \note The elements are not necessarily those used in the mesh_subsets.
+    /// Derive a LocalToGlobalIndexMap constrained to the mesh subset and mesh
+    /// subset's elements. A new mesh component map will be constructed using
+    /// the passed mesh_subset for the given variable and component ids.
     LocalToGlobalIndexMap* deriveBoundaryConstrainedMap(
         int const variable_id,
         std::vector<int> const& component_ids,
-        MeshLib::MeshSubset&& mesh_subset,
-        std::vector<MeshLib::Element*> const& elements) const;
+        MeshLib::MeshSubset&& mesh_subset) const;
 
     /// Returns total number of degrees of freedom including those located in
     /// the ghost nodes.

NumLib/DOF/MeshComponentMap.cpp

@@ -116,21 +116,65 @@ MeshComponentMap::MeshComponentMap(
 #endif // end of USE_PETSC
 
 MeshComponentMap MeshComponentMap::getSubset(
-    std::vector<int> const& component_ids,
-    MeshLib::MeshSubset const& mesh_subset) const
+    std::vector<MeshLib::MeshSubset> const& bulk_mesh_subsets,
+    MeshLib::MeshSubset const& new_mesh_subset,
+    std::vector<int> const& new_global_component_ids) const
 {
+    {  // Testing first an assumption met later in the code that the meshes for
+       // the all bulk_mesh_subsets are equal.
+        auto const first_mismatch =
+            std::adjacent_find(begin(bulk_mesh_subsets), end(bulk_mesh_subsets),
+                               [](auto const& a, auto const& b) {
+                                   return a.getMeshID() != b.getMeshID();
+                               });
+        if (first_mismatch != end(bulk_mesh_subsets))
+        {
+            OGS_FATAL(
+                "Assumption in the MeshComponentMap violated. Expecting "
+                "all of mesh ids to be the same, but it is not true for "
+                "the mesh '%s' with id %d.",
+                first_mismatch->getMesh().getName().c_str(),
+                first_mismatch->getMeshID());
+        }
+    }
+
+    // Mapping of the nodes in the new_mesh_subset to the bulk mesh nodes
+    auto const& new_mesh_properties = new_mesh_subset.getMesh().getProperties();
+    if (!new_mesh_properties.template existsPropertyVector<std::size_t>(
+            "bulk_node_ids"))
+    {
+        OGS_FATAL(
+            "Bulk node ids map expected in the construction of the mesh "
+            "subset.");
+    }
+    auto const& bulk_node_ids_map =
+        *new_mesh_properties.template getPropertyVector<std::size_t>(
+            "bulk_node_ids");
+
     // New dictionary for the subset.
     ComponentGlobalIndexDict subset_dict;
 
-    std::size_t const mesh_id = mesh_subset.getMeshID();
+    std::size_t const new_mesh_id = new_mesh_subset.getMeshID();
     // Lookup the locations in the current mesh component map and
-    // insert the full lines into the subset dictionary.
-    for (std::size_t j = 0; j < mesh_subset.getNumberOfNodes(); j++)
-        for (auto component_id : component_ids)
-            subset_dict.insert(
-                getLine(Location(mesh_id, MeshLib::MeshItemType::Node,
-                                 mesh_subset.getNodeID(j)),
-                        component_id));
+    // insert the full lines into the new subset dictionary.
+    for (auto* const node : new_mesh_subset.getNodes())
+    {
+        auto const node_id = node->getID();
+
+        MeshLib::Location const new_location{
+            new_mesh_id, MeshLib::MeshItemType::Node, node_id};
+
+        // Assuming the meshes for the all bulk_mesh_subsets are equal.
+        MeshLib::Location const bulk_location{
+            bulk_mesh_subsets.front().getMeshID(), MeshLib::MeshItemType::Node,
+            bulk_node_ids_map[node_id]};
+
+        for (auto component_id : new_global_component_ids)
+        {
+            subset_dict.insert({new_location, component_id,
+                                getGlobalIndex(bulk_location, component_id)});
+        }
+    }
 
     return MeshComponentMap(subset_dict);
 }

NumLib/DOF/MeshComponentMap.h

@@ -44,12 +44,18 @@ public:
     /// The order (BY_LOCATION/BY_COMPONENT) of components is the same as of the
     /// current map.
     ///
-    /// \note For each component the same mesh_subset will be used.
+    /// \attention For each component the same new_mesh_subset will be used.
     ///
-    /// \param component_ids  The vector of global components id.
-    /// \param components components that should remain in the created subset
-    MeshComponentMap getSubset(std::vector<int> const& component_ids,
-                               MeshLib::MeshSubset const& mesh_subset) const;
+    /// \param bulk_mesh_subsets components that should remain in the created
+    /// subset, one for each global component.
+    /// \param new_mesh_subset The constraining mesh subset with a mapping of
+    /// node ids to the bulk mesh nodes.
+    /// \param new_global_component_ids The components for which the
+    /// bulk_mesh_subsets should be intersected with the new_mesh_subset.
+    MeshComponentMap getSubset(
+        std::vector<MeshLib::MeshSubset> const& bulk_mesh_subsets,
+        MeshLib::MeshSubset const& new_mesh_subset,
+        std::vector<int> const& new_global_component_ids) const;
 
     /// The number of dofs including the those located in the ghost nodes.
     std::size_t dofSizeWithGhosts() const