diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
index ead1f9bbc2bf9820a653d60844c20dd91bc45737..d9a9ad8c7612c35c4640458a2dffd382c4cca5cf 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
@@ -106,16 +106,13 @@ BoundaryConditionBuilder::createDirichletBoundaryCondition(
     // Sorted ids of all mesh_subsets.
     std::vector<std::size_t> sorted_nodes_ids;
 
-    auto const& mesh_subsets =
-        dof_table.getMeshSubsets(variable_id, *config.component_id);
-    for (auto const& mesh_subset : mesh_subsets)
-    {
-        auto const& nodes = mesh_subset->getNodes();
-        sorted_nodes_ids.reserve(sorted_nodes_ids.size() + nodes.size());
-        std::transform(std::begin(nodes), std::end(nodes),
-                       std::back_inserter(sorted_nodes_ids),
-                       [](MeshLib::Node* const n) { return n->getID(); });
-    }
+    auto const& mesh_subset =
+        dof_table.getMeshSubset(variable_id, *config.component_id);
+    auto const& nodes = mesh_subset.getNodes();
+    sorted_nodes_ids.reserve(sorted_nodes_ids.size() + nodes.size());
+    std::transform(std::begin(nodes), std::end(nodes),
+                   std::back_inserter(sorted_nodes_ids),
+                   [](MeshLib::Node* const n) { return n->getID(); });
     std::sort(std::begin(sorted_nodes_ids), std::end(sorted_nodes_ids));
 
     auto ids_new_end_iterator = std::end(ids);
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
index 44f6ad81fa8c7a44d78a9db5ba5c365b98dc26db..a65f5a0de037a58a783cad943f9ecd9e86120f54 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
@@ -52,19 +52,16 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
     DBUG("Found %d nodes for Natural BCs for the variable %d and component %d",
          nodes.size(), variable_id, component_id);
 
-    auto const& mesh_subsets =
-        dof_table_bulk.getMeshSubsets(variable_id, component_id);
+    auto const& mesh_subset =
+        dof_table_bulk.getMeshSubset(variable_id, component_id);
 
-    // TODO extend the node intersection to all parts of mesh_subsets, i.e.
-    // to each of the MeshSubset in the mesh_subsets.
-    _mesh_subset_all_nodes.reset(
-        mesh_subsets.getMeshSubset(0).getIntersectionByNodes(nodes));
-    MeshLib::MeshSubsets all_mesh_subsets{_mesh_subset_all_nodes.get()};
+    MeshLib::MeshSubset bc_mesh_subset =
+        mesh_subset.getIntersectionByNodes(nodes);
 
     // 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(bc_mesh_subset), _elements));
 
     createLocalAssemblers<LocalAssemblerImplementation>(
         global_dim, _elements, *_dof_table_boundary, shapefunction_order,
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
index 1955f8579fccae82c9d35681212f42a516d791cf..4cc07be9aac4c39bde0c69ed7f2dd8cfeaa0376f 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
@@ -55,8 +55,6 @@ private:
     /// defined.
     std::vector<MeshLib::Element*> _elements;
 
-    std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_all_nodes;
-
     /// Local dof table, a subset of the global one restricted to the
     /// participating #_elements of the boundary condition.
     std::unique_ptr<NumLib::LocalToGlobalIndexMap> _dof_table_boundary;
diff --git a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h
index 6d262a0e052cb96c89ea5538d5dd63da3d805857..28bbbb894153402f8df2fb1bb6018b70f362838d 100644
--- a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h
@@ -72,8 +72,7 @@ void GenericNonuniformNaturalBoundaryCondition<
     _mesh_subset_all_nodes.reset(
         new MeshLib::MeshSubset(*_boundary_mesh, &_boundary_mesh->getNodes()));
 
-    std::vector<MeshLib::MeshSubsets> all_mesh_subsets{
-        _mesh_subset_all_nodes.get()};
+    std::vector<MeshLib::MeshSubset> all_mesh_subsets{*_mesh_subset_all_nodes};
 
     std::vector<int> vec_var_n_components{1};
 
diff --git a/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h
index be32b566921550e2da715362ff6389acd7ec7444..ce4475b2c8b4fd2b263d5fef5130fc923f926146 100644
--- a/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h
@@ -41,13 +41,9 @@ NormalTractionBoundaryCondition<LocalAssemblerImplementation>::
 
     // Assume that the mesh subsets are equal for all components of the
     // variable.
-    auto const& mesh_subsets = dof_table_bulk.getMeshSubsets(variable_id, 0);
+    auto const& mesh_subset = dof_table_bulk.getMeshSubset(variable_id, 0);
 
-    // TODO extend the node intersection to all parts of mesh_subsets, i.e.
-    // to each of the MeshSubset in the mesh_subsets.
-    _mesh_subset_all_nodes.reset(
-        mesh_subsets.getMeshSubset(0).getIntersectionByNodes(nodes));
-    MeshLib::MeshSubsets all_mesh_subsets{_mesh_subset_all_nodes.get()};
+    auto bc_mesh_subset = mesh_subset.getIntersectionByNodes(nodes);
 
     // Create component ids vector for the current variable.
     auto const& number_of_components =
@@ -58,7 +54,7 @@ NormalTractionBoundaryCondition<LocalAssemblerImplementation>::
     // 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_ids, std::move(all_mesh_subsets), _elements));
+        variable_id, component_ids, std::move(bc_mesh_subset), _elements));
 
     createLocalAssemblers<LocalAssemblerImplementation>(
         global_dim, _elements, *_dof_table_boundary, shapefunction_order,
diff --git a/ProcessLib/SourceTerms/SourceTermBuilder.cpp b/ProcessLib/SourceTerms/SourceTermBuilder.cpp
index 156bb4c9ecad7c74771126ae23099856f2233047..9f2bfeb5c0db3ab7b423771b7723859b08abde12 100644
--- a/ProcessLib/SourceTerms/SourceTermBuilder.cpp
+++ b/ProcessLib/SourceTerms/SourceTermBuilder.cpp
@@ -62,16 +62,13 @@ std::unique_ptr<NodalSourceTerm> SourceTermBuilder::createNodalSourceTerm(
     // Sorted ids of all mesh_subsets.
     std::vector<std::size_t> sorted_nodes_ids;
 
-    auto const& mesh_subsets =
-        dof_table.getMeshSubsets(variable_id, *config.component_id);
-    for (auto const& mesh_subset : mesh_subsets)
-    {
-        auto const& nodes = mesh_subset->getNodes();
-        sorted_nodes_ids.reserve(sorted_nodes_ids.size() + nodes.size());
-        std::transform(std::begin(nodes), std::end(nodes),
-                       std::back_inserter(sorted_nodes_ids),
-                       [](MeshLib::Node* const n) { return n->getID(); });
-    }
+    auto const& mesh_subset =
+        dof_table.getMeshSubset(variable_id, *config.component_id);
+    auto const& nodes = mesh_subset.getNodes();
+    sorted_nodes_ids.reserve(sorted_nodes_ids.size() + nodes.size());
+    std::transform(std::begin(nodes), std::end(nodes),
+                   std::back_inserter(sorted_nodes_ids),
+                   [](MeshLib::Node* const n) { return n->getID(); });
     std::sort(std::begin(sorted_nodes_ids), std::end(sorted_nodes_ids));
 
     auto ids_new_end_iterator = std::end(ids);