From d9cfcb3dab38ada8c6155a040d5b55ba4b8814a5 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Mon, 19 Dec 2016 17:41:05 +0100
Subject: [PATCH] [NL] Add multicomponents to LocalToGlobalIndexMap.

---
 NumLib/DOF/LocalToGlobalIndexMap.cpp | 69 ++++++++++++++++++++++++++++
 NumLib/DOF/LocalToGlobalIndexMap.h   | 28 ++++++++++-
 2 files changed, 96 insertions(+), 1 deletion(-)

diff --git a/NumLib/DOF/LocalToGlobalIndexMap.cpp b/NumLib/DOF/LocalToGlobalIndexMap.cpp
index e79c3b0cbf2..8ad7646007f 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.cpp
+++ b/NumLib/DOF/LocalToGlobalIndexMap.cpp
@@ -221,6 +221,41 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
     }
 }
 
+LocalToGlobalIndexMap::LocalToGlobalIndexMap(
+    std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
+    std::vector<std::size_t> const& global_component_ids,
+    std::vector<MeshLib::Element*> const& elements,
+    NumLib::MeshComponentMap&& mesh_component_map)
+    : _mesh_subsets(std::move(mesh_subsets)),
+      _mesh_component_map(std::move(mesh_component_map)),
+      _variable_component_offsets(
+          to_cumulative(std::vector<unsigned>(1, 1)))  // Single variable only.
+{
+    // Each subset in the mesh_subsets represents a single component.
+    if (_mesh_subsets.size() != global_component_ids.size())
+        OGS_FATAL(
+            "Number of mesh subsets is not equal to number of components. "
+            "There are %d mesh subsets and %d components.",
+            mesh_subsets.size(), global_component_ids.size());
+
+    for (std::size_t i = 0; i < global_component_ids.size(); ++i)
+    {
+        auto const& mss = *_mesh_subsets[i];
+
+        // For all MeshSubset in mesh_subsets and each element of that
+        // MeshSubset
+        // save a line of global indices.
+        for (MeshLib::MeshSubset const* const ms : mss)
+        {
+            std::size_t const mesh_id = ms->getMeshID();
+
+            findGlobalIndices(elements.cbegin(), elements.cend(),
+                              ms->getNodes(), mesh_id, global_component_ids[i],
+                              i);
+        }
+    }
+}
+
 LocalToGlobalIndexMap* LocalToGlobalIndexMap::deriveBoundaryConstrainedMap(
     int const variable_id,
     int const component_id,
@@ -242,6 +277,40 @@ LocalToGlobalIndexMap* LocalToGlobalIndexMap::deriveBoundaryConstrainedMap(
                                      std::move(mesh_component_map));
 }
 
+LocalToGlobalIndexMap* LocalToGlobalIndexMap::deriveBoundaryConstrainedMap(
+    int const variable_id,
+    std::vector<int> const& component_ids,
+    std::unique_ptr<MeshLib::MeshSubsets>&& mesh_subsets,
+    std::vector<MeshLib::Element*> const& elements) const
+{
+    DBUG("Construct reduced local to global index map.");
+
+    if (component_ids.empty())
+        OGS_FATAL("Expected non-empty vector of component ids.");
+
+    // Create a subset of the current mesh component map.
+    std::vector<std::size_t> global_component_ids;
+
+    for (auto component_id : component_ids)
+        global_component_ids.push_back(
+            getGlobalComponent(variable_id, component_id));
+
+    auto mesh_component_map =
+        _mesh_component_map.getSubset(global_component_ids, *mesh_subsets);
+
+    // Create copies of the mesh_subsets for each of the global components.
+    // The last component is moved after the for-loop.
+    std::vector<std::unique_ptr<MeshLib::MeshSubsets>> all_mesh_subsets;
+    for (std::size_t i = 0; i < global_component_ids.size() - 1; ++i)
+        all_mesh_subsets.emplace_back(
+            new MeshLib::MeshSubsets{*mesh_subsets});
+    all_mesh_subsets.emplace_back(std::move(mesh_subsets));
+
+    return new LocalToGlobalIndexMap(std::move(all_mesh_subsets),
+                                     global_component_ids, elements,
+                                     std::move(mesh_component_map));
+}
+
 std::size_t
 LocalToGlobalIndexMap::dofSizeWithGhosts() const
 {
diff --git a/NumLib/DOF/LocalToGlobalIndexMap.h b/NumLib/DOF/LocalToGlobalIndexMap.h
index 0fe84a9084b..fbbba84a825 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.h
+++ b/NumLib/DOF/LocalToGlobalIndexMap.h
@@ -80,7 +80,9 @@ public:
 
     /// 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.
+    /// mesh_subsets for the given variable and component id.
+    ///
+    /// This is single-component version.
     ///
     /// \note The elements are not necessarily those used in the mesh_subsets.
     LocalToGlobalIndexMap* deriveBoundaryConstrainedMap(
@@ -89,6 +91,19 @@ public:
         std::unique_ptr<MeshLib::MeshSubsets>&& mesh_subsets,
         std::vector<MeshLib::Element*> const& elements) const;
 
+    /// 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.
+    ///
+    /// This is multi-component version.
+    ///
+    /// \note The elements are not necessarily those used in the mesh_subsets.
+    LocalToGlobalIndexMap* deriveBoundaryConstrainedMap(
+        int const variable_id,
+        std::vector<int> const& component_ids,
+        std::unique_ptr<MeshLib::MeshSubsets>&& mesh_subsets,
+        std::vector<MeshLib::Element*> const& elements) const;
+
     /// Returns total number of degrees of freedom including those located in
     /// the ghost nodes.
     std::size_t dofSizeWithGhosts() const;
@@ -174,6 +189,17 @@ private:
         std::vector<MeshLib::Element*> const& elements,
         NumLib::MeshComponentMap&& mesh_component_map);
 
+    /// Private constructor used by internally created local-to-global index
+    /// maps. The mesh_component_map is passed as argument instead of being
+    /// created by the constructor.
+    /// \attention The passed mesh_component_map is in undefined state after
+    /// this construtor.
+    explicit LocalToGlobalIndexMap(
+        std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
+        std::vector<std::size_t> const& global_component_ids,
+        std::vector<MeshLib::Element*> const& elements,
+        NumLib::MeshComponentMap&& mesh_component_map);
+
     template <typename ElementIterator>
     void
     findGlobalIndices(ElementIterator first, ElementIterator last,
-- 
GitLab