From d9a24a2f70a0558e382f3e48d96fcf8e2a84d33c Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Thu, 11 Jan 2018 15:09:44 +0100
Subject: [PATCH] [Coupling] Passed a vector of DOF tables to
 VectorMatrixAssembler

for flexibility in coupling.
---
 .../ComponentTransportProcess.cpp             |  10 +-
 .../CoupledSolutionsForStaggeredScheme.cpp    |  12 +-
 .../CoupledSolutionsForStaggeredScheme.h      |   9 +-
 .../GroundwaterFlowProcess.cpp                |  10 +-
 ProcessLib/HT/HTProcess.cpp                   |  19 ++-
 ProcessLib/HT/StaggeredHTFEM-impl.h           |  15 +-
 .../HeatConduction/HeatConductionProcess.cpp  |  10 +-
 .../HydroMechanics/HydroMechanicsFEM-impl.h   |   2 +-
 .../HydroMechanicsProcess-impl.h              |  24 +--
 .../HydroMechanics/HydroMechanicsProcess.cpp  |  10 +-
 .../SmallDeformationProcess.cpp               |  10 +-
 ProcessLib/LiquidFlow/LiquidFlowProcess.cpp   |  11 +-
 .../PhaseField/PhaseFieldProcess-impl.h       |  10 +-
 .../RichardsComponentTransportProcess.cpp     |  10 +-
 .../RichardsFlow/RichardsFlowProcess.cpp      |  10 +-
 .../SmallDeformationProcess-impl.h            |  10 +-
 ProcessLib/TES/TESProcess.cpp                 |  11 +-
 .../ThermalTwoPhaseFlowWithPPProcess.cpp      |  11 +-
 .../ThermoMechanicsProcess-impl.h             |  10 +-
 .../TwoPhaseFlowWithPPProcess.cpp             |  11 +-
 .../TwoPhaseFlowWithPrhoProcess.cpp           |  11 +-
 ProcessLib/VectorMatrixAssembler.cpp          | 141 ++++++------------
 ProcessLib/VectorMatrixAssembler.h            |  24 ++-
 23 files changed, 209 insertions(+), 192 deletions(-)

diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index a00794cb1bf..0d8a1a4eb2d 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -64,10 +64,12 @@ void ComponentTransportProcess::assembleConcreteProcess(
 {
     DBUG("Assemble ComponentTransportProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b, _coupled_solutions);
+        dof_table, t, x, M, K, b, _coupled_solutions);
 }
 
 void ComponentTransportProcess::assembleWithJacobianConcreteProcess(
@@ -77,11 +79,13 @@ void ComponentTransportProcess::assembleWithJacobianConcreteProcess(
 {
     DBUG("AssembleWithJacobian ComponentTransportProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac, _coupled_solutions, nullptr);
+        _local_assemblers, dof_table, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupled_solutions);
 }
 
 }  // namespace ComponentTransport
diff --git a/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp b/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
index a594ca6f390..9706bcfeaf5 100644
--- a/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
+++ b/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
@@ -30,9 +30,7 @@ CoupledSolutionsForStaggeredScheme::CoupledSolutionsForStaggeredScheme(
 
 std::vector<std::vector<double>> getPreviousLocalSolutions(
     const CoupledSolutionsForStaggeredScheme& cpl_xs,
-    const std::vector<
-        std::reference_wrapper<const std::vector<GlobalIndexType>>>&
-        indices)
+    const std::vector<std::vector<GlobalIndexType>>& indices)
 {
     if (cpl_xs.coupled_xs_t0.empty())
         return {};
@@ -44,7 +42,7 @@ std::vector<std::vector<double>> getPreviousLocalSolutions(
     int coupling_id = 0;
     for (auto const& x_t0 : cpl_xs.coupled_xs_t0)
     {
-        local_xs_t0.emplace_back(x_t0->get(indices[coupling_id].get()));
+        local_xs_t0.emplace_back(x_t0->get(indices[coupling_id]));
         coupling_id++;
     }
     return local_xs_t0;
@@ -52,9 +50,7 @@ std::vector<std::vector<double>> getPreviousLocalSolutions(
 
 std::vector<std::vector<double>> getCurrentLocalSolutions(
     const CoupledSolutionsForStaggeredScheme& cpl_xs,
-    const std::vector<
-        std::reference_wrapper<const std::vector<GlobalIndexType>>>&
-        indices)
+    const std::vector<std::vector<GlobalIndexType>>& indices)
 {
     const auto number_of_coupled_solutions = cpl_xs.coupled_xs.size();
     std::vector<std::vector<double>> local_xs_t1;
@@ -63,7 +59,7 @@ std::vector<std::vector<double>> getCurrentLocalSolutions(
     int coupling_id = 0;
     for (auto const& x_t1 : cpl_xs.coupled_xs)
     {
-        local_xs_t1.emplace_back(x_t1.get().get(indices[coupling_id].get()));
+        local_xs_t1.emplace_back(x_t1.get().get(indices[coupling_id]));
         coupling_id++;
     }
     return local_xs_t1;
diff --git a/ProcessLib/CoupledSolutionsForStaggeredScheme.h b/ProcessLib/CoupledSolutionsForStaggeredScheme.h
index f38f9509563..ee800490f71 100644
--- a/ProcessLib/CoupledSolutionsForStaggeredScheme.h
+++ b/ProcessLib/CoupledSolutionsForStaggeredScheme.h
@@ -12,7 +12,6 @@
 
 #pragma once
 
-#include <functional>
 #include <utility>
 #include <vector>
 
@@ -83,9 +82,7 @@ struct LocalCoupledSolutions
  */
 std::vector<std::vector<double>> getPreviousLocalSolutions(
     const CoupledSolutionsForStaggeredScheme& cpl_xs,
-    const std::vector<
-        std::reference_wrapper<const std::vector<GlobalIndexType>>>&
-        indices);
+    const std::vector<std::vector<GlobalIndexType>>& indices);
 
 /**
  * Fetch the nodal solutions of all coupled processes of the current time step
@@ -96,7 +93,5 @@ std::vector<std::vector<double>> getPreviousLocalSolutions(
  */
 std::vector<std::vector<double>> getCurrentLocalSolutions(
     const CoupledSolutionsForStaggeredScheme& cpl_xs,
-    const std::vector<
-        std::reference_wrapper<const std::vector<GlobalIndexType>>>&
-        indices);
+    const std::vector<std::vector<GlobalIndexType>>& indices);
 }  // end of ProcessLib
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
index 4cc79a22647..24c16c7fba1 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
@@ -67,10 +67,12 @@ void GroundwaterFlowProcess::assembleConcreteProcess(const double t,
 {
     DBUG("Assemble GroundwaterFlowProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b, _coupled_solutions);
+        dof_table, t, x, M, K, b, _coupled_solutions);
 }
 
 void GroundwaterFlowProcess::assembleWithJacobianConcreteProcess(
@@ -80,11 +82,13 @@ void GroundwaterFlowProcess::assembleWithJacobianConcreteProcess(
 {
     DBUG("AssembleWithJacobian GroundwaterFlowProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac, _coupled_solutions, nullptr);
+        _local_assemblers, dof_table, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupled_solutions);
 }
 
 }   // namespace GroundwaterFlow
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index ac379c210ab..96db3414705 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -85,9 +85,12 @@ void HTProcess::assembleConcreteProcess(const double t,
                                         GlobalMatrix& K,
                                         GlobalVector& b)
 {
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_tables;
     if (_use_monolithic_scheme)
     {
         DBUG("Assemble HTProcess.");
+        dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
     }
     else
     {
@@ -104,12 +107,14 @@ void HTProcess::assembleConcreteProcess(const double t,
                 "fluid flow process within HTProcess.");
         }
         setCoupledSolutionsOfPreviousTimeStep();
+        dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
+        dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
     }
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b, _coupled_solutions);
+        dof_tables, t, x, M, K, b, _coupled_solutions);
 }
 
 void HTProcess::assembleWithJacobianConcreteProcess(
@@ -119,16 +124,24 @@ void HTProcess::assembleWithJacobianConcreteProcess(
 {
     DBUG("AssembleWithJacobian HTProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_tables;
     if (!_use_monolithic_scheme)
     {
         setCoupledSolutionsOfPreviousTimeStep();
+        dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
+    }
+    else
+    {
+        dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
+        dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
     }
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac, _coupled_solutions, nullptr);
+        _local_assemblers, dof_tables, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupled_solutions);
 }
 
 void HTProcess::preTimestepConcreteProcess(GlobalVector const& x,
diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h
index bdd0c2074bb..7c12558ec2b 100644
--- a/ProcessLib/HT/StaggeredHTFEM-impl.h
+++ b/ProcessLib/HT/StaggeredHTFEM-impl.h
@@ -13,8 +13,6 @@
 
 #include "StaggeredHTFEM.h"
 
-#include <functional> // for std::reference_wrapper
-
 #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
 
 namespace ProcessLib
@@ -25,10 +23,10 @@ template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
     assembleForStaggeredScheme(double const t,
-                            std::vector<double>& local_M_data,
-                            std::vector<double>& local_K_data,
-                            std::vector<double>& local_b_data,
-                            LocalCoupledSolutions const& coupled_xs)
+                               std::vector<double>& local_M_data,
+                               std::vector<double>& local_K_data,
+                               std::vector<double>& local_b_data,
+                               LocalCoupledSolutions const& coupled_xs)
 {
     if (coupled_xs.process_id == 0)
     {
@@ -279,9 +277,8 @@ StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
 {
     auto const indices = NumLib::getIndices(this->_element.getID(), dof_table);
     assert(!indices.empty());
-    std::vector<std::reference_wrapper<const std::vector<GlobalIndexType>>>
-        indices_of_all_coupled_processes = {std::ref(indices),
-                                            std::ref(indices)};
+    std::vector<std::vector<GlobalIndexType>> indices_of_all_coupled_processes =
+        {indices, indices};
     auto const local_xs = getCurrentLocalSolutions(
         *(this->_coupled_solutions), indices_of_all_coupled_processes);
 
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.cpp b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
index 09db5281523..6213e31ffc2 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.cpp
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
@@ -78,10 +78,12 @@ void HeatConductionProcess::assembleConcreteProcess(const double t,
 {
     DBUG("Assemble HeatConductionProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b, _coupled_solutions);
+        dof_table, t, x, M, K, b, _coupled_solutions);
 }
 
 void HeatConductionProcess::assembleWithJacobianConcreteProcess(
@@ -91,11 +93,13 @@ void HeatConductionProcess::assembleWithJacobianConcreteProcess(
 {
     DBUG("AssembleWithJacobian HeatConductionProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac, _coupled_solutions, nullptr);
+        _local_assemblers, dof_table, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupled_solutions);
 }
 
 void HeatConductionProcess::computeSecondaryVariableConcrete(
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
index a48f5a440ea..0fabb9d555b 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
@@ -1,6 +1,6 @@
 /**
  * \copyright
- * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
  *            Distributed under a Modified BSD License.
  *              See accompanying file LICENSE.txt or
  *              http://www.opengeosys.org/project/license
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsProcess-impl.h
index 07c3a78e3b1..65f9fd4657a 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess-impl.h
@@ -235,7 +235,7 @@ void HydroMechanicsProcess<DisplacementDim>::initializeBoundaryConditions()
     {
         const int process_id_of_up = 0;
         initializeProcessBoundaryCondition(*_local_to_global_index_map,
-                                          process_id_of_up);
+                                           process_id_of_up);
         return;
     }
 
@@ -248,7 +248,7 @@ void HydroMechanicsProcess<DisplacementDim>::initializeBoundaryConditions()
     // for the equations of deformation.
     const int process_id_of_u = 1;
     initializeProcessBoundaryCondition(*_local_to_global_index_map,
-                                      process_id_of_u);
+                                       process_id_of_u);
 }
 
 template <int DisplacementDim>
@@ -262,10 +262,12 @@ void HydroMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
     // only the Newton-Raphson method is employed to simulate coupled HM
     // processes in this class, this function is actually not used so far.
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b, _coupled_solutions);
+        dof_table, t, x, M, K, b, _coupled_solutions);
 }
 
 template <int DisplacementDim>
@@ -283,11 +285,13 @@ void HydroMechanicsProcess<DisplacementDim>::
         DBUG(
             "Assemble the Jacobian of HydroMechanics for the monolithic"
             " scheme.");
+        std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+            dof_table = {std::ref(*_local_to_global_index_map)};
         // Call global assembler for each local assembly item.
         GlobalExecutor::executeMemberDereferenced(
             _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-            _local_assemblers, *_local_to_global_index_map, t, x, xdot,
-            dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions, nullptr);
+            _local_assemblers, dof_table, t, x, xdot, dxdot_dx, dx_dx, M, K, b,
+            Jac, _coupled_solutions);
         return;
     }
 
@@ -305,13 +309,13 @@ void HydroMechanicsProcess<DisplacementDim>::
             "HydroMechanics for the staggered scheme.");
     }
 
-    // Note: _local_to_global_index_map_with_base_nodes is asserted in
-    //       constructDofTable().
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_tables = {std::ref(*_local_to_global_index_map_with_base_nodes),
+                      std::ref(*_local_to_global_index_map)};
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac, _coupled_solutions,
-        _local_to_global_index_map_with_base_nodes.get());
+        _local_assemblers, dof_tables, t, x, xdot, dxdot_dx, dx_dx, M, K, b,
+        Jac, _coupled_solutions);
 }
 
 template <int DisplacementDim>
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
index 13cfc632538..549ff8016ad 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
@@ -508,10 +508,12 @@ void HydroMechanicsProcess<GlobalDim>::assembleConcreteProcess(
 {
     DBUG("Assemble HydroMechanicsProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b, _coupled_solutions);
+        dof_table, t, x, M, K, b, _coupled_solutions);
 }
 
 template <int GlobalDim>
@@ -523,10 +525,12 @@ void HydroMechanicsProcess<GlobalDim>::assembleWithJacobianConcreteProcess(
     DBUG("AssembleWithJacobian HydroMechanicsProcess.");
 
     // Call global assembler for each local assembly item.
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac, _coupled_solutions, nullptr);
+        _local_assemblers, dof_table, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupled_solutions);
 }
 
 template <int GlobalDim>
diff --git a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
index 209fdcc943a..0a570a64bab 100644
--- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
+++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
@@ -417,10 +417,12 @@ void SmallDeformationProcess<DisplacementDim>::assembleConcreteProcess(
 {
     DBUG("Assemble SmallDeformationProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b, _coupled_solutions);
+        dof_table, t, x, M, K, b, _coupled_solutions);
 }
 template <int DisplacementDim>
 void SmallDeformationProcess<DisplacementDim>::
@@ -434,10 +436,12 @@ void SmallDeformationProcess<DisplacementDim>::
     DBUG("AssembleWithJacobian SmallDeformationProcess.");
 
     // Call global assembler for each local assembly item.
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac, _coupled_solutions, nullptr);
+        _local_assemblers, dof_table, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupled_solutions);
 }
 template <int DisplacementDim>
 void SmallDeformationProcess<DisplacementDim>::preTimestepConcreteProcess(
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index 3815f637d3b..8c04ad6d3c6 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -85,10 +85,13 @@ void LiquidFlowProcess::assembleConcreteProcess(const double t,
                                                 GlobalVector& b)
 {
     DBUG("Assemble LiquidFlowProcess.");
+
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b, _coupled_solutions);
+        dof_table, t, x, M, K, b, _coupled_solutions);
 }
 
 void LiquidFlowProcess::assembleWithJacobianConcreteProcess(
@@ -98,11 +101,13 @@ void LiquidFlowProcess::assembleWithJacobianConcreteProcess(
 {
     DBUG("AssembleWithJacobian LiquidFlowProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac, _coupled_solutions, nullptr);
+        _local_assemblers, dof_table, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupled_solutions);
 }
 
 void LiquidFlowProcess::computeSecondaryVariableConcrete(const double t,
diff --git a/ProcessLib/PhaseField/PhaseFieldProcess-impl.h b/ProcessLib/PhaseField/PhaseFieldProcess-impl.h
index 6ea0104b361..a12755269c4 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcess-impl.h
+++ b/ProcessLib/PhaseField/PhaseFieldProcess-impl.h
@@ -142,10 +142,12 @@ void PhaseFieldProcess<DisplacementDim>::assembleConcreteProcess(
 {
     DBUG("Assemble PhaseFieldProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b, _coupled_solutions);
+        dof_table, t, x, M, K, b, _coupled_solutions);
 }
 
 template <int DisplacementDim>
@@ -156,11 +158,13 @@ void PhaseFieldProcess<DisplacementDim>::assembleWithJacobianConcreteProcess(
 {
     // DBUG("AssembleJacobian PhaseFieldProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac, _coupled_solutions, nullptr);
+        _local_assemblers, dof_table, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupled_solutions);
 }
 
 template <int DisplacementDim>
diff --git a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp
index b555c827ac8..e15b9e5e26c 100644
--- a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp
+++ b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp
@@ -71,10 +71,12 @@ void RichardsComponentTransportProcess::assembleConcreteProcess(
 {
     DBUG("Assemble RichardsComponentTransportProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b, _coupled_solutions);
+        dof_table, t, x, M, K, b, _coupled_solutions);
 }
 
 void RichardsComponentTransportProcess::assembleWithJacobianConcreteProcess(
@@ -84,11 +86,13 @@ void RichardsComponentTransportProcess::assembleWithJacobianConcreteProcess(
 {
     DBUG("AssembleWithJacobian RichardsComponentTransportProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac, _coupled_solutions, nullptr);
+        _local_assemblers, dof_table, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupled_solutions);
 }
 
 }  // namespace RichardsComponentTransport
diff --git a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
index 81c24f18de3..a3a49cf9df3 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
+++ b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
@@ -72,10 +72,12 @@ void RichardsFlowProcess::assembleConcreteProcess(
 {
     DBUG("Assemble RichardsFlowProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b, _coupled_solutions);
+        dof_table, t, x, M, K, b, _coupled_solutions);
 }
 
 void RichardsFlowProcess::assembleWithJacobianConcreteProcess(
@@ -85,11 +87,13 @@ void RichardsFlowProcess::assembleWithJacobianConcreteProcess(
 {
     DBUG("AssembleWithJacobian RichardsFlowProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac, _coupled_solutions, nullptr);
+        _local_assemblers, dof_table, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupled_solutions);
 }
 
 }  // namespace RichardsFlow
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h b/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h
index 368e052e70b..cb19a67400a 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h
@@ -148,10 +148,12 @@ void SmallDeformationProcess<DisplacementDim>::assembleConcreteProcess(
 {
     DBUG("Assemble SmallDeformationProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b, _coupled_solutions);
+        dof_table, t, x, M, K, b, _coupled_solutions);
 }
 
 template <int DisplacementDim>
@@ -163,11 +165,13 @@ void SmallDeformationProcess<DisplacementDim>::
 {
     DBUG("AssembleWithJacobian SmallDeformationProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac, _coupled_solutions, nullptr);
+        _local_assemblers, dof_table, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupled_solutions);
 
     b.copyValues(*_nodal_forces);
     std::transform(_nodal_forces->begin(), _nodal_forces->end(),
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
index 45db3b7bad8..217001080fd 100644
--- a/ProcessLib/TES/TESProcess.cpp
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -231,10 +231,12 @@ void TESProcess::assembleConcreteProcess(const double t,
 {
     DBUG("Assemble TESProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b, _coupled_solutions);
+        dof_table, t, x, M, K, b, _coupled_solutions);
 }
 
 void TESProcess::assembleWithJacobianConcreteProcess(
@@ -242,10 +244,13 @@ void TESProcess::assembleWithJacobianConcreteProcess(
     const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
     GlobalVector& b, GlobalMatrix& Jac)
 {
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
+    // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac, _coupled_solutions, nullptr);
+        _local_assemblers, dof_table, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupled_solutions);
 }
 
 void TESProcess::preTimestepConcreteProcess(GlobalVector const& x,
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
index 2a425087d52..fcb75282219 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
@@ -79,10 +79,13 @@ void ThermalTwoPhaseFlowWithPPProcess::assembleConcreteProcess(
     GlobalVector& b)
 {
     DBUG("Assemble ThermalTwoPhaseFlowWithPPProcess.");
+
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b, _coupled_solutions);
+        dof_table, t, x, M, K, b, _coupled_solutions);
 }
 
 void ThermalTwoPhaseFlowWithPPProcess::assembleWithJacobianConcreteProcess(
@@ -92,11 +95,13 @@ void ThermalTwoPhaseFlowWithPPProcess::assembleWithJacobianConcreteProcess(
 {
     DBUG("AssembleWithJacobian ThermalTwoPhaseFlowWithPPProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac, _coupled_solutions, nullptr);
+        _local_assemblers, dof_table, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupled_solutions);
 }
 void ThermalTwoPhaseFlowWithPPProcess::preTimestepConcreteProcess(
     GlobalVector const& x, double const t, double const delta_t,
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess-impl.h b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess-impl.h
index 4e056c74ca5..120310b627c 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess-impl.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess-impl.h
@@ -154,10 +154,12 @@ void ThermoMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
 {
     DBUG("Assemble ThermoMechanicsProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b, _coupled_solutions);
+        dof_table, t, x, M, K, b, _coupled_solutions);
 }
 
 template <int DisplacementDim>
@@ -171,11 +173,13 @@ void ThermoMechanicsProcess<DisplacementDim>::
 {
     DBUG("AssembleJacobian ThermoMechanicsProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac, _coupled_solutions, nullptr);
+        _local_assemblers, dof_table, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupled_solutions);
 }
 
 template <int DisplacementDim>
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
index 83d84df8102..4a6b65f77e0 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
@@ -77,10 +77,13 @@ void TwoPhaseFlowWithPPProcess::assembleConcreteProcess(const double t,
                                                         GlobalVector& b)
 {
     DBUG("Assemble TwoPhaseFlowWithPPProcess.");
+
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b, _coupled_solutions);
+        dof_table, t, x, M, K, b, _coupled_solutions);
 }
 
 void TwoPhaseFlowWithPPProcess::assembleWithJacobianConcreteProcess(
@@ -90,11 +93,13 @@ void TwoPhaseFlowWithPPProcess::assembleWithJacobianConcreteProcess(
 {
     DBUG("AssembleWithJacobian TwoPhaseFlowWithPPProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac, _coupled_solutions, nullptr);
+        _local_assemblers, dof_table, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupled_solutions);
 }
 
 }  // end of namespace
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
index 4ff49bdb543..1b878fe50e3 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
@@ -77,10 +77,13 @@ void TwoPhaseFlowWithPrhoProcess::assembleConcreteProcess(const double t,
                                                           GlobalVector& b)
 {
     DBUG("Assemble TwoPhaseFlowWithPrhoProcess.");
+
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b, _coupled_solutions);
+        dof_table, t, x, M, K, b, _coupled_solutions);
 }
 
 void TwoPhaseFlowWithPrhoProcess::assembleWithJacobianConcreteProcess(
@@ -90,11 +93,13 @@ void TwoPhaseFlowWithPrhoProcess::assembleWithJacobianConcreteProcess(
 {
     DBUG("AssembleWithJacobian TwoPhaseFlowWithPrhoProcess.");
 
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+       dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac, _coupled_solutions, nullptr);
+        _local_assemblers, dof_table, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupled_solutions);
 }
 void TwoPhaseFlowWithPrhoProcess::preTimestepConcreteProcess(
     GlobalVector const& x, double const t, double const dt,
diff --git a/ProcessLib/VectorMatrixAssembler.cpp b/ProcessLib/VectorMatrixAssembler.cpp
index f5b7aaae341..e81523fec81 100644
--- a/ProcessLib/VectorMatrixAssembler.cpp
+++ b/ProcessLib/VectorMatrixAssembler.cpp
@@ -40,12 +40,22 @@ void VectorMatrixAssembler::preAssemble(
 
 void VectorMatrixAssembler::assemble(
     const std::size_t mesh_item_id, LocalAssemblerInterface& local_assembler,
-    const NumLib::LocalToGlobalIndexMap& dof_table, const double t,
-    const GlobalVector& x, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
-    CoupledSolutionsForStaggeredScheme const* const cpl_xs)
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
+        dof_tables,
+    const double t, const GlobalVector& x, GlobalMatrix& M, GlobalMatrix& K,
+    GlobalVector& b, CoupledSolutionsForStaggeredScheme const* const cpl_xs)
 {
-    auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
+    std::vector<std::vector<GlobalIndexType>> indices_of_processes;
+    indices_of_processes.reserve(dof_tables.size());
+    for (std::size_t i = 0; i < dof_tables.size(); i++)
+    {
+        indices_of_processes.emplace_back(
+            NumLib::getIndices(mesh_item_id, dof_tables[i].get()));
+    }
 
+    auto const& indices = (dof_tables.size() == 1 && cpl_xs == nullptr)
+                              ? indices_of_processes[0]
+                              : indices_of_processes[cpl_xs->process_id];
     _local_M_data.clear();
     _local_K_data.clear();
     _local_b_data.clear();
@@ -58,24 +68,10 @@ void VectorMatrixAssembler::assemble(
     }
     else
     {
-        // Different processes in a staggered loop are allowed to use different
-        // orders of element. That means that the global indices can be
-        // different among different processes. The following vector stores the
-        // reference of the vectors of the global indices of all processes, and
-        // it is used to fetch the nodal solutions of all processes of the
-        // current element.
-        std::vector<std::reference_wrapper<const std::vector<GlobalIndexType>>>
-            indices_of_all_coupled_processes;
-        indices_of_all_coupled_processes.reserve(cpl_xs->coupled_xs.size());
-        for (std::size_t i = 0; i < cpl_xs->coupled_xs.size(); i++)
-        {
-            indices_of_all_coupled_processes.emplace_back(std::ref(indices));
-        }
-
-        auto local_coupled_xs0 = getPreviousLocalSolutions(
-            *cpl_xs, indices_of_all_coupled_processes);
+        auto local_coupled_xs0 =
+            getPreviousLocalSolutions(*cpl_xs, indices_of_processes);
         auto local_coupled_xs =
-            getCurrentLocalSolutions(*cpl_xs, indices_of_all_coupled_processes);
+            getCurrentLocalSolutions(*cpl_xs, indices_of_processes);
 
         ProcessLib::LocalCoupledSolutions local_coupled_solutions(
             cpl_xs->dt, cpl_xs->process_id, std::move(local_coupled_xs0),
@@ -109,22 +105,24 @@ void VectorMatrixAssembler::assemble(
 
 void VectorMatrixAssembler::assembleWithJacobian(
     std::size_t const mesh_item_id, LocalAssemblerInterface& local_assembler,
-    NumLib::LocalToGlobalIndexMap const& dof_table, const double t,
-    GlobalVector const& x, GlobalVector const& xdot, const double dxdot_dx,
-    const double dx_dx, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
-    GlobalMatrix& Jac, CoupledSolutionsForStaggeredScheme const* const cpl_xs,
-    NumLib::LocalToGlobalIndexMap const* const base_dof_table)
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
+        dof_tables,
+    const double t, GlobalVector const& x, GlobalVector const& xdot,
+    const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
+    GlobalVector& b, GlobalMatrix& Jac,
+    CoupledSolutionsForStaggeredScheme const* const cpl_xs)
 {
-    // If base_dof_table != nullptr, it means that the staggered scheme is
-    // applied for coupling, meanwhile DOF tables of different are different
-    // as well.
-    auto const indices =
-        ((base_dof_table == nullptr) ||
-         (cpl_xs->process_id ==
-              static_cast<int>(cpl_xs->coupled_xs.size()) - 1 &&
-          cpl_xs != nullptr))
-            ? NumLib::getIndices(mesh_item_id, dof_table)
-            : NumLib::getIndices(mesh_item_id, *base_dof_table);
+    std::vector<std::vector<GlobalIndexType>> indices_of_processes;
+    indices_of_processes.reserve(dof_tables.size());
+    for (std::size_t i = 0; i < dof_tables.size(); i++)
+    {
+        indices_of_processes.emplace_back(
+            NumLib::getIndices(mesh_item_id, dof_tables[i].get()));
+    }
+
+    auto const& indices = (dof_tables.size() == 1 && cpl_xs == nullptr)
+                              ? indices_of_processes[0]
+                              : indices_of_processes[cpl_xs->process_id];
     auto const local_xdot = xdot.get(indices);
 
     _local_M_data.clear();
@@ -141,32 +139,19 @@ void VectorMatrixAssembler::assembleWithJacobian(
     }
     else
     {
-        if (base_dof_table == nullptr)
-        {
-            local_assembleWithJacobianForStaggeredScheme(
-                t, indices, indices, local_xdot, local_assembler, dxdot_dx,
-                dx_dx, cpl_xs);
-        }
-        else
-        {
-            if (cpl_xs->process_id ==
-                static_cast<int>(cpl_xs->coupled_xs.size()) - 1)
-            {
-                const auto base_indices =
-                    NumLib::getIndices(mesh_item_id, *base_dof_table);
-                local_assembleWithJacobianForStaggeredScheme(
-                    t, base_indices, indices, local_xdot, local_assembler,
-                    dxdot_dx, dx_dx, cpl_xs);
-            }
-            else
-            {
-                const auto full_indices =
-                    NumLib::getIndices(mesh_item_id, dof_table);
-                local_assembleWithJacobianForStaggeredScheme(
-                    t, indices, full_indices, local_xdot, local_assembler,
-                    dxdot_dx, dx_dx, cpl_xs);
-            }
-        }
+        auto local_coupled_xs0 =
+            getPreviousLocalSolutions(*cpl_xs, indices_of_processes);
+        auto local_coupled_xs =
+            getCurrentLocalSolutions(*cpl_xs, indices_of_processes);
+
+        ProcessLib::LocalCoupledSolutions local_coupled_solutions(
+            cpl_xs->dt, cpl_xs->process_id, std::move(local_coupled_xs0),
+            std::move(local_coupled_xs));
+
+        _jacobian_assembler->assembleWithJacobianForStaggeredScheme(
+            local_assembler, t, local_xdot, dxdot_dx, dx_dx, _local_M_data,
+            _local_K_data, _local_b_data, _local_Jac_data,
+            local_coupled_solutions);
     }
 
     auto const num_r_c = indices.size();
@@ -202,36 +187,4 @@ void VectorMatrixAssembler::assembleWithJacobian(
     }
 }
 
-void VectorMatrixAssembler::local_assembleWithJacobianForStaggeredScheme(
-    const double t, std::vector<GlobalIndexType> const& base_indices,
-    std::vector<GlobalIndexType> const& full_indices,
-    std::vector<double> const& local_xdot,
-    LocalAssemblerInterface& local_assembler, const double dxdot_dx,
-    const double dx_dx, CoupledSolutionsForStaggeredScheme const* const cpl_xs)
-{
-    // The vector has the same purpose as that in assemble(..) in this file.
-    // For the detailed description, please see the comment inside assemble(..).
-    std::vector<std::reference_wrapper<const std::vector<GlobalIndexType>>>
-        indices_of_all_coupled_processes;
-    indices_of_all_coupled_processes.reserve(cpl_xs->coupled_xs.size());
-    for (std::size_t i = 0; i < cpl_xs->coupled_xs.size() - 1; i++)
-    {
-        indices_of_all_coupled_processes.emplace_back(std::ref(base_indices));
-    }
-    indices_of_all_coupled_processes.emplace_back(std::ref(full_indices));
-
-    auto local_coupled_xs0 =
-        getPreviousLocalSolutions(*cpl_xs, indices_of_all_coupled_processes);
-    auto local_coupled_xs =
-        getCurrentLocalSolutions(*cpl_xs, indices_of_all_coupled_processes);
-
-    ProcessLib::LocalCoupledSolutions local_coupled_solutions(
-        cpl_xs->dt, cpl_xs->process_id, std::move(local_coupled_xs0),
-        std::move(local_coupled_xs));
-
-    _jacobian_assembler->assembleWithJacobianForStaggeredScheme(
-        local_assembler, t, local_xdot, dxdot_dx, dx_dx, _local_M_data,
-        _local_K_data, _local_b_data, _local_Jac_data, local_coupled_solutions);
-}
-
 }  // namespace ProcessLib
diff --git a/ProcessLib/VectorMatrixAssembler.h b/ProcessLib/VectorMatrixAssembler.h
index 7961e7256c1..69ef7eec01a 100644
--- a/ProcessLib/VectorMatrixAssembler.h
+++ b/ProcessLib/VectorMatrixAssembler.h
@@ -44,7 +44,8 @@ public:
     //! \remark Jacobian is not assembled here, see assembleWithJacobian().
     void assemble(std::size_t const mesh_item_id,
                   LocalAssemblerInterface& local_assembler,
-                  NumLib::LocalToGlobalIndexMap const& dof_table,
+                  std::vector<std::reference_wrapper<
+                      NumLib::LocalToGlobalIndexMap>> const& dof_tables,
                   double const t, GlobalVector const& x, GlobalMatrix& M,
                   GlobalMatrix& K, GlobalVector& b,
                   CoupledSolutionsForStaggeredScheme const* const cpl_xs);
@@ -54,12 +55,13 @@ public:
     void assembleWithJacobian(
         std::size_t const mesh_item_id,
         LocalAssemblerInterface& local_assembler,
-        NumLib::LocalToGlobalIndexMap const& dof_table, const double t,
-        GlobalVector const& x, GlobalVector const& xdot, const double dxdot_dx,
-        const double dx_dx, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
-        GlobalMatrix& Jac,
-        CoupledSolutionsForStaggeredScheme const* const cpl_xs,
-        NumLib::LocalToGlobalIndexMap const* const base_dof_table);
+        std::vector<
+            std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
+            dof_tables,
+        const double t, GlobalVector const& x, GlobalVector const& xdot,
+        const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        CoupledSolutionsForStaggeredScheme const* const cpl_xs);
 
 private:
     // temporary data only stored here in order to avoid frequent memory
@@ -71,14 +73,6 @@ private:
 
     //! Used to assemble the Jacobian.
     std::unique_ptr<AbstractJacobianAssembler> _jacobian_assembler;
-
-    void local_assembleWithJacobianForStaggeredScheme(
-        const double t, std::vector<GlobalIndexType> const& base_indices,
-        std::vector<GlobalIndexType> const& full_indices,
-        std::vector<double> const& local_xdot,
-        LocalAssemblerInterface& local_assembler, const double dxdot_dx,
-        const double dx_dx,
-        CoupledSolutionsForStaggeredScheme const* const cpl_xs);
 };
 
 }  // namespace ProcessLib
-- 
GitLab