diff --git a/NumLib/DOF/LocalToGlobalIndexMap.cpp b/NumLib/DOF/LocalToGlobalIndexMap.cpp
index e79c3b0cbf261b8f64bad9cac68a6b5d59c5e503..c435af4b98211ca7dd2c884bc85c8c30efe0e58a 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.cpp
+++ b/NumLib/DOF/LocalToGlobalIndexMap.cpp
@@ -194,51 +194,71 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
     }
 }
 
-
 LocalToGlobalIndexMap::LocalToGlobalIndexMap(
     std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
-    int const component_id,
+    std::vector<int> 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.
+      _variable_component_offsets(
+          to_cumulative(std::vector<unsigned>(1, 1)))  // Single variable only.
 {
-    // There is only one mesh_subsets in the vector _mesh_subsets.
-    assert(_mesh_subsets.size() == 1);
-    auto const mss = *_mesh_subsets.front();
-
-    // 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)
+    // 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 (int i = 0; i < global_component_ids.size(); ++i)
     {
-        std::size_t const mesh_id = ms->getMeshID();
+        auto const& mss = *_mesh_subsets[i];
 
-        findGlobalIndices(elements.cbegin(), elements.cend(), ms->getNodes(), mesh_id,
-                          component_id, 0);  // There is only one component to
-                                             // write out, therefore the zero
-                                             // parameter.
+        // 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,
+    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.
-    auto const global_component_id =
-        getGlobalComponent(variable_id, component_id);
+    std::vector<int> 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_id, *mesh_subsets);
+        _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 (int i = 0; i < static_cast<int>(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_id, elements,
+                                     global_component_ids, elements,
                                      std::move(mesh_component_map));
 }
 
diff --git a/NumLib/DOF/LocalToGlobalIndexMap.h b/NumLib/DOF/LocalToGlobalIndexMap.h
index 0fe84a9084b213a975c54a1056eaee0449043d25..a5357bbb66836f5c5e1e78ba4e84a29680cfd144 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.h
+++ b/NumLib/DOF/LocalToGlobalIndexMap.h
@@ -85,7 +85,7 @@ public:
     /// \note The elements are not necessarily those used in the mesh_subsets.
     LocalToGlobalIndexMap* deriveBoundaryConstrainedMap(
         int const variable_id,
-        int const component_id,
+        std::vector<int> const& component_ids,
         std::unique_ptr<MeshLib::MeshSubsets>&& mesh_subsets,
         std::vector<MeshLib::Element*> const& elements) const;
 
@@ -170,7 +170,7 @@ private:
     /// this construtor.
     explicit LocalToGlobalIndexMap(
         std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
-        int const component_id,
+        std::vector<int> const& global_component_ids,
         std::vector<MeshLib::Element*> const& elements,
         NumLib::MeshComponentMap&& mesh_component_map);
 
@@ -190,8 +190,7 @@ private:
 
     /// The global component id for the specific variable (like velocity) and a
     /// component (like x, or y, or z).
-    std::size_t getGlobalComponent(int const variable_id,
-                                   int const component_id) const
+    int getGlobalComponent(int const variable_id, int const component_id) const
     {
         return _variable_component_offsets[variable_id] + component_id;
     }
diff --git a/NumLib/DOF/MeshComponentMap.cpp b/NumLib/DOF/MeshComponentMap.cpp
index 0228453d281f1aedb1525e3cf129d7e4085f2901..401a383447657c947db61b007cc15e044fb9f04a 100644
--- a/NumLib/DOF/MeshComponentMap.cpp
+++ b/NumLib/DOF/MeshComponentMap.cpp
@@ -30,8 +30,6 @@ GlobalIndexType const MeshComponentMap::nop =
 MeshComponentMap::MeshComponentMap(
     const std::vector<std::unique_ptr<MeshLib::MeshSubsets>>& components,
     ComponentOrder order)
-    : _num_components(components.size())
-
 {
     // get number of unknows
     GlobalIndexType num_unknowns = 0;
@@ -121,7 +119,6 @@ MeshComponentMap::MeshComponentMap(
 MeshComponentMap::MeshComponentMap(
     const std::vector<std::unique_ptr<MeshLib::MeshSubsets>>& components,
     ComponentOrder order)
-    : _num_components(components.size())
 {
     // construct dict (and here we number global_index by component type)
     GlobalIndexType global_index = 0;
@@ -149,10 +146,9 @@ MeshComponentMap::MeshComponentMap(
 #endif // end of USE_PETSC
 
 MeshComponentMap MeshComponentMap::getSubset(
-    std::size_t const component_id,
+    std::vector<int> const& component_ids,
     MeshLib::MeshSubsets const& mesh_subsets) const
 {
-    assert(component_id < _num_components);
     // New dictionary for the subset.
     ComponentGlobalIndexDict subset_dict;
 
@@ -162,18 +158,20 @@ MeshComponentMap MeshComponentMap::getSubset(
         // 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++)
-            subset_dict.insert(
-                getLine(Location(mesh_id, MeshLib::MeshItemType::Node,
-                                 mesh_subset->getNodeID(j)),
-                        component_id));
+            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++)
-            subset_dict.insert(
-                getLine(Location(mesh_id, MeshLib::MeshItemType::Cell,
-                                 mesh_subset->getElementID(j)),
-                        component_id));
+            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, 1);
+    return MeshComponentMap(subset_dict);
 }
 
 void MeshComponentMap::renumberByLocation(GlobalIndexType offset)
diff --git a/NumLib/DOF/MeshComponentMap.h b/NumLib/DOF/MeshComponentMap.h
index 465cd3127ff2cd02f6cc285f363a60c604cee0e0..9248a7461582397a83c271ef48d4b3dfe1484ac4 100644
--- a/NumLib/DOF/MeshComponentMap.h
+++ b/NumLib/DOF/MeshComponentMap.h
@@ -41,13 +41,15 @@ public:
         std::vector<std::unique_ptr<MeshLib::MeshSubsets>> const& components,
         ComponentOrder order);
 
-    /// Creates a single-component subset of the current mesh component map.
+    /// 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.
     ///
-    /// \param component_id  The global components id.
+    /// \note For each component the given MeshSubsets object 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::size_t const component_id,
+    MeshComponentMap getSubset(std::vector<int> const& component_ids,
                                MeshLib::MeshSubsets const& components) const;
 
     /// The number of dofs including the those located in the ghost nodes.
@@ -158,9 +160,8 @@ public:
 
 private:
     /// Private constructor used by internally created mesh component maps.
-    MeshComponentMap(detail::ComponentGlobalIndexDict& dict,
-                     unsigned const num_components)
-        : _dict(dict), _num_components(num_components)
+    explicit MeshComponentMap(detail::ComponentGlobalIndexDict& dict)
+        : _dict(dict)
     { }
 
     /// Looks up if a line is already stored in the dictionary.
@@ -182,10 +183,6 @@ private:
     std::size_t _num_global_dof = 0;
 #endif
 
-    /// Number of components
-    /// introduced mainly for error checking
-    unsigned const _num_components;
-
     /// Global ID for ghost entries
     std::vector<GlobalIndexType> _ghosts_indices;
 };
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
index 7613787b65e25914edf4b93512eb3f07ef1168f7..fe4949de469d88ac8051f343ac94978f5bfec567 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
@@ -55,7 +55,7 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
     // Create local DOF table from intersected mesh subsets for the given
     // variable and component ids.
     _dof_table_boundary.reset(dof_table_bulk.deriveBoundaryConstrainedMap(
-        variable_id, component_id, std::move(all_mesh_subsets), _elements));
+        variable_id, {component_id}, std::move(all_mesh_subsets), _elements));
 
     createLocalAssemblers<LocalAssemblerImplementation>(
         global_dim, _elements, *_dof_table_boundary, shapefunction_order, _local_assemblers,
diff --git a/ProcessLib/LIE/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h b/ProcessLib/LIE/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
index 7fe1d52b7c749971d9530ae945a05608adb3d84d..0509fe0d2427ebb6686d725af2c1da649f0bada5 100644
--- a/ProcessLib/LIE/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/LIE/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
@@ -62,7 +62,7 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
     // Create local DOF table from intersected mesh subsets for the given
     // variable and component ids.
     _dof_table_boundary.reset(dof_table_bulk.deriveBoundaryConstrainedMap(
-        variable_id, component_id, std::move(all_mesh_subsets), _elements));
+        variable_id, {component_id}, std::move(all_mesh_subsets), _elements));
 
     createLocalAssemblers<LocalAssemblerImplementation>(
         global_dim, _elements, *_dof_table_boundary, shapefunction_order, _local_assemblers,
diff --git a/Tests/NumLib/LocalToGlobalIndexMap.cpp b/Tests/NumLib/LocalToGlobalIndexMap.cpp
index f26ed7563f059509d9dd96b28ef95359fedc9b48..11ce39c3e04fc23bf16c840aef674218d2006324 100644
--- a/Tests/NumLib/LocalToGlobalIndexMap.cpp
+++ b/Tests/NumLib/LocalToGlobalIndexMap.cpp
@@ -107,8 +107,8 @@ TEST_F(NumLibLocalToGlobalIndexMapTest, DISABLED_SubsetByComponent)
         new MeshLib::MeshSubsets{selected_subset.get()}};
 
     auto dof_map_subset = std::unique_ptr<NumLib::LocalToGlobalIndexMap>{
-        dof_map->deriveBoundaryConstrainedMap(1,  // variable id
-                                              0,  // component id
+        dof_map->deriveBoundaryConstrainedMap(1,    // variable id
+                                              {0},  // component id
                                               std::move(selected_component),
                                               some_elements)};
 
diff --git a/Tests/NumLib/LocalToGlobalIndexMapMultiComponent.cpp b/Tests/NumLib/LocalToGlobalIndexMapMultiComponent.cpp
index 756fdf3c229771ee7b3c22a866a6d7b664f8b174..6777e610a70191f831bcafd9c205f8bd55fc88ee 100644
--- a/Tests/NumLib/LocalToGlobalIndexMapMultiComponent.cpp
+++ b/Tests/NumLib/LocalToGlobalIndexMapMultiComponent.cpp
@@ -77,7 +77,8 @@ public:
             delete e;
     }
 
-    void initComponents(const unsigned num_components, const unsigned selected_component,
+    void initComponents(const unsigned num_components,
+                        const int selected_component,
                         const NL::ComponentOrder order)
     {
         assert(selected_component < num_components);
@@ -97,7 +98,34 @@ public:
 
         dof_map_boundary.reset(dof_map->deriveBoundaryConstrainedMap(
             0,  // variable id
-            selected_component,
+            {selected_component},
+            std::move(components_boundary),
+            boundary_elements));
+    }
+
+    // Multi-component version.
+    void initComponents(unsigned const num_components,
+                        std::vector<int> const& selected_components,
+                        NL::ComponentOrder const order)
+    {
+        assert(selected_components.size() <= num_components);
+
+        std::vector<std::unique_ptr<MeshLib::MeshSubsets>> components;
+        for (unsigned i = 0; i < num_components; ++i)
+        {
+            components.emplace_back(
+                new MeL::MeshSubsets{mesh_items_all_nodes.get()});
+        }
+        std::vector<unsigned> vec_var_n_components(1, num_components);
+        dof_map.reset(new NL::LocalToGlobalIndexMap(
+            std::move(components), vec_var_n_components, order));
+
+        auto components_boundary = std::unique_ptr<MeshLib::MeshSubsets>{
+            new MeL::MeshSubsets{mesh_items_boundary.get()}};
+
+        dof_map_boundary.reset(dof_map->deriveBoundaryConstrainedMap(
+            0,  // variable id
+            selected_components,
             std::move(components_boundary),
             boundary_elements));
     }
@@ -107,6 +135,13 @@ public:
               std::function<std::size_t(std::size_t, std::size_t)> const&
                   compute_global_index);
 
+    // Multicomponent version
+    template <NL::ComponentOrder order>
+    void test(const unsigned num_components,
+              std::vector<int> const& selected_component,
+              std::function<std::size_t(std::size_t, std::size_t)> const&
+                  compute_global_index);
+
     std::unique_ptr<const MeshLib::Mesh> mesh;
     std::unique_ptr<const MeL::MeshSubset> mesh_items_all_nodes;
 
@@ -189,7 +224,10 @@ void NumLibLocalToGlobalIndexMapMultiDOFTest::test(
 
             ASSERT_EQ(2, global_idcs.size()); // boundary of quad is line with two nodes
 
-            for (unsigned n=0; n<2; ++n) // boundary of quad is line with two nodes
+            // e is the number of a boundary element (which is a line, so two
+            // nodes) from 0 to something. e+n must be the node ids along the
+            // x-axis of the mesh;
+            for (unsigned n = 0; n < 2; ++n)
             {
                 auto const node = e + n;
                 auto const glob_idx =
@@ -200,8 +238,70 @@ void NumLibLocalToGlobalIndexMapMultiDOFTest::test(
     }
 }
 
+// Multicomponent version
+template <NL::ComponentOrder ComponentOrder>
+void NumLibLocalToGlobalIndexMapMultiDOFTest::test(
+    unsigned const num_components,
+    std::vector<int> const& selected_components,
+    std::function<std::size_t(std::size_t, std::size_t)> const&
+        compute_global_index)
+{
+    initComponents(num_components, selected_components, ComponentOrder);
 
+    ASSERT_EQ(dof_map->getNumberOfComponents(), num_components);
+    ASSERT_EQ(dof_map->size(), mesh->getNumberOfElements());
 
+    ASSERT_EQ(dof_map_boundary->getNumberOfComponents(),
+              selected_components.size());
+    ASSERT_EQ(dof_map_boundary->size(), boundary_elements.size());
+
+    // check mesh elements
+    for (unsigned e = 0; e < dof_map->size(); ++e)
+    {
+        auto const element_nodes_size = mesh->getElement(e)->getNumberOfNodes();
+        auto const ptr_element_nodes = mesh->getElement(e)->getNodes();
+
+        for (unsigned c = 0; c < dof_map->getNumberOfComponents(); ++c)
+        {
+            auto const& global_idcs = (*dof_map)(e, c).rows;
+            ASSERT_EQ(element_nodes_size, global_idcs.size());
+
+            for (unsigned n = 0; n < element_nodes_size; ++n)
+            {
+                auto const node_id = ptr_element_nodes[n]->getID();
+                auto const glob_idx = compute_global_index(node_id, c);
+                EXPECT_EQ(glob_idx, global_idcs[n]);
+            }
+        }
+    }
+
+    // check boundary elements
+    for (unsigned e = 0; e < dof_map_boundary->size(); ++e)
+    {
+        ASSERT_EQ(selected_components.size(),
+                  dof_map_boundary->getNumberOfComponents());
+
+        for (unsigned c = 0; c < selected_components.size(); ++c)
+        {
+            auto const& global_idcs = (*dof_map_boundary)(e, c).rows;
+
+            ASSERT_EQ(
+                2,
+                global_idcs.size());  // boundary of quad is line with two nodes
+
+            // e is the number of a boundary element (which is a line, so two
+            // nodes) from 0 to something. e+n must be the node ids along the
+            // x-axis of the mesh;
+            for (unsigned n = 0; n < 2; ++n)
+            {
+                auto const node = e + n;
+                auto const glob_idx =
+                    compute_global_index(node, selected_components[c]);
+                EXPECT_EQ(glob_idx, global_idcs[n]);
+            }
+        }
+    }
+}
 
 void assert_equal(NL::LocalToGlobalIndexMap const& dof1, NL::LocalToGlobalIndexMap const& dof2)
 {
@@ -259,8 +359,42 @@ TEST_F(NumLibLocalToGlobalIndexMapMultiDOFTest, TestMultiCompByLocation)
 TEST_F(NumLibLocalToGlobalIndexMapMultiDOFTest, DISABLED_TestMultiCompByLocation)
 #endif
 {
-    unsigned const num_components = 5;
-    for (unsigned c = 0; c < num_components; ++c)
+    int const num_components = 5;
+    for (int c = 0; c < num_components; ++c)
         test<NL::ComponentOrder::BY_LOCATION>(
             num_components, c, ComputeGlobalIndexByLocation{num_components});
 }
+
+#ifndef USE_PETSC
+TEST_F(NumLibLocalToGlobalIndexMapMultiDOFTest, TestMultiComp2ByComponent)
+#else
+TEST_F(NumLibLocalToGlobalIndexMapMultiDOFTest,
+       DISABLED_TestMultiComp2ByComponent)
+#endif
+{
+    int const num_components = 5;
+    for (int c = 0; c < num_components; ++c)
+    {
+        std::vector<int> const components = {c, (c + 1) % num_components};
+        test<NL::ComponentOrder::BY_COMPONENT>(
+            num_components, components,
+            ComputeGlobalIndexByComponent{(mesh_subdivs + 1) *
+                                          (mesh_subdivs + 1)});
+    }
+}
+#ifndef USE_PETSC
+TEST_F(NumLibLocalToGlobalIndexMapMultiDOFTest, TestMultiComp2ByLocation)
+#else
+TEST_F(NumLibLocalToGlobalIndexMapMultiDOFTest,
+       DISABLED_TestMultiComp2ByLocation)
+#endif
+{
+    int const num_components = 5;
+    for (int c = 0; c < num_components; ++c)
+    {
+        std::vector<int> const components = {c, (c + 1) % num_components};
+        test<NL::ComponentOrder::BY_LOCATION>(
+            num_components, components,
+            ComputeGlobalIndexByLocation{num_components});
+    }
+}
diff --git a/Tests/NumLib/TestMeshComponentMap.cpp b/Tests/NumLib/TestMeshComponentMap.cpp
index 3cd9b5030149a25d1c42f5524268bd168bc0320c..213cc41a5103c2fc9ed0ad50ce2264a837fcf626 100644
--- a/Tests/NumLib/TestMeshComponentMap.cpp
+++ b/Tests/NumLib/TestMeshComponentMap.cpp
@@ -157,15 +157,14 @@ TEST_F(NumLibMeshComponentMapTest, DISABLED_SubsetOfNodesByComponent)
     for (std::size_t id : ids)
         some_nodes.push_back(const_cast<MeshLib::Node*>(mesh->getNode(id)));
 
-    MeshLib::MeshSubset some_nodes_mesh_subset(*mesh, &some_nodes);
+    MeshLib::MeshSubset const some_nodes_mesh_subset(*mesh, &some_nodes);
+    MeshLib::MeshSubsets const selected_component{&some_nodes_mesh_subset};
 
-    std::size_t const selected_component_id = 1;
-    auto selected_component = std::unique_ptr<MeshLib::MeshSubsets>{
-        new MeshLib::MeshSubsets{&some_nodes_mesh_subset}};
+    int const selected_component_id = 1;
 
     // Subset the original cmap.
     MeshComponentMap cmap_subset =
-        cmap->getSubset(selected_component_id, *selected_component);
+        cmap->getSubset({selected_component_id}, selected_component);
 
     // Check number of components as selected
     ASSERT_EQ(ids.size(), cmap_subset.dofSizeWithGhosts());
@@ -194,15 +193,14 @@ TEST_F(NumLibMeshComponentMapTest, DISABLED_SubsetOfNodesByLocation)
     for (std::size_t id : ids)
         some_nodes.push_back(const_cast<MeshLib::Node*>(mesh->getNode(id)));
 
-    MeshLib::MeshSubset some_nodes_mesh_subset(*mesh, &some_nodes);
+    MeshLib::MeshSubset const some_nodes_mesh_subset(*mesh, &some_nodes);
+    MeshLib::MeshSubsets const selected_component{&some_nodes_mesh_subset};
 
-    std::size_t const selected_component_id = 1;
-    auto selected_component = std::unique_ptr<MeshLib::MeshSubsets>{
-        new MeshLib::MeshSubsets{&some_nodes_mesh_subset}};
+    int const selected_component_id = 1;
 
     // Subset the original cmap.
     MeshComponentMap cmap_subset =
-        cmap->getSubset(selected_component_id, *selected_component);
+        cmap->getSubset({selected_component_id}, selected_component);
 
     // Check number of components as selected
     ASSERT_EQ(ids.size(), cmap_subset.dofSizeWithGhosts());
@@ -215,3 +213,78 @@ TEST_F(NumLibMeshComponentMapTest, DISABLED_SubsetOfNodesByLocation)
             cmap_subset.getGlobalIndex(l, comp1_id));
     }
 }
+
+#ifndef USE_PETSC
+TEST_F(NumLibMeshComponentMapTest, MulticomponentVariable)
+#else
+TEST_F(NumLibMeshComponentMapTest, DISABLED_MulticomponentVariable)
+#endif
+{
+    cmap =
+        new MeshComponentMap(components, NumLib::ComponentOrder::BY_LOCATION);
+
+    // Select some nodes from the full mesh.
+    std::array<std::size_t, 3> const ids = {{0, 5, 9}};
+    std::vector<MeshLib::Node*> some_nodes;
+    for (std::size_t id : ids)
+        some_nodes.push_back(const_cast<MeshLib::Node*>(mesh->getNode(id)));
+
+    MeshLib::MeshSubset const some_nodes_mesh_subset(*mesh, &some_nodes);
+    MeshLib::MeshSubsets const selected_component{&some_nodes_mesh_subset};
+
+    // Subset the original cmap.
+    std::vector<int> const selected_component_ids = {0, 1};
+    MeshComponentMap cmap_subset =
+        cmap->getSubset(selected_component_ids, selected_component);
+
+    // Check number of components as selected
+    ASSERT_EQ(ids.size() * selected_component_ids.size(),
+              cmap_subset.dofSizeWithGhosts());
+
+    // .. and the content of the subset.
+    for (std::size_t id : ids)
+    {
+        Location const l(mesh->getID(), MeshItemType::Node, id);
+        for (auto const& c : selected_component_ids)
+            EXPECT_EQ(cmap->getGlobalIndex(l, c),
+                      cmap_subset.getGlobalIndex(l, c));
+    }
+}
+
+#ifndef USE_PETSC
+TEST_F(NumLibMeshComponentMapTest, MulticomponentVariableSingleComponent)
+#else
+TEST_F(NumLibMeshComponentMapTest,
+       DISABLED_MulticomponentVariableSingleComponent)
+#endif
+{
+    cmap =
+        new MeshComponentMap(components, NumLib::ComponentOrder::BY_LOCATION);
+
+    // Select some nodes from the full mesh.
+    std::array<std::size_t, 3> const ids = {{0, 5, 9}};
+    std::vector<MeshLib::Node*> some_nodes;
+    for (std::size_t id : ids)
+        some_nodes.push_back(const_cast<MeshLib::Node*>(mesh->getNode(id)));
+
+    MeshLib::MeshSubset const some_nodes_mesh_subset(*mesh, &some_nodes);
+    MeshLib::MeshSubsets const selected_component{&some_nodes_mesh_subset};
+
+    // Subset the original cmap.
+    std::vector<int> const selected_component_ids = {1};
+    MeshComponentMap cmap_subset =
+        cmap->getSubset(selected_component_ids, selected_component);
+
+    // Check number of components as selected
+    ASSERT_EQ(ids.size() * selected_component_ids.size(),
+              cmap_subset.dofSizeWithGhosts());
+
+    // .. and the content of the subset.
+    for (std::size_t id : ids)
+    {
+        Location const l(mesh->getID(), MeshItemType::Node, id);
+        for (auto const& c : selected_component_ids)
+            EXPECT_EQ(cmap->getGlobalIndex(l, c),
+                      cmap_subset.getGlobalIndex(l, c));
+    }
+}