diff --git a/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp b/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
index 1da175880d096f59540cb93f4a6a2ffd560eae98..1c5c1fbc2f2239eaf01f4a3e5bd5aa31748f07dc 100644
--- a/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
+++ b/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
@@ -30,30 +30,38 @@ CoupledSolutionsForStaggeredScheme::CoupledSolutionsForStaggeredScheme(
 
 std::vector<std::vector<double>> getPreviousLocalSolutions(
     const CoupledSolutionsForStaggeredScheme& cpl_xs,
-    const std::vector<GlobalIndexType>& indices)
+    const std::vector<
+        std::reference_wrapper<const std::vector<GlobalIndexType>>>&
+        indices)
 {
     const auto number_of_coupled_solutions = cpl_xs.coupled_xs.size();
     std::vector<std::vector<double>> local_xs_t0;
     local_xs_t0.reserve(number_of_coupled_solutions);
 
+    int coupling_id = 0;
     for (auto const& x_t0 : cpl_xs.coupled_xs_t0)
     {
-        local_xs_t0.emplace_back(x_t0->get(indices));
+        local_xs_t0.emplace_back(x_t0->get(indices[coupling_id].get()));
+        coupling_id++;
     }
     return local_xs_t0;
 }
 
 std::vector<std::vector<double>> getCurrentLocalSolutions(
     const CoupledSolutionsForStaggeredScheme& cpl_xs,
-    const std::vector<GlobalIndexType>& indices)
+    const std::vector<
+        std::reference_wrapper<const std::vector<GlobalIndexType>>>&
+        indices)
 {
     const auto number_of_coupled_solutions = cpl_xs.coupled_xs.size();
     std::vector<std::vector<double>> local_xs_t1;
     local_xs_t1.reserve(number_of_coupled_solutions);
 
+    int coupling_id = 0;
     for (auto const& x_t1 : cpl_xs.coupled_xs)
     {
-        local_xs_t1.emplace_back(x_t1.get().get(indices));
+        local_xs_t1.emplace_back(x_t1.get().get(indices[coupling_id].get()));
+        coupling_id++;
     }
     return local_xs_t1;
 }
diff --git a/ProcessLib/CoupledSolutionsForStaggeredScheme.h b/ProcessLib/CoupledSolutionsForStaggeredScheme.h
index 272ceb44efe62afc1feb7d86cbf58caaa8d47ceb..91e1ab6755802da0554414b52b003decad3b7e9e 100644
--- a/ProcessLib/CoupledSolutionsForStaggeredScheme.h
+++ b/ProcessLib/CoupledSolutionsForStaggeredScheme.h
@@ -76,9 +76,13 @@ struct LocalCoupledSolutions
 
 std::vector<std::vector<double>> getPreviousLocalSolutions(
     const CoupledSolutionsForStaggeredScheme& cpl_xs,
-    const std::vector<GlobalIndexType>& indices);
+    const std::vector<
+        std::reference_wrapper<const std::vector<GlobalIndexType>>>&
+        indices);
 
 std::vector<std::vector<double>> getCurrentLocalSolutions(
     const CoupledSolutionsForStaggeredScheme& cpl_xs,
-    const std::vector<GlobalIndexType>& indices);
+    const std::vector<
+        std::reference_wrapper<const std::vector<GlobalIndexType>>>&
+        indices);
 }  // end of ProcessLib
diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h
index f13e73e70b834fb33683b49e335bbdbbf76d8701..9bd6806f4b475f06c401e25c0fb94af2e3a9adb9 100644
--- a/ProcessLib/HT/StaggeredHTFEM-impl.h
+++ b/ProcessLib/HT/StaggeredHTFEM-impl.h
@@ -13,6 +13,8 @@
 
 #include "StaggeredHTFEM.h"
 
+#include <functional> // for std::reference_wrapper
+
 #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
 
 namespace ProcessLib
@@ -277,8 +279,11 @@ StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
 {
     auto const indices = NumLib::getIndices(this->_element.getID(), dof_table);
     assert(!indices.empty());
-    auto const local_xs =
-        getCurrentLocalSolutions(*(this->_coupled_solutions), indices);
+    std::vector<std::reference_wrapper<const std::vector<GlobalIndexType>>>
+        indices_of_all_coupled_processes = {std::ref(indices),
+                                            std::ref(indices)};
+    auto const local_xs = getCurrentLocalSolutions(
+        *(this->_coupled_solutions), indices_of_all_coupled_processes);
 
     return this->getIntPtDarcyVelocityLocal(t, local_xs[0], local_xs[1], cache);
 }
diff --git a/ProcessLib/LocalAssemblerInterface.cpp b/ProcessLib/LocalAssemblerInterface.cpp
index b039bae02c7e53c588e1f875a0bc7a4e18a027a3..ab89bd751baec57fc6815f19ff1647951ee9ac4f 100644
--- a/ProcessLib/LocalAssemblerInterface.cpp
+++ b/ProcessLib/LocalAssemblerInterface.cpp
@@ -71,17 +71,11 @@ void LocalAssemblerInterface::computeSecondaryVariable(
 {
     auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
 
-    if (coupled_xs == nullptr)
-    {
-        auto const local_x = x.get(indices);
-        computeSecondaryVariableConcrete(t, local_x);
-    }
-    else
-    {
-        auto const local_coupled_xs =
-            getCurrentLocalSolutions(*coupled_xs, indices);
-        computeSecondaryVariableWithCoupledProcessConcrete(t, local_coupled_xs);
-    }
+    if (coupled_xs != nullptr)
+        return;
+
+    auto const local_x = x.get(indices);
+    computeSecondaryVariableConcrete(t, local_x);
 }
 
 void LocalAssemblerInterface::preTimestep(