diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
index a65f5a0de037a58a783cad943f9ecd9e86120f54..7f38b3d06593f1a098e98f714d0e2776b85c1831 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
@@ -55,8 +55,8 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
     auto const& mesh_subset =
         dof_table_bulk.getMeshSubset(variable_id, component_id);
 
-    MeshLib::MeshSubset bc_mesh_subset =
-        mesh_subset.getIntersectionByNodes(nodes);
+    _nodes_subset = nodesNodesIntersection(mesh_subset.getNodes(), nodes);
+    MeshLib::MeshSubset bc_mesh_subset(mesh_subset.getMesh(), &_nodes_subset);
 
     // Create local DOF table from intersected mesh subsets for the given
     // variable and component ids.
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
index 4cc07be9aac4c39bde0c69ed7f2dd8cfeaa0376f..62e452d51f47df750dd997d52cf3cf0f8bfd67db 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
@@ -55,6 +55,10 @@ private:
     /// defined.
     std::vector<MeshLib::Element*> _elements;
 
+    /// Intersection of boundary nodes and bulk mesh subset for the
+    /// variable_id/component_id pair.
+    std::vector<MeshLib::Node*> _nodes_subset;
+
     /// 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/NormalTractionBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h
index ce4475b2c8b4fd2b263d5fef5130fc923f926146..d4e85f390a0cac79e6c578e633243b8e1fa68892 100644
--- a/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h
@@ -43,7 +43,8 @@ NormalTractionBoundaryCondition<LocalAssemblerImplementation>::
     // variable.
     auto const& mesh_subset = dof_table_bulk.getMeshSubset(variable_id, 0);
 
-    auto bc_mesh_subset = mesh_subset.getIntersectionByNodes(nodes);
+    _nodes_subset = nodesNodesIntersection(mesh_subset.getNodes(), nodes);
+    MeshLib::MeshSubset bc_mesh_subset(mesh_subset.getMesh(), &_nodes_subset);
 
     // Create component ids vector for the current variable.
     auto const& number_of_components =
diff --git a/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition.h b/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition.h
index a9fed867a698fb0f55a5291cb0f5467aeb367417..1f59b291be495f09a6bb676c553b7c6c441e9fbb 100644
--- a/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition.h
@@ -57,6 +57,10 @@ private:
     /// defined.
     std::vector<MeshLib::Element*> _elements;
 
+    /// Intersection of boundary nodes and bulk mesh subset for the
+    /// variable_id/component_id pair.
+    std::vector<MeshLib::Node*> _nodes_subset;
+
     std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_all_nodes;
 
     /// Local dof table, a subset of the global one restricted to the
diff --git a/Tests/NumLib/LocalToGlobalIndexMap.cpp b/Tests/NumLib/LocalToGlobalIndexMap.cpp
index 7decd74106dbefc2c47cb2e29b0480aa5d73283d..c7b28fd212a8f33c52210e442401c82e7c5fcd9a 100644
--- a/Tests/NumLib/LocalToGlobalIndexMap.cpp
+++ b/Tests/NumLib/LocalToGlobalIndexMap.cpp
@@ -101,8 +101,10 @@ TEST_F(NumLibLocalToGlobalIndexMapTest, DISABLED_SubsetByComponent)
     // Find unique node ids of the selected elements for testing.
     std::vector<MeshLib::Node*> selected_nodes = MeshLib::getUniqueNodes(some_elements);
 
-    auto selected_component =
-        nodesSubset->getIntersectionByNodes(selected_nodes);
+    std::vector<MeshLib::Node*> nodes_intersection =
+        nodesNodesIntersection(nodesSubset->getNodes(), selected_nodes);
+    MeshLib::MeshSubset selected_component(nodesSubset->getMesh(),
+                                           &nodes_intersection);
 
     auto dof_map_subset = std::unique_ptr<NumLib::LocalToGlobalIndexMap>{
         dof_map->deriveBoundaryConstrainedMap(1,    // variable id
diff --git a/Tests/NumLib/LocalToGlobalIndexMapMultiComponent.cpp b/Tests/NumLib/LocalToGlobalIndexMapMultiComponent.cpp
index 329385f6b83cb7fefc8c18e678124d7faed98b53..4aa7ad8ef15e6631fcd1682088cf91a0775e71fa 100644
--- a/Tests/NumLib/LocalToGlobalIndexMapMultiComponent.cpp
+++ b/Tests/NumLib/LocalToGlobalIndexMapMultiComponent.cpp
@@ -69,8 +69,10 @@ public:
 
         std::vector<MeL::Node*> nodes = MeL::getUniqueNodes(boundary_elements);
 
+        nodes_subset =
+            nodesNodesIntersection(mesh_items_all_nodes->getNodes(), nodes);
         mesh_items_boundary = std::make_unique<MeshLib::MeshSubset>(
-            mesh_items_all_nodes->getIntersectionByNodes(nodes));
+            mesh_items_all_nodes->getMesh(), &nodes_subset);
     }
 
     ~NumLibLocalToGlobalIndexMapMultiDOFTest() override
@@ -150,6 +152,9 @@ public:
 
     std::unique_ptr<MeL::MeshSubset const> mesh_items_boundary;
     std::vector<MeL::Element*> boundary_elements;
+
+    /// Intersection of boundary nodes and bulk mesh subset.
+    std::vector<MeshLib::Node*> nodes_subset;
 };