diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 99b6af896dad227f9f089f18bfde597e854144d8..4e19d84b500630a61a7ac9cddf5cd046b1e07698 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -424,13 +424,12 @@ public:
                                    LocalCoupledSolutions const& coupled_xs)
     {
         auto local_p = Eigen::Map<const NodalVectorType>(
-            coupled_xs.local_coupled_xs[hydraulic_process_id].data(),
-            pressure_size);
+            &coupled_xs.local_coupled_xs[pressure_index], pressure_size);
         auto local_C = Eigen::Map<const NodalVectorType>(
-            coupled_xs.local_coupled_xs[first_transport_process_id].data(),
+            &coupled_xs.local_coupled_xs[first_concentration_index],
             concentration_size);
         auto local_C0 = Eigen::Map<const NodalVectorType>(
-            coupled_xs.local_coupled_xs0[first_transport_process_id].data(),
+            &coupled_xs.local_coupled_xs0[first_concentration_index],
             concentration_size);
 
         auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
@@ -535,14 +534,14 @@ public:
         LocalCoupledSolutions const& coupled_xs, int const transport_process_id)
     {
         auto local_C = Eigen::Map<const NodalVectorType>(
-            coupled_xs.local_coupled_xs[transport_process_id].data(),
+            &coupled_xs.local_coupled_xs[first_concentration_index +
+                                         (transport_process_id - 1) *
+                                             concentration_size],
             concentration_size);
         auto local_p = Eigen::Map<const NodalVectorType>(
-            coupled_xs.local_coupled_xs[hydraulic_process_id].data(),
-            pressure_size);
+            &coupled_xs.local_coupled_xs[pressure_index], pressure_size);
         auto local_p0 = Eigen::Map<const NodalVectorType>(
-            coupled_xs.local_coupled_xs0[hydraulic_process_id].data(),
-            pressure_size);
+            &coupled_xs.local_coupled_xs0[pressure_index], pressure_size);
 
         auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
             local_M_data, concentration_size, concentration_size);
@@ -731,12 +730,14 @@ public:
                 indices_of_all_coupled_processes(
                     _coupled_solutions->coupled_xs.size(), indices);
 
-            auto const local_xs = getCurrentLocalSolutions(
-                *(this->_coupled_solutions), indices_of_all_coupled_processes);
+            auto const local_xs =
+                getCoupledLocalSolutions(_coupled_solutions->coupled_xs,
+                                         indices_of_all_coupled_processes);
 
-            auto const local_p = MathLib::toVector(local_xs[hydraulic_process_id]);
-            auto const local_C =
-                MathLib::toVector(local_xs[first_transport_process_id]);
+            auto const local_p = Eigen::Map<const NodalVectorType>(
+                &local_xs[pressure_index], pressure_size);
+            auto const local_C = Eigen::Map<const NodalVectorType>(
+                &local_xs[first_concentration_index], concentration_size);
 
             return calculateIntPtDarcyVelocity(t, local_p, local_C, cache);
         }
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index b0776379699329ebb4bbe76fb1315f549af88d8f..6ed8579335b7c9b7960fb8eac0f8113897f28be5 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -167,8 +167,8 @@ Eigen::Vector3d ComponentTransportProcess::getFlux(std::size_t const element_id,
         std::vector<std::vector<GlobalIndexType>>
             indices_of_all_coupled_processes{
                 _coupled_solutions->coupled_xs.size(), r_c_indices.rows};
-        auto const local_xs = getCurrentLocalSolutions(
-            *(this->_coupled_solutions), indices_of_all_coupled_processes);
+        auto const local_xs = getCoupledLocalSolutions(
+            _coupled_solutions->coupled_xs, indices_of_all_coupled_processes);
 
         return _local_assemblers[element_id]->getFlux(p, t, local_xs);
 }
diff --git a/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp b/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
index 2e780eb6117a0c4388a764e643f8b76cff5866ba..de688cf32f66ecf8f3ed78dac4f34e44ba4a2e5c 100644
--- a/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
+++ b/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
@@ -12,6 +12,8 @@
 
 #include "CoupledSolutionsForStaggeredScheme.h"
 
+#include <numeric>
+
 #include "MathLib/LinAlg/LinAlg.h"
 #include "Process.h"
 
@@ -27,48 +29,34 @@ CoupledSolutionsForStaggeredScheme::CoupledSolutionsForStaggeredScheme(
     }
 }
 
-std::vector<std::vector<double>> getPreviousLocalSolutions(
-    const CoupledSolutionsForStaggeredScheme& cpl_xs,
-    const std::vector<std::vector<GlobalIndexType>>& indices)
+std::vector<double> getCoupledLocalSolutions(
+    std::vector<GlobalVector*> const& global_solutions,
+    std::vector<std::vector<GlobalIndexType>> const& indices)
 {
-    if (cpl_xs.coupled_xs_t0.empty())
+    if (global_solutions.empty())
     {
         return {};
     }
 
-    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)
+    std::size_t const local_solutions_size = std::accumulate(
+        cbegin(indices),
+        cend(indices),
+        std::size_t(0),
+        [](GlobalIndexType const size,
+           std::vector<GlobalIndexType> const& process_indices) {
+            return size + process_indices.size();
+        });
+    std::vector<double> local_solutions;
+    local_solutions.reserve(local_solutions_size);
+
+    int number_of_processes = static_cast<int>(global_solutions.size());
+    for (int process_id = 0; process_id < number_of_processes; ++process_id)
     {
-        local_xs_t0.emplace_back(x_t0->get(indices[coupling_id]));
-        coupling_id++;
+        auto values = global_solutions[process_id]->get(indices[process_id]);
+        local_solutions.insert(cend(local_solutions),
+                               std::make_move_iterator(begin(values)),
+                               std::make_move_iterator(end(values)));
     }
-    return local_xs_t0;
+    return local_solutions;
 }
-
-std::vector<std::vector<double>> getCurrentLocalSolutions(
-    const CoupledSolutionsForStaggeredScheme& cpl_xs,
-    const std::vector<std::vector<GlobalIndexType>>& indices)
-{
-    if (cpl_xs.coupled_xs.empty())
-    {
-        return {};
-    }
-
-    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(indices[coupling_id]));
-        coupling_id++;
-    }
-    return local_xs_t1;
-}
-
 }  // namespace ProcessLib
diff --git a/ProcessLib/CoupledSolutionsForStaggeredScheme.h b/ProcessLib/CoupledSolutionsForStaggeredScheme.h
index 3fe07a282e78be40b7ba2f2aaa84c054752703b8..cfce3363a5ecb74622f48c0d9acc458e5841126d 100644
--- a/ProcessLib/CoupledSolutionsForStaggeredScheme.h
+++ b/ProcessLib/CoupledSolutionsForStaggeredScheme.h
@@ -49,38 +49,25 @@ struct CoupledSolutionsForStaggeredScheme
  */
 struct LocalCoupledSolutions
 {
-    LocalCoupledSolutions(std::vector<std::vector<double>>&& local_coupled_xs0_,
-                          std::vector<std::vector<double>>&& local_coupled_xs_)
+    LocalCoupledSolutions(std::vector<double>&& local_coupled_xs0_,
+                          std::vector<double>&& local_coupled_xs_)
         : local_coupled_xs0(std::move(local_coupled_xs0_)),
           local_coupled_xs(std::move(local_coupled_xs_))
     {
     }
 
     /// Local solutions of the previous time step.
-    std::vector<std::vector<double>> const local_coupled_xs0;
+    std::vector<double> const local_coupled_xs0;
     /// Local solutions of the current time step.
-    std::vector<std::vector<double>> const local_coupled_xs;
+    std::vector<double> const local_coupled_xs;
 };
 
 /**
- * Fetch the nodal solutions of all coupled processes of the previous time step
- * of an element.
- * @param cpl_xs  Solutions of all coupled equations.
- * @param indices Nodal indices of an element.
- * @return Nodal solutions of the previous time step of an element
+ * Fetch the nodal solutions of all coupled processes from the given vector of
+ * global solutions for each process into a flat vector.
  */
-std::vector<std::vector<double>> getPreviousLocalSolutions(
-    const CoupledSolutionsForStaggeredScheme& cpl_xs,
-    const std::vector<std::vector<GlobalIndexType>>& indices);
+std::vector<double> getCoupledLocalSolutions(
+    std::vector<GlobalVector*> const& global_solutions,
+    std::vector<std::vector<GlobalIndexType>> const& indices);
 
-/**
- * Fetch the nodal solutions of all coupled processes of the current time step
- * of an element.
- * @param cpl_xs  Solutions of all coupled equations.
- * @param indices Nodal indices of an element.
- * @return Nodal solutions of the current time step of an element
- */
-std::vector<std::vector<double>> getCurrentLocalSolutions(
-    const CoupledSolutionsForStaggeredScheme& cpl_xs,
-    const std::vector<std::vector<GlobalIndexType>>& indices);
 }  // namespace ProcessLib
diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index 52d85f23ffb4c8b25a1e1f780b9a92048ba716b2..a748da5cdb68db8c364e2e77fac19f5cdf1df875 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -251,9 +251,16 @@ protected:
     }
 
     std::vector<double> const& getIntPtDarcyVelocityLocal(
-        const double t, std::vector<double> const& local_p,
-        std::vector<double> const& local_T, std::vector<double>& cache) const
+        const double t, std::vector<double> const& local_x,
+        std::vector<double>& cache) const
     {
+        std::vector<double> local_p{
+            local_x.data() + pressure_index,
+            local_x.data() + pressure_index + pressure_size};
+        std::vector<double> local_T{
+            local_x.data() + temperature_index,
+            local_x.data() + temperature_index + temperature_size};
+
         auto const n_integration_points =
             _integration_method.getNumberOfPoints();
 
diff --git a/ProcessLib/HT/MonolithicHTFEM.h b/ProcessLib/HT/MonolithicHTFEM.h
index 67423f05ff7986bbdb0b4538222747ae705a9520..d1b5db88c6fc0c93a2954945b312bdbb4b7e01e8 100644
--- a/ProcessLib/HT/MonolithicHTFEM.h
+++ b/ProcessLib/HT/MonolithicHTFEM.h
@@ -211,15 +211,9 @@ public:
         auto const indices =
             NumLib::getIndices(this->_element.getID(), dof_table);
         assert(!indices.empty());
-        auto local_x = current_solution.get(indices);
+        auto const& local_x = current_solution.get(indices);
 
-        std::vector<double> local_p(
-            std::make_move_iterator(local_x.begin() + local_x.size() / 2),
-            std::make_move_iterator(local_x.end()));
-        // only T is kept in local_x
-        local_x.erase(local_x.begin() + local_x.size() / 2, local_x.end());
-
-        return this->getIntPtDarcyVelocityLocal(t, local_p, local_x, cache);
+        return this->getIntPtDarcyVelocityLocal(t, local_x, cache);
     }
 
 private:
diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h
index b5f4865a33f54ddaa38c71b1e0813df8e167d1a4..71a10acc432809c5edc57638f8943bb600dd9313 100644
--- a/ProcessLib/HT/StaggeredHTFEM-impl.h
+++ b/ProcessLib/HT/StaggeredHTFEM-impl.h
@@ -51,23 +51,25 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
                               std::vector<double>& local_b_data,
                               LocalCoupledSolutions const& coupled_xs)
 {
-    auto const& local_p = coupled_xs.local_coupled_xs[_hydraulic_process_id];
-    auto const local_matrix_size = local_p.size();
-    // This assertion is valid only if all nodal d.o.f. use the same shape
-    // matrices.
-    assert(local_matrix_size == ShapeFunction::NPOINTS);
-
-    auto const& local_T1 =
-        coupled_xs.local_coupled_xs[_heat_transport_process_id];
-    auto const& local_T0 =
-        coupled_xs.local_coupled_xs0[_heat_transport_process_id];
+    auto const local_p = Eigen::Map<
+        typename ShapeMatricesType::template VectorType<pressure_size> const>(
+        &coupled_xs.local_coupled_xs[pressure_index], pressure_size);
+
+    auto const local_T1 =
+        Eigen::Map<typename ShapeMatricesType::template VectorType<
+            temperature_size> const>(
+            &coupled_xs.local_coupled_xs[temperature_index], temperature_size);
+    auto const local_T0 =
+        Eigen::Map<typename ShapeMatricesType::template VectorType<
+            temperature_size> const>(
+            &coupled_xs.local_coupled_xs0[temperature_index], temperature_size);
 
     auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
-        local_M_data, local_matrix_size, local_matrix_size);
+        local_M_data, pressure_size, pressure_size);
     auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
-        local_K_data, local_matrix_size, local_matrix_size);
-    auto local_b = MathLib::createZeroedVector<LocalVectorType>(
-        local_b_data, local_matrix_size);
+        local_K_data, pressure_size, pressure_size);
+    auto local_b = MathLib::createZeroedVector<LocalVectorType>(local_b_data,
+                                                                pressure_size);
 
     ParameterLib::SpatialPosition pos;
     pos.setElementID(this->_element.getID());
@@ -186,22 +188,19 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
                                   std::vector<double>& /*local_b_data*/,
                                   LocalCoupledSolutions const& coupled_xs)
 {
-    auto const& local_p = coupled_xs.local_coupled_xs[_hydraulic_process_id];
-    auto const local_matrix_size = local_p.size();
-    // This assertion is valid only if all nodal d.o.f. use the same shape
-    // matrices.
-    assert(local_matrix_size == ShapeFunction::NPOINTS);
-
-    auto local_p_Eigen_type =
-        Eigen::Map<const NodalVectorType>(&local_p[0], local_matrix_size);
+    auto const local_p = Eigen::Map<
+        typename ShapeMatricesType::template VectorType<pressure_size> const>(
+        &coupled_xs.local_coupled_xs[pressure_index], pressure_size);
 
-    auto const& local_T1 =
-        coupled_xs.local_coupled_xs[_heat_transport_process_id];
+    auto const local_T1 =
+        Eigen::Map<typename ShapeMatricesType::template VectorType<
+            temperature_size> const>(
+            &coupled_xs.local_coupled_xs[temperature_index], temperature_size);
 
     auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
-        local_M_data, local_matrix_size, local_matrix_size);
+        local_M_data, temperature_size, temperature_size);
     auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
-        local_K_data, local_matrix_size, local_matrix_size);
+        local_K_data, temperature_size, temperature_size);
 
     ParameterLib::SpatialPosition pos;
     pos.setElementID(this->_element.getID());
@@ -274,9 +273,9 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
             intrinsic_permeability / viscosity;
         GlobalDimVectorType const velocity =
             process_data.has_gravity
-                ? GlobalDimVectorType(-K_over_mu * (dNdx * local_p_Eigen_type -
-                                                    fluid_density * b))
-                : GlobalDimVectorType(-K_over_mu * dNdx * local_p_Eigen_type);
+                ? GlobalDimVectorType(-K_over_mu *
+                                      (dNdx * local_p - fluid_density * b))
+                : GlobalDimVectorType(-K_over_mu * dNdx * local_p);
 
         GlobalDimMatrixType const thermal_conductivity_dispersivity =
             this->getThermalConductivityDispersivity(
@@ -303,12 +302,10 @@ StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
     assert(!indices.empty());
     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);
+    auto const local_xs = getCoupledLocalSolutions(
+        this->_coupled_solutions->coupled_xs, indices_of_all_coupled_processes);
 
-    return this->getIntPtDarcyVelocityLocal(
-        t, local_xs[_hydraulic_process_id],
-        local_xs[_heat_transport_process_id], cache);
+    return this->getIntPtDarcyVelocityLocal(t, local_xs, cache);
 }
 }  // namespace HT
 }  // namespace ProcessLib
diff --git a/ProcessLib/HT/StaggeredHTFEM.h b/ProcessLib/HT/StaggeredHTFEM.h
index b1aa5d15fcb169a8f8c8676e10e14321bbcd25d1..76a84c5a1be529ecc06baef9730c0184c58834a3 100644
--- a/ProcessLib/HT/StaggeredHTFEM.h
+++ b/ProcessLib/HT/StaggeredHTFEM.h
@@ -50,6 +50,11 @@ class StaggeredHTFEM : public HTFEM<ShapeFunction, IntegrationMethod, GlobalDim>
         typename ShapeMatricesType::GlobalDimNodalMatrixType;
     using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType;
 
+    using HTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::pressure_index;
+    using HTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::pressure_size;
+    using HTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::temperature_index;
+    using HTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::temperature_size;
+
 public:
     StaggeredHTFEM(MeshLib::Element const& element,
                    std::size_t const local_matrix_size,
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
index 49372e36f3ac102deb467ac6a138770c1219ca22..5102fd192d5409888b3738016e3abccf9ce5f335 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
@@ -357,25 +357,24 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         const LocalCoupledSolutions& local_coupled_solutions)
 {
-    auto const& local_p = local_coupled_solutions.local_coupled_xs[0];
-    auto const& local_u = local_coupled_solutions.local_coupled_xs[1];
-    assert(local_p.size() == pressure_size);
-    assert(local_u.size() == displacement_size);
-
-    auto const local_matrix_size = local_p.size();
     auto local_rhs =
         MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
                                         template VectorType<pressure_size>>(
-            local_b_data, local_matrix_size);
+            local_b_data, pressure_size);
 
     ParameterLib::SpatialPosition pos;
     pos.setElementID(this->_element.getID());
 
-    auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
-        pressure_size> const>(local_p.data(), pressure_size);
-    auto u =
+    auto const p =
+        Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+            pressure_size> const>(
+            &local_coupled_solutions.local_coupled_xs[pressure_index],
+            pressure_size);
+    auto const u =
         Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
-            displacement_size> const>(local_u.data(), displacement_size);
+            displacement_size> const>(
+            &local_coupled_solutions.local_coupled_xs[displacement_index],
+            displacement_size);
 
     auto p_dot =
         Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
@@ -469,14 +468,16 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         const LocalCoupledSolutions& local_coupled_solutions)
 {
-    auto const& local_p = local_coupled_solutions.local_coupled_xs[0];
-    auto const& local_u = local_coupled_solutions.local_coupled_xs[1];
-    assert(local_p.size() == pressure_size);
-    assert(local_u.size() == displacement_size);
-
-    auto u =
+    auto const p =
+        Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+            pressure_size> const>(
+            &local_coupled_solutions.local_coupled_xs[pressure_index],
+            pressure_size);
+    auto const u =
         Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
-            displacement_size> const>(local_u.data(), displacement_size);
+            displacement_size> const>(
+            &local_coupled_solutions.local_coupled_xs[displacement_index],
+            displacement_size);
 
     auto local_Jac = MathLib::createZeroedMatrix<
         typename ShapeMatricesTypeDisplacement::template MatrixType<
@@ -534,7 +535,7 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
         local_Jac.noalias() += B.transpose() * C * B * w;
 
         double p_at_xi = 0.;
-        NumLib::shapeFunctionInterpolate(local_p, N_p, p_at_xi);
+        NumLib::shapeFunctionInterpolate(p, N_p, p_at_xi);
 
         double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
         local_rhs.noalias() -=
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
index 5098caa42eac07306e43f22fdf3fbbe641a197f4..2441d8c21cc777b7948dc3fc5ccc3fe58d02ce87 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
@@ -55,31 +55,25 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions)
 {
-    using DeformationVector =
-        typename ShapeMatricesType::template VectorType<displacement_size>;
     using DeformationMatrix =
         typename ShapeMatricesType::template MatrixType<displacement_size,
                                                         displacement_size>;
-    using PhaseFieldVector =
-        typename ShapeMatricesType::template VectorType<phasefield_size>;
 
-    auto const& local_u = local_coupled_solutions.local_coupled_xs[0];
-    auto const& local_d = local_coupled_solutions.local_coupled_xs[1];
-    assert(local_u.size() == displacement_size);
-    assert(local_d.size() == phasefield_size);
+    assert(local_coupled_solutions.local_coupled_xs.size() ==
+           phasefield_size + displacement_size);
 
-    auto const local_matrix_size = local_u.size();
-    auto d =
-        Eigen::Map<PhaseFieldVector const>(local_d.data(), phasefield_size);
-
-    auto u =
-        Eigen::Map<DeformationVector const>(local_u.data(), displacement_size);
+    auto const d = Eigen::Map<PhaseFieldVector const>(
+        &local_coupled_solutions.local_coupled_xs[phasefield_index],
+        phasefield_size);
+    auto const u = Eigen::Map<DeformationVector const>(
+        &local_coupled_solutions.local_coupled_xs[displacement_index],
+        displacement_size);
 
     auto local_Jac = MathLib::createZeroedMatrix<DeformationMatrix>(
-        local_Jac_data, local_matrix_size, local_matrix_size);
+        local_Jac_data, displacement_size, displacement_size);
 
     auto local_rhs = MathLib::createZeroedVector<DeformationVector>(
-        local_b_data, local_matrix_size);
+        local_b_data, displacement_size);
 
     ParameterLib::SpatialPosition x_position;
     x_position.setElementID(_element.getID());
@@ -158,29 +152,17 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions)
 {
-    using DeformationVector =
-        typename ShapeMatricesType::template VectorType<displacement_size>;
-    using PhaseFieldVector =
-        typename ShapeMatricesType::template VectorType<phasefield_size>;
-    using PhaseFieldMatrix =
-        typename ShapeMatricesType::template MatrixType<phasefield_size,
-                                                        phasefield_size>;
-
-    auto const& local_u = local_coupled_solutions.local_coupled_xs[0];
-    auto const& local_d = local_coupled_solutions.local_coupled_xs[1];
-    assert(local_u.size() == displacement_size);
-    assert(local_d.size() == phasefield_size);
-
-    auto const local_matrix_size = local_d.size();
-    auto d =
-        Eigen::Map<PhaseFieldVector const>(local_d.data(), phasefield_size);
-    auto u =
-        Eigen::Map<DeformationVector const>(local_u.data(), displacement_size);
+    auto const d = Eigen::Map<PhaseFieldVector const>(
+        &local_coupled_solutions.local_coupled_xs[phasefield_index],
+        phasefield_size);
+    auto const u = Eigen::Map<DeformationVector const>(
+        &local_coupled_solutions.local_coupled_xs[displacement_index],
+        displacement_size);
 
     auto local_Jac = MathLib::createZeroedMatrix<PhaseFieldMatrix>(
-        local_Jac_data, local_matrix_size, local_matrix_size);
+        local_Jac_data, phasefield_size, phasefield_size);
     auto local_rhs = MathLib::createZeroedVector<PhaseFieldVector>(
-        local_b_data, local_matrix_size);
+        local_b_data, phasefield_size);
 
     ParameterLib::SpatialPosition x_position;
     x_position.setElementID(_element.getID());
@@ -280,21 +262,13 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
                    });
 
     auto local_coupled_xs =
-        getCurrentLocalSolutions(*cpl_xs, indices_of_processes);
-    assert(local_coupled_xs.size() == 2);
-
-    auto const& local_u = local_coupled_xs[0];
-    auto const& local_d = local_coupled_xs[1];
-
-    assert(local_u.size() == displacement_size);
-    assert(local_d.size() == phasefield_size);
-
-    auto d = Eigen::Map<
-        typename ShapeMatricesType::template VectorType<phasefield_size> const>(
-        local_d.data(), phasefield_size);
+        getCoupledLocalSolutions(cpl_xs->coupled_xs, indices_of_processes);
+    assert(local_coupled_xs.size() == displacement_size + phasefield_size);
 
-    auto u = Eigen::Map<typename ShapeMatricesType::template VectorType<
-        displacement_size> const>(local_u.data(), displacement_size);
+    auto const d = Eigen::Map<PhaseFieldVector const>(
+        &local_coupled_xs[phasefield_index], phasefield_size);
+    auto const u = Eigen::Map<DeformationVector const>(
+        &local_coupled_xs[displacement_index], displacement_size);
 
     ParameterLib::SpatialPosition x_position;
     x_position.setElementID(_element.getID());
@@ -346,22 +320,14 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
                        return NumLib::getIndices(mesh_item_id, dof_table);
                    });
 
-    auto local_coupled_xs =
-        getCurrentLocalSolutions(*cpl_xs, indices_of_processes);
-    assert(local_coupled_xs.size() == 2);
-
-    auto const& local_u = local_coupled_xs[0];
-    auto const& local_d = local_coupled_xs[1];
-
-    assert(local_u.size() == displacement_size);
-    assert(local_d.size() == phasefield_size);
-
-    auto d = Eigen::Map<
-        typename ShapeMatricesType::template VectorType<phasefield_size> const>(
-        local_d.data(), phasefield_size);
+    auto const local_coupled_xs =
+        getCoupledLocalSolutions(cpl_xs->coupled_xs, indices_of_processes);
+    assert(local_coupled_xs.size() == displacement_size + phasefield_size);
 
-    auto u = Eigen::Map<typename ShapeMatricesType::template VectorType<
-        displacement_size> const>(local_u.data(), displacement_size);
+    auto const d = Eigen::Map<PhaseFieldVector const>(
+        &local_coupled_xs[phasefield_index], phasefield_size);
+    auto const u = Eigen::Map<DeformationVector const>(
+        &local_coupled_xs[displacement_index], displacement_size);
 
     ParameterLib::SpatialPosition x_position;
     x_position.setElementID(_element.getID());
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM.h b/ProcessLib/PhaseField/PhaseFieldFEM.h
index 9727fd3346cc9c50574c0cea406e0dcdfa2d0123..14bd3e546c1bac9d6c3905a215a76a0a6bcf501a 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM.h
@@ -107,6 +107,14 @@ template <typename ShapeFunction, typename IntegrationMethod,
           int DisplacementDim>
 class PhaseFieldLocalAssembler : public PhaseFieldLocalAssemblerInterface
 {
+private:
+    static constexpr int displacement_index = 0;
+    static constexpr int displacement_size =
+        ShapeFunction::NPOINTS * DisplacementDim;
+    static constexpr int phasefield_index =
+        displacement_index + displacement_size;
+    static constexpr int phasefield_size = ShapeFunction::NPOINTS;
+
 public:
     using ShapeMatricesType =
         ShapeMatrixPolicyType<ShapeFunction, DisplacementDim>;
@@ -117,6 +125,14 @@ public:
 
     using NodalForceVectorType = typename BMatricesType::NodalForceVectorType;
 
+    using DeformationVector =
+        typename ShapeMatricesType::template VectorType<displacement_size>;
+    using PhaseFieldVector =
+        typename ShapeMatricesType::template VectorType<phasefield_size>;
+    using PhaseFieldMatrix =
+        typename ShapeMatricesType::template MatrixType<phasefield_size,
+                                                        phasefield_size>;
+
     PhaseFieldLocalAssembler(PhaseFieldLocalAssembler const&) = delete;
     PhaseFieldLocalAssembler(PhaseFieldLocalAssembler&&) = delete;
 
@@ -337,12 +353,6 @@ private:
     MeshLib::Element const& _element;
     SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
     bool const _is_axially_symmetric;
-
-    static const int phasefield_index = 0;
-    static const int phasefield_size = ShapeFunction::NPOINTS;
-    static const int displacement_index = ShapeFunction::NPOINTS;
-    static const int displacement_size =
-        ShapeFunction::NPOINTS * DisplacementDim;
 };
 
 }  // namespace PhaseField
diff --git a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM-impl.h b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM-impl.h
index 5e2cdf6b12a51a70e3de79cfb95abeae997c0530..32584dc29b4bd08ace4fc7f2ed3cf1e3abb069b3 100644
--- a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM-impl.h
+++ b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM-impl.h
@@ -64,25 +64,18 @@ void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions)
 {
-    auto const& local_d =
-        local_coupled_solutions.local_coupled_xs[_phase_field_process_id];
-    auto const& local_u =
-        local_coupled_solutions.local_coupled_xs[_mechanics_related_process_id];
-    auto const& local_T =
-        local_coupled_solutions.local_coupled_xs[_heat_conduction_process_id];
-    assert(local_T.size() == temperature_size);
-    assert(local_d.size() == phasefield_size);
-    assert(local_u.size() == displacement_size);
-
-    auto d = Eigen::Map<
-        typename ShapeMatricesType::template VectorType<phasefield_size> const>(
-        local_d.data(), phasefield_size);
-
-    auto u = Eigen::Map<typename ShapeMatricesType::template VectorType<
-        displacement_size> const>(local_u.data(), displacement_size);
-
-    auto T = Eigen::Map<typename ShapeMatricesType::template VectorType<
-        temperature_size> const>(local_T.data(), temperature_size);
+    assert(local_coupled_solutions.local_coupled_xs.size() ==
+           phasefield_size + displacement_size + temperature_size);
+
+    auto const d = Eigen::Map<PhaseFieldVector const>(
+        &local_coupled_solutions.local_coupled_xs[phasefield_index],
+        phasefield_size);
+    auto const u = Eigen::Map<DeformationVector const>(
+        &local_coupled_solutions.local_coupled_xs[displacement_index],
+        displacement_size);
+    auto const T = Eigen::Map<TemperatureVector const>(
+        &local_coupled_solutions.local_coupled_xs[temperature_index],
+        temperature_size);
 
     auto local_Jac = MathLib::createZeroedMatrix<
         typename ShapeMatricesType::template MatrixType<displacement_size,
@@ -172,19 +165,15 @@ void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions)
 {
-    auto const& local_d =
-        local_coupled_solutions.local_coupled_xs[_phase_field_process_id];
-    auto const& local_T =
-        local_coupled_solutions.local_coupled_xs[_heat_conduction_process_id];
-    assert(local_T.size() == temperature_size);
-    assert(local_d.size() == phasefield_size);
+    assert(local_coupled_solutions.local_coupled_xs.size() ==
+           phasefield_size + displacement_size + temperature_size);
 
-    auto d = Eigen::Map<
-        typename ShapeMatricesType::template VectorType<phasefield_size> const>(
-        local_d.data(), phasefield_size);
-
-    auto T = Eigen::Map<typename ShapeMatricesType::template VectorType<
-        temperature_size> const>(local_T.data(), temperature_size);
+    auto const d = Eigen::Map<PhaseFieldVector const>(
+        &local_coupled_solutions.local_coupled_xs[phasefield_index],
+        phasefield_size);
+    auto const T = Eigen::Map<TemperatureVector const>(
+        &local_coupled_solutions.local_coupled_xs[temperature_index],
+        temperature_size);
 
     auto T_dot = Eigen::Map<typename ShapeMatricesType::template VectorType<
         temperature_size> const>(local_xdot.data(), temperature_size);
@@ -276,13 +265,12 @@ void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions)
 {
-    auto const& local_d =
-        local_coupled_solutions.local_coupled_xs[_phase_field_process_id];
-    assert(local_d.size() == phasefield_size);
+    assert(local_coupled_solutions.local_coupled_xs.size() ==
+           phasefield_size + displacement_size + temperature_size);
 
-    auto d = Eigen::Map<
-        typename ShapeMatricesType::template VectorType<phasefield_size> const>(
-        local_d.data(), phasefield_size);
+    auto const d = Eigen::Map<PhaseFieldVector const>(
+        &local_coupled_solutions.local_coupled_xs[phasefield_index],
+        phasefield_size);
 
     auto local_Jac = MathLib::createZeroedMatrix<
         typename ShapeMatricesType::template MatrixType<phasefield_size,
diff --git a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
index fff7bf8ecafb28171437546d242b1489d4cdd246..e3dcf0e5bcfbe9b59d45c78e275e0f1d42b01049 100644
--- a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
+++ b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
@@ -111,6 +111,17 @@ template <typename ShapeFunction, typename IntegrationMethod,
 class ThermoMechanicalPhaseFieldLocalAssembler
     : public ThermoMechanicalPhaseFieldLocalAssemblerInterface
 {
+private:
+    static constexpr int temperature_index = 0;
+    static constexpr int temperature_size = ShapeFunction::NPOINTS;
+    static constexpr int displacement_index =
+        temperature_index + temperature_size;
+    static constexpr int displacement_size =
+        ShapeFunction::NPOINTS * DisplacementDim;
+    static constexpr int phasefield_index =
+        displacement_index + displacement_size;
+    static constexpr int phasefield_size = ShapeFunction::NPOINTS;
+
 public:
     using ShapeMatricesType =
         ShapeMatrixPolicyType<ShapeFunction, DisplacementDim>;
@@ -123,6 +134,13 @@ public:
 
     using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
 
+    using TemperatureVector =
+        typename ShapeMatricesType::template VectorType<temperature_size>;
+    using DeformationVector =
+        typename ShapeMatricesType::template VectorType<displacement_size>;
+    using PhaseFieldVector =
+        typename ShapeMatricesType::template VectorType<phasefield_size>;
+
     ThermoMechanicalPhaseFieldLocalAssembler(
         ThermoMechanicalPhaseFieldLocalAssembler const&) = delete;
     ThermoMechanicalPhaseFieldLocalAssembler(
@@ -371,14 +389,6 @@ private:
     SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
     bool const _is_axially_symmetric;
 
-    static const int temperature_index = 0;
-    static const int temperature_size = ShapeFunction::NPOINTS;
-    static const int phasefield_index = ShapeFunction::NPOINTS;
-    static const int phasefield_size = ShapeFunction::NPOINTS;
-    static const int displacement_index = 2 * ShapeFunction::NPOINTS;
-    static const int displacement_size =
-        ShapeFunction::NPOINTS * DisplacementDim;
-
     /// ID of the processes that contains mechanical process.
     int const _mechanics_related_process_id;
 
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h
index ac1af0754c992afd31085af45e7e1a015d074678..7e51944ba7538578f21e7df462cf6dc2731dc15e 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h
@@ -338,26 +338,22 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         const LocalCoupledSolutions& local_coupled_solutions)
 {
-    auto const& local_T_vector =
-        local_coupled_solutions
-            .local_coupled_xs[_process_data.heat_conduction_process_id];
-    assert(local_T_vector.size() == temperature_size);
     auto const local_T =
         Eigen::Map<typename ShapeMatricesType::template VectorType<
-            temperature_size> const>(local_T_vector.data(), temperature_size);
+            temperature_size> const>(
+            &local_coupled_solutions.local_coupled_xs[temperature_index],
+            temperature_size);
 
-    auto const& local_T0_vector = local_coupled_solutions.local_coupled_xs0[0];
-    assert(local_T0_vector.size() == temperature_size);
     auto const local_T0 =
         Eigen::Map<typename ShapeMatricesType::template VectorType<
-            temperature_size> const>(local_T0_vector.data(), temperature_size);
+            temperature_size> const>(
+            &local_coupled_solutions.local_coupled_xs0[temperature_index],
+            temperature_size);
 
-    auto const& local_u_vector =
-        local_coupled_solutions
-            .local_coupled_xs[_process_data.mechanics_process_id];
-    assert(local_u_vector.size() == displacement_size);
     auto const u = Eigen::Map<typename ShapeMatricesType::template VectorType<
-        displacement_size> const>(local_u_vector.data(), displacement_size);
+        displacement_size> const>(
+        &local_coupled_solutions.local_coupled_xs[displacement_index],
+        displacement_size);
 
     auto local_Jac = MathLib::createZeroedMatrix<
         typename ShapeMatricesType::template MatrixType<displacement_size,
@@ -476,20 +472,19 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         const LocalCoupledSolutions& local_coupled_solutions)
 {
-    auto const& local_T_vector =
-        local_coupled_solutions
-            .local_coupled_xs[_process_data.heat_conduction_process_id];
-    assert(local_T_vector.size() == temperature_size);
     auto const local_T =
         Eigen::Map<typename ShapeMatricesType::template VectorType<
-            temperature_size> const>(local_T_vector.data(), temperature_size);
+            temperature_size> const>(
+            &local_coupled_solutions.local_coupled_xs[temperature_index],
+            temperature_size);
 
-    auto const& local_T0_vector = local_coupled_solutions.local_coupled_xs0[0];
-    assert(local_T0_vector.size() == temperature_size);
-    auto const local_dT =
-        local_T -
+    auto const local_T0 =
         Eigen::Map<typename ShapeMatricesType::template VectorType<
-            temperature_size> const>(local_T0_vector.data(), temperature_size);
+            temperature_size> const>(
+            &local_coupled_solutions.local_coupled_xs0[temperature_index],
+            temperature_size);
+
+    auto const local_dT = local_T - local_T0;
 
     auto local_Jac = MathLib::createZeroedMatrix<
         typename ShapeMatricesType::template MatrixType<temperature_size,
diff --git a/ProcessLib/VectorMatrixAssembler.cpp b/ProcessLib/VectorMatrixAssembler.cpp
index 68268ff6096c3ab5f3eb21e723216c4bf6612716..1a05594abcf5dabc84b18a9edf4788c9b3c95d3d 100644
--- a/ProcessLib/VectorMatrixAssembler.cpp
+++ b/ProcessLib/VectorMatrixAssembler.cpp
@@ -68,10 +68,10 @@ void VectorMatrixAssembler::assemble(
     }
     else
     {
-        auto local_coupled_xs0 =
-            getPreviousLocalSolutions(*cpl_xs, indices_of_processes);
+        auto local_coupled_xs0 = getCoupledLocalSolutions(cpl_xs->coupled_xs_t0,
+                                                          indices_of_processes);
         auto local_coupled_xs =
-            getCurrentLocalSolutions(*cpl_xs, indices_of_processes);
+            getCoupledLocalSolutions(cpl_xs->coupled_xs, indices_of_processes);
 
         ProcessLib::LocalCoupledSolutions local_coupled_solutions(
             std::move(local_coupled_xs0), std::move(local_coupled_xs));
@@ -136,10 +136,10 @@ void VectorMatrixAssembler::assembleWithJacobian(
     }
     else
     {
-        auto local_coupled_xs0 =
-            getPreviousLocalSolutions(*cpl_xs, indices_of_processes);
+        auto local_coupled_xs0 = getCoupledLocalSolutions(cpl_xs->coupled_xs_t0,
+                                                          indices_of_processes);
         auto local_coupled_xs =
-            getCurrentLocalSolutions(*cpl_xs, indices_of_processes);
+            getCoupledLocalSolutions(cpl_xs->coupled_xs, indices_of_processes);
 
         ProcessLib::LocalCoupledSolutions local_coupled_solutions(
             std::move(local_coupled_xs0), std::move(local_coupled_xs));