diff --git a/NumLib/DOF/MeshComponentMap.cpp b/NumLib/DOF/MeshComponentMap.cpp
index 0228453d281f1aedb1525e3cf129d7e4085f2901..341db59f6da77e9553524f3049b9219869b16f95 100644
--- a/NumLib/DOF/MeshComponentMap.cpp
+++ b/NumLib/DOF/MeshComponentMap.cpp
@@ -176,6 +176,36 @@ MeshComponentMap MeshComponentMap::getSubset(
     return MeshComponentMap(subset_dict, 1);
 }
 
+MeshComponentMap MeshComponentMap::getSubset(
+    std::vector<std::size_t> const& component_ids,
+    MeshLib::MeshSubsets const& mesh_subsets) const
+{
+    assert(component_ids.size() <= _num_components);
+    // 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 (std::size_t 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 (std::size_t component_id : component_ids)
+                subset_dict.insert(
+                    getLine(Location(mesh_id, MeshLib::MeshItemType::Cell,
+                                     mesh_subset->getElementID(j)),
+                            component_id));
+    }
+
+    return MeshComponentMap(subset_dict, 1);
+}
+
 void MeshComponentMap::renumberByLocation(GlobalIndexType offset)
 {
     GlobalIndexType global_index = offset;
diff --git a/NumLib/DOF/MeshComponentMap.h b/NumLib/DOF/MeshComponentMap.h
index 465cd3127ff2cd02f6cc285f363a60c604cee0e0..bdb9692056ead2a463852b780661c576e8acc885 100644
--- a/NumLib/DOF/MeshComponentMap.h
+++ b/NumLib/DOF/MeshComponentMap.h
@@ -45,9 +45,24 @@ public:
     /// The order (BY_LOCATION/BY_COMPONENT) of components is the same as of the
     /// current map.
     ///
+    /// This is single-component version.
+    ///
     /// \param component_id  The global components id.
     /// \param components components that should remain in the created subset
-    MeshComponentMap getSubset(std::size_t const component_id,
+    MeshComponentMap getSubset(std::size_t const component_ids,
+                               MeshLib::MeshSubsets const& components) const;
+
+    /// 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.
+    ///
+    /// This is multi-component version.
+    ///
+    /// \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<std::size_t> const& component_ids,
                                MeshLib::MeshSubsets const& components) const;
 
     /// The number of dofs including the those located in the ghost nodes.