diff --git a/NumLib/DOF/DOFTableUtil.cpp b/NumLib/DOF/DOFTableUtil.cpp
index 67a961c0dbb243a5352b047b524a099db5d9fc0c..ceb39a6c51b7435354459a8bb66758e4eec919a8 100644
--- a/NumLib/DOF/DOFTableUtil.cpp
+++ b/NumLib/DOF/DOFTableUtil.cpp
@@ -24,20 +24,14 @@ double norm(GlobalVector const& x, unsigned const global_component,
 #endif
 
     double res = 0.0;
-    MeshLib::MeshSubsets const& mss =
-        dof_table.getMeshSubsets(global_component);
+    auto const& ms = dof_table.getMeshSubset(global_component);
 
-    for (unsigned i = 0; i < mss.size(); i++)
+    assert(ms.getMeshID() == mesh.getID());
+    for (MeshLib::Node const* node : ms.getNodes())
     {
-        MeshLib::MeshSubset const& ms = mss.getMeshSubset(i);
-        if (ms.getMeshID() != mesh.getID())
-            continue;
-        for (MeshLib::Node const* node : ms.getNodes())
-        {
-            auto const value = getNonGhostNodalValue(
-                x, mesh, dof_table, node->getID(), global_component);
-            res = calculate_norm(res, value);
-        }
+        auto const value = getNonGhostNodalValue(
+            x, mesh, dof_table, node->getID(), global_component);
+        res = calculate_norm(res, value);
     }
     return res;
 }
diff --git a/NumLib/DOF/DOFTableUtil.h b/NumLib/DOF/DOFTableUtil.h
index 633f227eff3bb7d8c84ef9dd8e251a2d917148a7..a880d1aa6ddeee0f4bfe5abb162afa08d0c50019 100644
--- a/NumLib/DOF/DOFTableUtil.h
+++ b/NumLib/DOF/DOFTableUtil.h
@@ -70,22 +70,17 @@ void transformVariableFromGlobalVector(
         local_to_global_index_map.getNumberOfVariableComponents(variable_id);
     for (int component = 0; component < n_components; ++component)
     {
-        auto const& mesh_subsets =
-            local_to_global_index_map.getMeshSubsets(variable_id, component);
-        assert(mesh_subsets.size() ==
-               1);  // Multiple meshes are not supported by the output_vector.
-        for (auto const& ms : mesh_subsets)
+        auto const& mesh_subset =
+            local_to_global_index_map.getMeshSubset(variable_id, component);
+        auto const mesh_id = mesh_subset.getMeshID();
+        for (auto const& node : mesh_subset.getNodes())
         {
-            auto const mesh_id = ms->getMeshID();
-            for (auto const& node : ms->getNodes())
-            {
-                auto const node_id = node->getID();
-                MeshLib::Location const l(mesh_id, MeshLib::MeshItemType::Node,
-                                          node_id);
-                output_vector.getComponent(node_id, component) = mapFunction(
-                    input_vector[local_to_global_index_map.getGlobalIndex(
-                        l, variable_id, component)]);
-            }
+            auto const node_id = node->getID();
+            MeshLib::Location const l(mesh_id, MeshLib::MeshItemType::Node,
+                                      node_id);
+            output_vector.getComponent(node_id, component) = mapFunction(
+                input_vector[local_to_global_index_map.getGlobalIndex(
+                    l, variable_id, component)]);
         }
     }
 }
diff --git a/NumLib/DOF/LocalToGlobalIndexMap.cpp b/NumLib/DOF/LocalToGlobalIndexMap.cpp
index a1375d44e58c5321a1fc8c51c3fbe307ce3801ab..c9fb761692c6b7988fbfc535a78bc5a3d3e54763 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.cpp
+++ b/NumLib/DOF/LocalToGlobalIndexMap.cpp
@@ -102,7 +102,7 @@ void LocalToGlobalIndexMap::findGlobalIndices(
 }
 
 LocalToGlobalIndexMap::LocalToGlobalIndexMap(
-    std::vector<MeshLib::MeshSubsets>&& mesh_subsets,
+    std::vector<MeshLib::MeshSubset>&& mesh_subsets,
     NumLib::ComponentOrder const order)
     : LocalToGlobalIndexMap(std::move(mesh_subsets),
                             std::vector<int>(mesh_subsets.size(), 1), order)
@@ -110,17 +110,14 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
 }
 
 LocalToGlobalIndexMap::LocalToGlobalIndexMap(
-    std::vector<MeshLib::MeshSubsets>&& mesh_subsets,
+    std::vector<MeshLib::MeshSubset>&& mesh_subsets,
     std::vector<int> const& vec_var_n_components,
     NumLib::ComponentOrder const order)
     : _mesh_subsets(std::move(mesh_subsets)),
       _mesh_component_map(_mesh_subsets, order),
       _variable_component_offsets(to_cumulative(vec_var_n_components))
 {
-    // For all MeshSubsets and each of their MeshSubset's and each element
-    // of that MeshSubset save a line of global indices.
-
-
+    // For each element of that MeshSubset save a line of global indices.
     std::size_t offset = 0;
     for (int variable_id = 0; variable_id < static_cast<int>(vec_var_n_components.size());
          ++variable_id)
@@ -131,24 +128,20 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
             auto const global_component_id =
                 getGlobalComponent(variable_id, component_id);
 
-            auto const& mss = _mesh_subsets[global_component_id];
-            for (int subset_id = 0; subset_id < static_cast<int>(mss.size());
-                 ++subset_id)
-            {
-                MeshLib::MeshSubset const& ms = mss.getMeshSubset(subset_id);
-                std::size_t const mesh_id = ms.getMeshID();
+            auto const& ms = _mesh_subsets[global_component_id];
+            std::size_t const mesh_id = ms.getMeshID();
 
-                findGlobalIndices(ms.elementsBegin(), ms.elementsEnd(), ms.getNodes(),
-                                  mesh_id, global_component_id, global_component_id);
-            }
+            findGlobalIndices(ms.elementsBegin(), ms.elementsEnd(),
+                              ms.getNodes(), mesh_id, global_component_id,
+                              global_component_id);
             // increase by number of components of that variable
-            offset += mss.size();
+            offset += _mesh_subsets.size();
         }
     }
 }
 
 LocalToGlobalIndexMap::LocalToGlobalIndexMap(
-    std::vector<MeshLib::MeshSubsets>&& mesh_subsets,
+    std::vector<MeshLib::MeshSubset>&& mesh_subsets,
     std::vector<int> const& vec_var_n_components,
     std::vector<std::vector<MeshLib::Element*> const*> const& vec_var_elements,
     NumLib::ComponentOrder const order)
@@ -158,8 +151,7 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
 {
     assert(vec_var_n_components.size() == vec_var_elements.size());
 
-    // For all MeshSubsets and each of their MeshSubset's and each element
-    // of that MeshSubset save a line of global indices.
+    // For each element of that MeshSubset save a line of global indices.
 
     // _rows should be resized based on an element ID
     std::size_t max_elem_id = 0;
@@ -181,24 +173,20 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
             auto const global_component_id =
                 getGlobalComponent(variable_id, component_id);
 
-            auto const& mss = _mesh_subsets[global_component_id];
-            for (int subset_id = 0; subset_id < static_cast<int>(mss.size());
-                 ++subset_id)
-            {
-                MeshLib::MeshSubset const& ms = mss.getMeshSubset(subset_id);
-                std::size_t const mesh_id = ms.getMeshID();
+            auto const& ms = _mesh_subsets[global_component_id];
+            std::size_t const mesh_id = ms.getMeshID();
 
-                findGlobalIndicesWithElementID(var_elements.cbegin(), var_elements.cend(), ms.getNodes(),
-                                               mesh_id, global_component_id, global_component_id);
-            }
+            findGlobalIndicesWithElementID(
+                var_elements.cbegin(), var_elements.cend(), ms.getNodes(),
+                mesh_id, global_component_id, global_component_id);
             // increase by number of components of that variable
-            offset += mss.size();
+            offset += _mesh_subsets.size();
         }
     }
 }
 
 LocalToGlobalIndexMap::LocalToGlobalIndexMap(
-    std::vector<MeshLib::MeshSubsets>&& mesh_subsets,
+    std::vector<MeshLib::MeshSubset>&& mesh_subsets,
     std::vector<int> const& global_component_ids,
     std::vector<MeshLib::Element*> const& elements,
     NumLib::MeshComponentMap&& mesh_component_map)
@@ -212,30 +200,25 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
         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());
+            _mesh_subsets.size(), global_component_ids.size());
 
     for (int i = 0; i < static_cast<int>(global_component_ids.size()); ++i)
     {
-        auto const& mss = _mesh_subsets[i];
+        auto const& ms = _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();
+        // MeshSubset save a line of global indices.
+        std::size_t const mesh_id = ms.getMeshID();
 
-            findGlobalIndices(elements.cbegin(), elements.cend(),
-                              ms->getNodes(), mesh_id, global_component_ids[i],
-                              i);
-        }
+        findGlobalIndices(elements.cbegin(), elements.cend(), ms.getNodes(),
+                          mesh_id, global_component_ids[i], i);
     }
 }
 
 LocalToGlobalIndexMap* LocalToGlobalIndexMap::deriveBoundaryConstrainedMap(
     int const variable_id,
     std::vector<int> const& component_ids,
-    MeshLib::MeshSubsets&& mesh_subsets,
+    MeshLib::MeshSubset&& mesh_subset,
     std::vector<MeshLib::Element*> const& elements) const
 {
     DBUG("Construct reduced local to global index map.");
@@ -251,14 +234,14 @@ LocalToGlobalIndexMap* LocalToGlobalIndexMap::deriveBoundaryConstrainedMap(
             getGlobalComponent(variable_id, component_id));
 
     auto mesh_component_map =
-        _mesh_component_map.getSubset(global_component_ids, mesh_subsets);
+        _mesh_component_map.getSubset(global_component_ids, mesh_subset);
 
-    // Create copies of the mesh_subsets for each of the global components.
+    // Create copies of the mesh_subset for each of the global components.
     // The last component is moved after the for-loop.
-    std::vector<MeshLib::MeshSubsets> all_mesh_subsets;
+    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_subsets);
-    all_mesh_subsets.emplace_back(std::move(mesh_subsets));
+        all_mesh_subsets.emplace_back(mesh_subset);
+    all_mesh_subsets.emplace_back(std::move(mesh_subset));
 
     return new LocalToGlobalIndexMap(std::move(all_mesh_subsets),
                                      global_component_ids, elements,
@@ -387,13 +370,13 @@ GlobalIndexType LocalToGlobalIndexMap::getLocalIndex(
                                              range_end);
 }
 
-MeshLib::MeshSubsets const& LocalToGlobalIndexMap::getMeshSubsets(
+MeshLib::MeshSubset const& LocalToGlobalIndexMap::getMeshSubset(
     int const variable_id, int const component_id) const
 {
-    return getMeshSubsets(getGlobalComponent(variable_id, component_id));
+    return getMeshSubset(getGlobalComponent(variable_id, component_id));
 }
 
-MeshLib::MeshSubsets const& LocalToGlobalIndexMap::getMeshSubsets(
+MeshLib::MeshSubset const& LocalToGlobalIndexMap::getMeshSubset(
     int const global_component_id) const
 {
     return _mesh_subsets[global_component_id];
diff --git a/NumLib/DOF/LocalToGlobalIndexMap.h b/NumLib/DOF/LocalToGlobalIndexMap.h
index 7677ccdbad7ef3f1f69a76cfa6912a3aca322693..e2c52d8400c653b5f864d7fb3403029e2b6c827c 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.h
+++ b/NumLib/DOF/LocalToGlobalIndexMap.h
@@ -18,10 +18,15 @@
 #include <Eigen/Dense>
 
 #include "MathLib/LinAlg/RowColumnIndices.h"
-#include "MeshLib/MeshSubsets.h"
 
 #include "MeshComponentMap.h"
 
+namespace MeshLib
+{
+class Node;
+class Element;
+}  // namespace MeshLib
+
 namespace NumLib
 {
 
@@ -46,9 +51,8 @@ public:
     ///
     /// \attention This constructor assumes the number of the given mesh subsets
     /// is equal to the number of variables, i.e. every variable has a single component.
-    LocalToGlobalIndexMap(
-        std::vector<MeshLib::MeshSubsets>&& mesh_subsets,
-        NumLib::ComponentOrder const order);
+    LocalToGlobalIndexMap(std::vector<MeshLib::MeshSubset>&& mesh_subsets,
+                          NumLib::ComponentOrder const order);
 
     /// Creates a MeshComponentMap internally and stores the global indices for
     /// each mesh element of the given mesh_subsets.
@@ -58,7 +62,7 @@ public:
     /// The size of the vector should be equal to the number of variables. Sum of the entries
     /// should be equal to the size of the mesh_subsets.
     /// \param order  type of ordering values in a vector
-    LocalToGlobalIndexMap(std::vector<MeshLib::MeshSubsets>&& mesh_subsets,
+    LocalToGlobalIndexMap(std::vector<MeshLib::MeshSubset>&& mesh_subsets,
                           std::vector<int> const& vec_var_n_components,
                           NumLib::ComponentOrder const order);
 
@@ -72,7 +76,7 @@ public:
     /// \param vec_var_elements  a vector of active mesh elements for each variable.
     /// \param order  type of ordering values in a vector
     LocalToGlobalIndexMap(
-        std::vector<MeshLib::MeshSubsets>&& mesh_subsets,
+        std::vector<MeshLib::MeshSubset>&& mesh_subsets,
         std::vector<int> const& vec_var_n_components,
         std::vector<std::vector<MeshLib::Element*> const*> const&
             vec_var_elements,
@@ -86,7 +90,7 @@ public:
     LocalToGlobalIndexMap* deriveBoundaryConstrainedMap(
         int const variable_id,
         std::vector<int> const& component_ids,
-        MeshLib::MeshSubsets&& mesh_subsets,
+        MeshLib::MeshSubset&& mesh_subset,
         std::vector<MeshLib::Element*> const& elements) const;
 
     /// Returns total number of degrees of freedom including those located in
@@ -137,10 +141,10 @@ public:
                                   std::size_t const range_begin,
                                   std::size_t const range_end) const;
 
-    MeshLib::MeshSubsets const& getMeshSubsets(int const variable_id,
-                                               int const component_id) const;
+    MeshLib::MeshSubset const& getMeshSubset(int const variable_id,
+                                             int const component_id) const;
 
-    MeshLib::MeshSubsets const& getMeshSubsets(
+    MeshLib::MeshSubset const& getMeshSubset(
         int const global_component_id) const;
 
 private:
@@ -150,7 +154,7 @@ private:
     /// \attention The passed mesh_component_map is in undefined state after
     /// this construtor.
     explicit LocalToGlobalIndexMap(
-        std::vector<MeshLib::MeshSubsets>&& mesh_subsets,
+        std::vector<MeshLib::MeshSubset>&& mesh_subsets,
         std::vector<int> const& global_component_ids,
         std::vector<MeshLib::Element*> const& elements,
         NumLib::MeshComponentMap&& mesh_component_map);
@@ -173,7 +177,7 @@ private:
 
 private:
     /// A vector of mesh subsets for each process variables' components.
-    std::vector<MeshLib::MeshSubsets> const _mesh_subsets;
+    std::vector<MeshLib::MeshSubset> _mesh_subsets;
     NumLib::MeshComponentMap _mesh_component_map;
 
     using Table = Eigen::Matrix<LineIndex, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
diff --git a/NumLib/DOF/MeshComponentMap.cpp b/NumLib/DOF/MeshComponentMap.cpp
index 066eef992bd4a8918b0f7ddcac0f712e35800f9e..815f8491de90bfe45dc90263337dbf361e1e7eec 100644
--- a/NumLib/DOF/MeshComponentMap.cpp
+++ b/NumLib/DOF/MeshComponentMap.cpp
@@ -13,7 +13,7 @@
 #include "MeshComponentMap.h"
 
 #include "BaseLib/Error.h"
-#include "MeshLib/MeshSubsets.h"
+#include "MeshLib/MeshSubset.h"
 
 #ifdef USE_PETSC
 #include "MeshLib/NodePartitionedMesh.h"
@@ -29,23 +29,16 @@ GlobalIndexType const MeshComponentMap::nop =
 
 #ifdef USE_PETSC
 MeshComponentMap::MeshComponentMap(
-    const std::vector<MeshLib::MeshSubsets>& components, ComponentOrder order)
+    std::vector<MeshLib::MeshSubset> const& components, ComponentOrder order)
 {
     // get number of unknows
     GlobalIndexType num_unknowns = 0;
     for (auto const& c : components)
     {
-        for (unsigned mesh_subset_index = 0; mesh_subset_index < c.size();
-             mesh_subset_index++)
-        {
-            MeshLib::MeshSubset const& mesh_subset =
-                c.getMeshSubset(mesh_subset_index);
             // PETSc always works with MeshLib::NodePartitionedMesh.
             const MeshLib::NodePartitionedMesh& mesh =
-                static_cast<const MeshLib::NodePartitionedMesh&>(
-                    mesh_subset.getMesh());
+                static_cast<const MeshLib::NodePartitionedMesh&>(c.getMesh());
             num_unknowns += mesh.getNumberOfGlobalNodes();
-        }
     }
 
     // construct dict (and here we number global_index by component type)
@@ -55,84 +48,75 @@ MeshComponentMap::MeshComponentMap(
     _num_local_dof = 0;
     for (auto const& c : components)
     {
-        for (unsigned mesh_subset_index = 0; mesh_subset_index < c.size();
-             mesh_subset_index++)
+        assert(dynamic_cast<MeshLib::NodePartitionedMesh const*>(
+                   &c.getMesh()) != nullptr);
+        std::size_t const mesh_id = c.getMeshID();
+        const MeshLib::NodePartitionedMesh& mesh =
+            static_cast<const MeshLib::NodePartitionedMesh&>(c.getMesh());
+
+        // mesh items are ordered first by node, cell, ....
+        for (std::size_t j = 0; j < c.getNumberOfNodes(); j++)
         {
-            MeshLib::MeshSubset const& mesh_subset =
-                c.getMeshSubset(mesh_subset_index);
-            assert(dynamic_cast<MeshLib::NodePartitionedMesh const*>(
-                       &mesh_subset.getMesh()) != nullptr);
-            std::size_t const mesh_id = mesh_subset.getMeshID();
-            const MeshLib::NodePartitionedMesh& mesh =
-                static_cast<const MeshLib::NodePartitionedMesh&>(
-                    mesh_subset.getMesh());
-
-            // mesh items are ordered first by node, cell, ....
-            for (std::size_t j = 0; j < mesh_subset.getNumberOfNodes(); j++)
+            GlobalIndexType global_id = 0;
+            if (order != ComponentOrder::BY_LOCATION)
             {
-                GlobalIndexType global_id = 0;
-                if (order != ComponentOrder::BY_LOCATION)
-                {
-                    // Deactivated since this case is not suitable to
-                    // arrange non ghost entries of a partition within
-                    // a rank in the parallel computing.
-                    OGS_FATAL("Global index in the system of equations"
-                              " can only be numbered by the order type"
-                              " of ComponentOrder::BY_LOCATION");
-                }
-                global_id = static_cast<GlobalIndexType>(
-                    components.size() * mesh.getGlobalNodeID(j)
-                        + comp_id);
-                const bool is_ghost =
-                    mesh.isGhostNode(mesh.getNode(j)->getID());
-                if (is_ghost)
-                {
-                    _ghosts_indices.push_back(global_id);
-                    global_id = -global_id;
-                    // If the ghost entry has an index of 0,
-                    // its index is set to the negative value of unknowns.
-                    if (global_id == 0) global_id = -num_unknowns;
-                }
-                else
-                    _num_local_dof++;
-
-                _dict.insert(
-                    Line(Location(mesh_id, MeshLib::MeshItemType::Node, j),
-                         comp_id, global_id));
+                // Deactivated since this case is not suitable to
+                // arrange non ghost entries of a partition within
+                // a rank in the parallel computing.
+                OGS_FATAL(
+                    "Global index in the system of equations"
+                    " can only be numbered by the order type"
+                    " of ComponentOrder::BY_LOCATION");
             }
+            global_id = static_cast<GlobalIndexType>(
+                components.size() * mesh.getGlobalNodeID(j) + comp_id);
+            const bool is_ghost = mesh.isGhostNode(mesh.getNode(j)->getID());
+            if (is_ghost)
+            {
+                _ghosts_indices.push_back(global_id);
+                global_id = -global_id;
+                // If the ghost entry has an index of 0,
+                // its index is set to the negative value of unknowns.
+                if (global_id == 0)
+                    global_id = -num_unknowns;
+            }
+            else
+                _num_local_dof++;
 
-            // Note: If the cells are really used (e.g. for the mixed FEM),
-            // the following global cell index must be reconsidered
-            // according to the employed cell indexing method.
-            for (std::size_t j = 0; j < mesh_subset.getNumberOfElements(); j++)
-                _dict.insert(
-                    Line(Location(mesh_id, MeshLib::MeshItemType::Cell, j),
-                         comp_id, cell_index++));
-
-            _num_global_dof += mesh.getNumberOfGlobalNodes();
+            _dict.insert(Line(Location(mesh_id, MeshLib::MeshItemType::Node, j),
+                              comp_id, global_id));
         }
+
+        // Note: If the cells are really used (e.g. for the mixed FEM),
+        // the following global cell index must be reconsidered
+        // according to the employed cell indexing method.
+        for (std::size_t j = 0; j < mesh_subset.getNumberOfElements(); j++)
+            _dict.insert(Line(Location(mesh_id, MeshLib::MeshItemType::Cell, j),
+                              comp_id, cell_index++));
+
+        _num_global_dof += mesh.getNumberOfGlobalNodes();
         comp_id++;
     }
 }
 #else
 MeshComponentMap::MeshComponentMap(
-    const std::vector<MeshLib::MeshSubsets>& components, ComponentOrder order)
+    std::vector<MeshLib::MeshSubset> const& components, ComponentOrder order)
 {
     // construct dict (and here we number global_index by component type)
     GlobalIndexType global_index = 0;
     std::size_t comp_id = 0;
     for (auto const& c : components)
     {
-        for (std::size_t mesh_subset_index = 0; mesh_subset_index < c.size(); mesh_subset_index++)
-        {
-            MeshLib::MeshSubset const& mesh_subset = c.getMeshSubset(mesh_subset_index);
-            std::size_t const mesh_id = mesh_subset.getMeshID();
-            // mesh items are ordered first by node, cell, ....
-            for (std::size_t j=0; j<mesh_subset.getNumberOfNodes(); j++)
-                _dict.insert(Line(Location(mesh_id, MeshLib::MeshItemType::Node, mesh_subset.getNodeID(j)), comp_id, global_index++));
-            for (std::size_t j=0; j<mesh_subset.getNumberOfElements(); j++)
-                _dict.insert(Line(Location(mesh_id, MeshLib::MeshItemType::Cell, mesh_subset.getElementID(j)), comp_id, global_index++));
-        }
+        std::size_t const mesh_id = c.getMeshID();
+        // mesh items are ordered first by node, cell, ....
+        for (std::size_t j = 0; j < c.getNumberOfNodes(); j++)
+            _dict.insert(Line(
+                Location(mesh_id, MeshLib::MeshItemType::Node, c.getNodeID(j)),
+                              comp_id, global_index++));
+        for (std::size_t j = 0; j < c.getNumberOfElements(); j++)
+            _dict.insert(Line(Location(mesh_id, MeshLib::MeshItemType::Cell,
+                                       c.getElementID(j)),
+                comp_id, global_index++));
         comp_id++;
     }
     _num_local_dof = _dict.size();
@@ -144,29 +128,26 @@ MeshComponentMap::MeshComponentMap(
 
 MeshComponentMap MeshComponentMap::getSubset(
     std::vector<int> const& component_ids,
-    MeshLib::MeshSubsets const& mesh_subsets) const
+    MeshLib::MeshSubset const& mesh_subset) const
 {
     // New dictionary for the subset.
     ComponentGlobalIndexDict subset_dict;
 
-    for (auto const& mesh_subset : mesh_subsets)
-    {
-        std::size_t const mesh_id = 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));
-        for (std::size_t j = 0; j < mesh_subset->getNumberOfElements(); j++)
-            for (auto component_id : component_ids)
-                subset_dict.insert(
-                    getLine(Location(mesh_id, MeshLib::MeshItemType::Cell,
-                                     mesh_subset->getElementID(j)),
-                            component_id));
-    }
+    std::size_t const mesh_id = 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));
+    for (std::size_t j = 0; j < mesh_subset.getNumberOfElements(); j++)
+        for (auto component_id : component_ids)
+            subset_dict.insert(
+                getLine(Location(mesh_id, MeshLib::MeshItemType::Cell,
+                                 mesh_subset.getElementID(j)),
+                        component_id));
 
     return MeshComponentMap(subset_dict);
 }
diff --git a/NumLib/DOF/MeshComponentMap.h b/NumLib/DOF/MeshComponentMap.h
index 78072e729c119d8de80a63785b38b99a1c4d2d52..926831e2e3ef45cc21516991a0d7358272252140 100644
--- a/NumLib/DOF/MeshComponentMap.h
+++ b/NumLib/DOF/MeshComponentMap.h
@@ -14,12 +14,9 @@
 
 #include "numlib_export.h"
 
-#include "ComponentGlobalIndexDict.h"
+#include "MeshLib/MeshSubset.h"
 
-namespace MeshLib
-{
-    class MeshSubsets;
-}
+#include "ComponentGlobalIndexDict.h"
 
 namespace NumLib
 {
@@ -40,19 +37,19 @@ public:
 public:
     /// \param components   a vector of components
     /// \param order        type of ordering values in a vector
-    MeshComponentMap(std::vector<MeshLib::MeshSubsets> const& components,
+    MeshComponentMap(std::vector<MeshLib::MeshSubset> const& components,
                      ComponentOrder order);
 
     /// Creates a multi-component subset of the current mesh component map.
     /// The order (BY_LOCATION/BY_COMPONENT) of components is the same as of the
     /// current map.
     ///
-    /// \note For each component the given MeshSubsets object will be used.
+    /// \note For each component the same 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::MeshSubsets const& components) const;
+                               MeshLib::MeshSubset const& mesh_subset) const;
 
     /// The number of dofs including the those located in the ghost nodes.
     std::size_t dofSizeWithGhosts() const