diff --git a/NumLib/DOF/LocalToGlobalIndexMap.cpp b/NumLib/DOF/LocalToGlobalIndexMap.cpp
index 100cb8402f3430ee54db9b30a87933b203bd4a85..3be86d468591a29b2a892e0b171e0bea4f1f031c 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.cpp
+++ b/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,
diff --git a/NumLib/DOF/LocalToGlobalIndexMap.h b/NumLib/DOF/LocalToGlobalIndexMap.h
index e2c52d8400c653b5f864d7fb3403029e2b6c827c..7fe89ade313f426ab2a0cd502e4245b7acc196cc 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.h
+++ b/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.
diff --git a/NumLib/DOF/MeshComponentMap.cpp b/NumLib/DOF/MeshComponentMap.cpp
index c75cc86bf72b0f3b8e66d17c00daff3a40c3dbc6..d7793b299da79213bd2471d2bed8184953387aee 100644
--- a/NumLib/DOF/MeshComponentMap.cpp
+++ b/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);
}
diff --git a/NumLib/DOF/MeshComponentMap.h b/NumLib/DOF/MeshComponentMap.h
index 8a77ec4986160727f55172ea57bc95c410f06152..85558969d2f9d3911f1a2b43e5d6e3e471611a68 100644
--- a/NumLib/DOF/MeshComponentMap.h
+++ b/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