diff --git a/ProcessLib/AbstractJacobianAssembler.h b/ProcessLib/AbstractJacobianAssembler.h
index e11cb1917c75945caa7930cad3ba2dc7c8609ea5..b148e7cf952fe7f354eb1088b2f59aea268bf814 100644
--- a/ProcessLib/AbstractJacobianAssembler.h
+++ b/ProcessLib/AbstractJacobianAssembler.h
@@ -10,7 +10,9 @@
 
 #pragma once
 
+#include <Eigen/Dense>
 #include <vector>
+
 #include "BaseLib/Error.h"
 
 namespace ProcessLib
@@ -38,9 +40,10 @@ public:
     //! \f$b\f$ with coupling.
     virtual void assembleWithJacobianForStaggeredScheme(
         LocalAssemblerInterface& /*local_assembler*/, double const /*t*/,
-        double const /*dt*/, std::vector<double> const& /*local_xdot*/,
-        const double /*dxdot_dx*/, const double /*dx_dx*/,
-        int const /*process_id*/, std::vector<double>& /*local_M_data*/,
+        double const /*dt*/, Eigen::VectorXd const& /*local_x*/,
+        std::vector<double> const& /*local_xdot*/, const double /*dxdot_dx*/,
+        const double /*dx_dx*/, int const /*process_id*/,
+        std::vector<double>& /*local_M_data*/,
         std::vector<double>& /*local_K_data*/,
         std::vector<double>& /*local_b_data*/,
         std::vector<double>& /*local_Jac_data*/,
diff --git a/ProcessLib/AnalyticalJacobianAssembler.cpp b/ProcessLib/AnalyticalJacobianAssembler.cpp
index 1959c2ef151139493faf15ef3ea8516982fd1c06..991f5b8edc305f9cf201ec9d0a4009c5799d2ba7 100644
--- a/ProcessLib/AnalyticalJacobianAssembler.cpp
+++ b/ProcessLib/AnalyticalJacobianAssembler.cpp
@@ -28,14 +28,14 @@ void AnalyticalJacobianAssembler::assembleWithJacobian(
 
 void AnalyticalJacobianAssembler::assembleWithJacobianForStaggeredScheme(
     LocalAssemblerInterface& local_assembler, double const t, double const dt,
-    std::vector<double> const& local_xdot, const double dxdot_dx,
-    const double dx_dx, int const process_id, std::vector<double>& local_M_data,
-    std::vector<double>& local_K_data, std::vector<double>& local_b_data,
-    std::vector<double>& local_Jac_data,
+    Eigen::VectorXd const& local_x, std::vector<double> const& local_xdot,
+    const double dxdot_dx, const double dx_dx, int const process_id,
+    std::vector<double>& local_M_data, std::vector<double>& local_K_data,
+    std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
     LocalCoupledSolutions const& local_coupled_solutions)
 {
     local_assembler.assembleWithJacobianForStaggeredScheme(
-        t, dt, local_xdot, dxdot_dx, dx_dx, process_id, local_M_data,
+        t, dt, local_x, local_xdot, dxdot_dx, dx_dx, process_id, local_M_data,
         local_K_data, local_b_data, local_Jac_data, local_coupled_solutions);
 }
 
diff --git a/ProcessLib/AnalyticalJacobianAssembler.h b/ProcessLib/AnalyticalJacobianAssembler.h
index 814f1166f0937c7d084005b786f9b5878bc70095..c8d7feef0a402fd18e41a368100549e7e9aeae0b 100644
--- a/ProcessLib/AnalyticalJacobianAssembler.h
+++ b/ProcessLib/AnalyticalJacobianAssembler.h
@@ -43,8 +43,9 @@ public:
 
     void assembleWithJacobianForStaggeredScheme(
         LocalAssemblerInterface& local_assembler, double const t,
-        double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx, int const process_id,
+        double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, int const process_id,
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions) override;
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index cce5776b2852c4b5350b50b147d7d65b2afaf25d..8f461ad47a1303d8e5d949720243edf57197726d 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -397,20 +397,20 @@ public:
     }
 
     void assembleForStaggeredScheme(
-        double const t, double const dt, int const process_id,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
-        std::vector<double>& local_b_data,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        int const process_id, std::vector<double>& local_M_data,
+        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
         LocalCoupledSolutions const& coupled_xs) override
     {
         if (process_id == hydraulic_process_id)
         {
-            assembleHydraulicEquation(t, dt, local_M_data, local_K_data,
-                                      local_b_data, coupled_xs);
+            assembleHydraulicEquation(t, dt, local_x, local_M_data,
+                                      local_K_data, local_b_data, coupled_xs);
         }
         else
         {
             // Go for assembling in an order of transport process id.
-            assembleComponentTransportEquation(t, dt, local_M_data,
+            assembleComponentTransportEquation(t, dt, local_x, local_M_data,
                                                local_K_data, local_b_data,
                                                coupled_xs, process_id);
         }
@@ -418,16 +418,15 @@ public:
 
     void assembleHydraulicEquation(double const t,
                                    double const dt,
+                                   Eigen::VectorXd const& local_x,
                                    std::vector<double>& local_M_data,
                                    std::vector<double>& local_K_data,
                                    std::vector<double>& local_b_data,
                                    LocalCoupledSolutions const& coupled_xs)
     {
-        auto local_p = Eigen::Map<const NodalVectorType>(
-            &coupled_xs.local_coupled_xs[pressure_index], pressure_size);
-        auto local_C = Eigen::Map<const NodalVectorType>(
-            &coupled_xs.local_coupled_xs[first_concentration_index],
-            concentration_size);
+        auto local_p = local_x.template segment<pressure_size>(pressure_index);
+        auto local_C = local_x.template segment<concentration_size>(
+            first_concentration_index);
         auto local_C0 = Eigen::Map<const NodalVectorType>(
             &coupled_xs.local_coupled_xs0[first_concentration_index],
             concentration_size);
@@ -530,18 +529,15 @@ public:
     }
 
     void assembleComponentTransportEquation(
-        double const t, double const dt, std::vector<double>& local_M_data,
-        std::vector<double>& local_K_data,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& /*local_b_data*/,
         LocalCoupledSolutions const& coupled_xs, int const transport_process_id)
     {
-        auto local_C = Eigen::Map<const NodalVectorType>(
-            &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[pressure_index], pressure_size);
+        auto local_p = local_x.template segment<pressure_size>(pressure_index);
+        auto local_C = local_x.template segment<concentration_size>(
+            first_concentration_index +
+            (transport_process_id - 1) * concentration_size);
         auto local_p0 = Eigen::Map<const NodalVectorType>(
             &coupled_xs.local_coupled_xs0[pressure_index], pressure_size);
 
diff --git a/ProcessLib/CoupledSolutionsForStaggeredScheme.h b/ProcessLib/CoupledSolutionsForStaggeredScheme.h
index cfce3363a5ecb74622f48c0d9acc458e5841126d..a799cecacd98bc2a83d16f4bb28287f01699f830 100644
--- a/ProcessLib/CoupledSolutionsForStaggeredScheme.h
+++ b/ProcessLib/CoupledSolutionsForStaggeredScheme.h
@@ -49,17 +49,13 @@ struct CoupledSolutionsForStaggeredScheme
  */
 struct LocalCoupledSolutions
 {
-    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_))
+    LocalCoupledSolutions(std::vector<double>&& local_coupled_xs0_)
+        : local_coupled_xs0(std::move(local_coupled_xs0_))
     {
     }
 
     /// Local solutions of the previous time step.
     std::vector<double> const local_coupled_xs0;
-    /// Local solutions of the current time step.
-    std::vector<double> const local_coupled_xs;
 };
 
 /**
diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h
index fb6f271c2a9f2089dd005cc25eef01829f210363..06d557212e3750b69e4630d2542568d65f4bb76b 100644
--- a/ProcessLib/HT/StaggeredHTFEM-impl.h
+++ b/ProcessLib/HT/StaggeredHTFEM-impl.h
@@ -25,6 +25,7 @@ template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
     assembleForStaggeredScheme(double const t, double const dt,
+                               Eigen::VectorXd const& local_x,
                                int const process_id,
                                std::vector<double>& local_M_data,
                                std::vector<double>& local_K_data,
@@ -33,32 +34,31 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
 {
     if (process_id == _heat_transport_process_id)
     {
-        assembleHeatTransportEquation(t, local_M_data, local_K_data,
-                                      local_b_data, coupled_xs);
+        assembleHeatTransportEquation(t, local_x, local_M_data, local_K_data,
+                                      local_b_data);
         return;
     }
 
-    assembleHydraulicEquation(t, dt, local_M_data, local_K_data, local_b_data,
-                              coupled_xs);
+    assembleHydraulicEquation(t, dt, local_x, local_M_data, local_K_data,
+                              local_b_data, coupled_xs);
 }
 
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
     assembleHydraulicEquation(double const t, double const dt,
+                              Eigen::VectorXd const& local_x,
                               std::vector<double>& local_M_data,
                               std::vector<double>& local_K_data,
                               std::vector<double>& local_b_data,
                               LocalCoupledSolutions const& coupled_xs)
 {
-    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_p =
+        local_x.template segment<pressure_size>(pressure_index);
 
     auto const local_T1 =
-        Eigen::Map<typename ShapeMatricesType::template VectorType<
-            temperature_size> const>(
-            &coupled_xs.local_coupled_xs[temperature_index], temperature_size);
+        local_x.template segment<temperature_size>(temperature_index);
+
     auto const local_T0 =
         Eigen::Map<typename ShapeMatricesType::template VectorType<
             temperature_size> const>(
@@ -183,19 +183,16 @@ template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
     assembleHeatTransportEquation(double const t,
+                                  Eigen::VectorXd const& local_x,
                                   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_b_data*/)
 {
-    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_p =
+        local_x.template segment<pressure_size>(pressure_index);
 
     auto const local_T1 =
-        Eigen::Map<typename ShapeMatricesType::template VectorType<
-            temperature_size> const>(
-            &coupled_xs.local_coupled_xs[temperature_index], temperature_size);
+        local_x.template segment<temperature_size>(temperature_index);
 
     auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
         local_M_data, temperature_size, temperature_size);
diff --git a/ProcessLib/HT/StaggeredHTFEM.h b/ProcessLib/HT/StaggeredHTFEM.h
index 8ce6c973de54acc96d8a08480617c70d8b3b9a26..7253962fd48c897d61a79a56783013f2dbaef500 100644
--- a/ProcessLib/HT/StaggeredHTFEM.h
+++ b/ProcessLib/HT/StaggeredHTFEM.h
@@ -72,9 +72,9 @@ public:
     }
 
     void assembleForStaggeredScheme(
-        double const t, double const dt, int const process_id,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
-        std::vector<double>& local_b_data,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        int const process_id, std::vector<double>& local_M_data,
+        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
         LocalCoupledSolutions const& coupled_xs) override;
 
     std::vector<double> const& getIntPtDarcyVelocity(
@@ -85,15 +85,17 @@ public:
 
 private:
     void assembleHydraulicEquation(double const t, double const dt,
+                                   Eigen::VectorXd const& local_x,
                                    std::vector<double>& local_M_data,
                                    std::vector<double>& local_K_data,
                                    std::vector<double>& local_b_data,
                                    LocalCoupledSolutions const& coupled_xs);
 
-    void assembleHeatTransportEquation(
-        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);
+    void assembleHeatTransportEquation(double const t,
+                                       Eigen::VectorXd const& local_x,
+                                       std::vector<double>& local_M_data,
+                                       std::vector<double>& local_K_data,
+                                       std::vector<double>& local_b_data);
     const int _heat_transport_process_id;
     const int _hydraulic_process_id;
 };
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
index 8e9e177562de606254bfea781776ba26ddc3edd2..a9b1d442f68f92b26a41b0d6dfafdcd1f82b81a6 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
@@ -351,12 +351,11 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
                                   ShapeFunctionPressure, IntegrationMethod,
                                   DisplacementDim>::
     assembleWithJacobianForPressureEquations(
-        const double t, double const dt, const std::vector<double>& local_xdot,
-        const double /*dxdot_dx*/, const double /*dx_dx*/,
-        std::vector<double>& /*local_M_data*/,
+        const double t, double const dt, Eigen::VectorXd const& local_x,
+        const std::vector<double>& local_xdot, const double /*dxdot_dx*/,
+        const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
         std::vector<double>& /*local_K_data*/,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        const LocalCoupledSolutions& local_coupled_solutions)
+        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
 {
     auto local_rhs =
         MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
@@ -366,16 +365,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
     ParameterLib::SpatialPosition pos;
     pos.setElementID(this->_element.getID());
 
-    auto const p =
-        Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
-            pressure_size> const>(
-            &local_coupled_solutions.local_coupled_xs[pressure_index],
-            pressure_size);
+    auto const p = local_x.template segment<pressure_size>(pressure_index);
     auto const u =
-        Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
-            displacement_size> const>(
-            &local_coupled_solutions.local_coupled_xs[displacement_index],
-            displacement_size);
+        local_x.template segment<displacement_size>(displacement_index);
 
     auto p_dot =
         Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
@@ -462,23 +454,15 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
                                   ShapeFunctionPressure, IntegrationMethod,
                                   DisplacementDim>::
     assembleWithJacobianForDeformationEquations(
-        const double t, double const dt,
+        const double t, double const dt, Eigen::VectorXd const& local_x,
         const std::vector<double>& /*local_xdot*/, const double /*dxdot_dx*/,
         const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
         std::vector<double>& /*local_K_data*/,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        const LocalCoupledSolutions& local_coupled_solutions)
+        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
 {
-    auto const p =
-        Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
-            pressure_size> const>(
-            &local_coupled_solutions.local_coupled_xs[pressure_index],
-            pressure_size);
+    auto const p = local_x.template segment<pressure_size>(pressure_index);
     auto const u =
-        Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
-            displacement_size> const>(
-            &local_coupled_solutions.local_coupled_xs[displacement_index],
-            displacement_size);
+        local_x.template segment<displacement_size>(displacement_index);
 
     auto local_Jac = MathLib::createZeroedMatrix<
         typename ShapeMatricesTypeDisplacement::template MatrixType<
@@ -552,25 +536,26 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
                                   ShapeFunctionPressure, IntegrationMethod,
                                   DisplacementDim>::
     assembleWithJacobianForStaggeredScheme(
-        const double t, double const dt, const std::vector<double>& local_xdot,
-        const double dxdot_dx, const double dx_dx, int const process_id,
+        const double t, double const dt, Eigen::VectorXd const& local_x,
+        const std::vector<double>& local_xdot, const double dxdot_dx,
+        const double dx_dx, int const process_id,
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        const LocalCoupledSolutions& local_coupled_solutions)
+        const LocalCoupledSolutions& /*local_coupled_solutions*/)
 {
     // For the equations with pressure
     if (process_id == 0)
     {
         assembleWithJacobianForPressureEquations(
-            t, dt, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
-            local_b_data, local_Jac_data, local_coupled_solutions);
+            t, dt, local_x, local_xdot, dxdot_dx, dx_dx, local_M_data,
+            local_K_data, local_b_data, local_Jac_data);
         return;
     }
 
     // For the equations with deformation
     assembleWithJacobianForDeformationEquations(
-        t, dt, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
-        local_b_data, local_Jac_data, local_coupled_solutions);
+        t, dt, local_x, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
+        local_b_data, local_Jac_data);
 }
 
 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index 217a7689ef45b4a46f7427eec5822145f902e4ec..ed902323a1da058284ad655ba519686a47a059de 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -158,8 +158,9 @@ public:
                               std::vector<double>& local_Jac_data) override;
 
     void assembleWithJacobianForStaggeredScheme(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx, int const process_id,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, int const process_id,
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions) override;
@@ -391,15 +392,13 @@ private:
      * @param local_b_data    Right hand side vector of an element.
      * @param local_Jac_data  Element Jacobian matrix for the Newton-Raphson
      *                        method.
-     * @param local_coupled_solutions Nodal values of solutions of the coupled
-     *                                processes of an element.
      */
     void assembleWithJacobianForDeformationEquations(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        LocalCoupledSolutions const& local_coupled_solutions);
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, std::vector<double>& local_M_data,
+        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
+        std::vector<double>& local_Jac_data);
 
     /**
      * Assemble local matrices and vectors arise from the linearized discretized
@@ -427,15 +426,13 @@ private:
      * @param local_b_data    Right hand side vector of an element.
      * @param local_Jac_data  Element Jacobian matrix for the Newton-Raphson
      *                        method.
-     * @param local_coupled_solutions Nodal values of solutions of the coupled
-     *                                processes of an element.
      */
     void assembleWithJacobianForPressureEquations(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        LocalCoupledSolutions const& local_coupled_solutions);
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, std::vector<double>& local_M_data,
+        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
+        std::vector<double>& local_Jac_data);
 
     unsigned getNumberOfIntegrationPoints() const override;
 
diff --git a/ProcessLib/LocalAssemblerInterface.cpp b/ProcessLib/LocalAssemblerInterface.cpp
index be82ee86765cb793d403a9724c7a2d61eaea65b2..65dee4cc301a18e2b58986b570ca92c5b5678c27 100644
--- a/ProcessLib/LocalAssemblerInterface.cpp
+++ b/ProcessLib/LocalAssemblerInterface.cpp
@@ -28,8 +28,8 @@ void LocalAssemblerInterface::assemble(double const /*t*/,
 }
 
 void LocalAssemblerInterface::assembleForStaggeredScheme(
-    double const /*t*/, double const /*dt*/, int const /*process_id*/,
-    std::vector<double>& /*local_M_data*/,
+    double const /*t*/, double const /*dt*/, Eigen::VectorXd const& /*local_x*/,
+    int const /*process_id*/, std::vector<double>& /*local_M_data*/,
     std::vector<double>& /*local_K_data*/,
     std::vector<double>& /*local_b_data*/,
     LocalCoupledSolutions const& /*coupled_solutions*/)
@@ -54,7 +54,7 @@ void LocalAssemblerInterface::assembleWithJacobian(
 }
 
 void LocalAssemblerInterface::assembleWithJacobianForStaggeredScheme(
-    double const /*t*/, double const /*dt*/,
+    double const /*t*/, double const /*dt*/, Eigen::VectorXd const& /*local_x*/,
     std::vector<double> const& /*local_xdot*/, const double /*dxdot_dx*/,
     const double /*dx_dx*/, int const /*process_id*/,
     std::vector<double>& /*local_M_data*/,
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index 7d447480c6ff3534dae6d8424c8da076bdf15119..b6bdd65e5f0e75747ceff3e4c2d43458b2e5b473 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -53,9 +53,9 @@ public:
                           std::vector<double>& local_b_data);
 
     virtual void assembleForStaggeredScheme(
-        double const t, double const dt, int const process_id,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
-        std::vector<double>& local_b_data,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        int const process_id, std::vector<double>& local_M_data,
+        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
         LocalCoupledSolutions const& coupled_solutions);
 
     virtual void assembleWithJacobian(double const t, double const dt,
@@ -68,8 +68,9 @@ public:
                                       std::vector<double>& local_Jac_data);
 
     virtual void assembleWithJacobianForStaggeredScheme(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx, int const process_id,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, int const process_id,
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions);
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
index 2441d8c21cc777b7948dc3fc5ccc3fe58d02ce87..acba68d7b685cf089416c8e15ec96bfcd511cb2b 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
@@ -22,25 +22,26 @@ template <typename ShapeFunction, typename IntegrationMethod,
 void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
                               DisplacementDim>::
     assembleWithJacobianForStaggeredScheme(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx, int const process_id,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, int const process_id,
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        LocalCoupledSolutions const& local_coupled_solutions)
+        LocalCoupledSolutions const& /*local_coupled_solutions*/)
 {
     // For the equations with phase field.
     if (process_id == 1)
     {
         assembleWithJacobianPhaseFiledEquations(
-            t, dt, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
-            local_b_data, local_Jac_data, local_coupled_solutions);
+            t, dt, local_x, local_xdot, dxdot_dx, dx_dx, local_M_data,
+            local_K_data, local_b_data, local_Jac_data);
         return;
     }
 
     // For the equations with deformation
     assembleWithJacobianForDeformationEquations(
-        t, dt, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
-        local_b_data, local_Jac_data, local_coupled_solutions);
+        t, dt, local_x, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
+        local_b_data, local_Jac_data);
 }
 
 template <typename ShapeFunction, typename IntegrationMethod,
@@ -48,26 +49,19 @@ template <typename ShapeFunction, typename IntegrationMethod,
 void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
                               DisplacementDim>::
     assembleWithJacobianForDeformationEquations(
-        double const t, double const dt,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
         std::vector<double> const& /*local_xdot*/, const double /*dxdot_dx*/,
         const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
         std::vector<double>& /*local_K_data*/,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        LocalCoupledSolutions const& local_coupled_solutions)
+        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
 {
     using DeformationMatrix =
         typename ShapeMatricesType::template MatrixType<displacement_size,
                                                         displacement_size>;
 
-    assert(local_coupled_solutions.local_coupled_xs.size() ==
-           phasefield_size + 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 const d = local_x.template segment<phasefield_size>(phasefield_index);
+    auto const u =
+        local_x.template segment<displacement_size>(displacement_index);
 
     auto local_Jac = MathLib::createZeroedMatrix<DeformationMatrix>(
         local_Jac_data, displacement_size, displacement_size);
@@ -145,19 +139,15 @@ template <typename ShapeFunction, typename IntegrationMethod,
 void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
                               DisplacementDim>::
     assembleWithJacobianPhaseFiledEquations(
-        double const t, double const dt,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
         std::vector<double> const& /*local_xdot*/, const double /*dxdot_dx*/,
         const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
         std::vector<double>& /*local_K_data*/,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        LocalCoupledSolutions const& local_coupled_solutions)
+        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
 {
-    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 d = local_x.template segment<phasefield_size>(phasefield_index);
+    auto const u =
+        local_x.template segment<displacement_size>(displacement_index);
 
     auto local_Jac = MathLib::createZeroedMatrix<PhaseFieldMatrix>(
         local_Jac_data, phasefield_size, phasefield_size);
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM.h b/ProcessLib/PhaseField/PhaseFieldFEM.h
index 4e4ce9cdccbcef00f6f320a5ce4a66acf0c77a3f..7ebc5de7870230c868c6cadfcd87c10cae5589c7 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM.h
@@ -219,8 +219,9 @@ public:
     }
 
     void assembleWithJacobianForStaggeredScheme(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx, int const process_id,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, int const process_id,
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions) override;
@@ -328,18 +329,18 @@ private:
     }
 
     void assembleWithJacobianPhaseFiledEquations(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        LocalCoupledSolutions const& local_coupled_solutions);
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, std::vector<double>& local_M_data,
+        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
+        std::vector<double>& local_Jac_data);
 
     void assembleWithJacobianForDeformationEquations(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        LocalCoupledSolutions const& local_coupled_solutions);
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, std::vector<double>& local_M_data,
+        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
+        std::vector<double>& local_Jac_data);
 
     PhaseFieldProcessData<DisplacementDim>& _process_data;
 
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index d4450f4f295d16d82dca4e385f7ca26ef6c54f35..f3a9250c75647270472672056c2d62fc5905231b 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -741,12 +741,12 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
                                      DisplacementDim>::
     assembleWithJacobianForPressureEquations(
         const double /*t*/, double const /*dt*/,
+        Eigen::VectorXd const& /*local_x*/,
         const std::vector<double>& /*local_xdot*/, const double /*dxdot_dx*/,
         const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
         std::vector<double>& /*local_K_data*/,
         std::vector<double>& /*local_b_data*/,
-        std::vector<double>& /*local_Jac_data*/,
-        const LocalCoupledSolutions& /*local_coupled_solutions*/)
+        std::vector<double>& /*local_Jac_data*/)
 {
     OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
 }
@@ -758,12 +758,12 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
                                      DisplacementDim>::
     assembleWithJacobianForDeformationEquations(
         const double /*t*/, double const /*dt*/,
+        Eigen::VectorXd const& /*local_x*/,
         const std::vector<double>& /*local_xdot*/, const double /*dxdot_dx*/,
         const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
         std::vector<double>& /*local_K_data*/,
         std::vector<double>& /*local_b_data*/,
-        std::vector<double>& /*local_Jac_data*/,
-        const LocalCoupledSolutions& /*local_coupled_solutions*/)
+        std::vector<double>& /*local_Jac_data*/)
 {
     OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
 }
@@ -774,25 +774,26 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
                                      ShapeFunctionPressure, IntegrationMethod,
                                      DisplacementDim>::
     assembleWithJacobianForStaggeredScheme(
-        const double t, double const dt, const std::vector<double>& local_xdot,
-        const double dxdot_dx, const double dx_dx, int const process_id,
+        const double t, double const dt, Eigen::VectorXd const& local_x,
+        const std::vector<double>& local_xdot, const double dxdot_dx,
+        const double dx_dx, int const process_id,
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        const LocalCoupledSolutions& local_coupled_solutions)
+        const LocalCoupledSolutions& /*local_coupled_solutions*/)
 {
     // For the equations with pressure
     if (process_id == 0)
     {
         assembleWithJacobianForPressureEquations(
-            t, dt, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
-            local_b_data, local_Jac_data, local_coupled_solutions);
+            t, dt, local_x, local_xdot, dxdot_dx, dx_dx, local_M_data,
+            local_K_data, local_b_data, local_Jac_data);
         return;
     }
 
     // For the equations with deformation
     assembleWithJacobianForDeformationEquations(
-        t, dt, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
-        local_b_data, local_Jac_data, local_coupled_solutions);
+        t, dt, local_x, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
+        local_b_data, local_Jac_data);
 }
 
 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
index 89ce9cfa5068b6afe727e64e5d9d17e395712217..f93fdd22b2b8115ac527f5311b50d43a13a76fac 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
@@ -85,8 +85,9 @@ public:
                               std::vector<double>& local_Jac_data) override;
 
     void assembleWithJacobianForStaggeredScheme(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx, int const process_id,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, int const process_id,
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions) override;
@@ -199,15 +200,13 @@ private:
      * @param local_b_data    Right hand side vector of an element.
      * @param local_Jac_data  Element Jacobian matrix for the Newton-Raphson
      *                        method.
-     * @param local_coupled_solutions Nodal values of solutions of the coupled
-     *                                processes of an element.
      */
     void assembleWithJacobianForDeformationEquations(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        LocalCoupledSolutions const& local_coupled_solutions);
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, std::vector<double>& local_M_data,
+        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
+        std::vector<double>& local_Jac_data);
 
     /**
      * Assemble local matrices and vectors arise from the linearized discretized
@@ -235,15 +234,13 @@ private:
      * @param local_b_data    Right hand side vector of an element.
      * @param local_Jac_data  Element Jacobian matrix for the Newton-Raphson
      *                        method.
-     * @param local_coupled_solutions Nodal values of solutions of the coupled
-     *                                processes of an element.
      */
     void assembleWithJacobianForPressureEquations(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        LocalCoupledSolutions const& local_coupled_solutions);
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, std::vector<double>& local_M_data,
+        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
+        std::vector<double>& local_Jac_data);
 
     unsigned getNumberOfIntegrationPoints() const override;
 
diff --git a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM-impl.h b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM-impl.h
index 32584dc29b4bd08ace4fc7f2ed3cf1e3abb069b3..9fc59d04d8bb4294502ce849847395e82f18995a 100644
--- a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM-impl.h
+++ b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM-impl.h
@@ -24,32 +24,33 @@ template <typename ShapeFunction, typename IntegrationMethod,
 void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
                                               DisplacementDim>::
     assembleWithJacobianForStaggeredScheme(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx, int const process_id,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, int const process_id,
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        LocalCoupledSolutions const& local_coupled_solutions)
+        LocalCoupledSolutions const& /*local_coupled_solutions*/)
 {
     if (process_id == _phase_field_process_id)
     {
         assembleWithJacobianForPhaseFieldEquations(
-            t, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
-            local_b_data, local_Jac_data, local_coupled_solutions);
+            t, local_x, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
+            local_b_data, local_Jac_data);
         return;
     }
 
     if (process_id == _heat_conduction_process_id)
     {
         assembleWithJacobianForHeatConductionEquations(
-            t, dt, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
-            local_b_data, local_Jac_data, local_coupled_solutions);
+            t, dt, local_x, local_xdot, dxdot_dx, dx_dx, local_M_data,
+            local_K_data, local_b_data, local_Jac_data);
         return;
     }
 
     // For the equations with deformation
     assembleWithJacobianForDeformationEquations(
-        t, dt, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
-        local_b_data, local_Jac_data, local_coupled_solutions);
+        t, dt, local_x, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
+        local_b_data, local_Jac_data);
 }
 
 template <typename ShapeFunction, typename IntegrationMethod,
@@ -57,25 +58,20 @@ template <typename ShapeFunction, typename IntegrationMethod,
 void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
                                               DisplacementDim>::
     assembleWithJacobianForDeformationEquations(
-        double const t, double const dt,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
         std::vector<double> const& /*local_xdot*/, const double /*dxdot_dx*/,
         const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
         std::vector<double>& /*local_K_data*/,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        LocalCoupledSolutions const& local_coupled_solutions)
+        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
 {
-    assert(local_coupled_solutions.local_coupled_xs.size() ==
+    assert(local_x.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 const d = local_x.template segment<phasefield_size>(phasefield_index);
+    auto const u =
+        local_x.template segment<displacement_size>(displacement_index);
+    auto const T =
+        local_x.template segment<temperature_size>(temperature_index);
 
     auto local_Jac = MathLib::createZeroedMatrix<
         typename ShapeMatricesType::template MatrixType<displacement_size,
@@ -158,22 +154,18 @@ template <typename ShapeFunction, typename IntegrationMethod,
 void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
                                               DisplacementDim>::
     assembleWithJacobianForHeatConductionEquations(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double /*dxdot_dx*/, const double /*dx_dx*/,
-        std::vector<double>& /*local_M_data*/,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double /*dxdot_dx*/,
+        const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
         std::vector<double>& /*local_K_data*/,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        LocalCoupledSolutions const& local_coupled_solutions)
+        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
 {
-    assert(local_coupled_solutions.local_coupled_xs.size() ==
+    assert(local_x.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 T = Eigen::Map<TemperatureVector const>(
-        &local_coupled_solutions.local_coupled_xs[temperature_index],
-        temperature_size);
+    auto const d = local_x.template segment<phasefield_size>(phasefield_index);
+    auto const T =
+        local_x.template segment<temperature_size>(temperature_index);
 
     auto T_dot = Eigen::Map<typename ShapeMatricesType::template VectorType<
         temperature_size> const>(local_xdot.data(), temperature_size);
@@ -258,19 +250,16 @@ template <typename ShapeFunction, typename IntegrationMethod,
 void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
                                               DisplacementDim>::
     assembleWithJacobianForPhaseFieldEquations(
-        double const t, std::vector<double> const& /*local_xdot*/,
-        const double /*dxdot_dx*/, const double /*dx_dx*/,
-        std::vector<double>& /*local_M_data*/,
+        double const t, Eigen::VectorXd const& local_x,
+        std::vector<double> const& /*local_xdot*/, const double /*dxdot_dx*/,
+        const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
         std::vector<double>& /*local_K_data*/,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        LocalCoupledSolutions const& local_coupled_solutions)
+        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
 {
-    assert(local_coupled_solutions.local_coupled_xs.size() ==
+    assert(local_x.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 d = local_x.template segment<phasefield_size>(phasefield_index);
 
     auto local_Jac = MathLib::createZeroedMatrix<
         typename ShapeMatricesType::template MatrixType<phasefield_size,
diff --git a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
index 0308107da7fc2fa3914802de542193066125c8a1..63fb82c1f3e69e27b124a0cb820e889429f3d1db 100644
--- a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
+++ b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
@@ -234,8 +234,9 @@ public:
     }
 
     void assembleWithJacobianForStaggeredScheme(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx, int const process_id,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, int const process_id,
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions) override;
@@ -356,25 +357,25 @@ private:
     }
 
     void assembleWithJacobianForDeformationEquations(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        LocalCoupledSolutions const& local_coupled_solutions);
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, std::vector<double>& local_M_data,
+        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
+        std::vector<double>& local_Jac_data);
 
     void assembleWithJacobianForHeatConductionEquations(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        LocalCoupledSolutions const& local_coupled_solutions);
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, std::vector<double>& local_M_data,
+        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
+        std::vector<double>& local_Jac_data);
 
     void assembleWithJacobianForPhaseFieldEquations(
-        double const t, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
-        LocalCoupledSolutions const& local_coupled_solutions);
+        double const t, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, std::vector<double>& local_M_data,
+        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
+        std::vector<double>& local_Jac_data);
 
     ThermoMechanicalPhaseFieldProcessData<DisplacementDim>& _process_data;
 
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h
index 44499e06acf7b89900fc10a861fa767fe76f5090..cc716057f214c88934b6e7e3778567d357badf90 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h
@@ -305,8 +305,9 @@ template <typename ShapeFunction, typename IntegrationMethod,
 void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
                                    DisplacementDim>::
     assembleWithJacobianForStaggeredScheme(
-        const double t, double const dt, const std::vector<double>& local_xdot,
-        const double dxdot_dx, const double dx_dx, int const process_id,
+        const double t, double const dt, Eigen::VectorXd const& local_x,
+        const std::vector<double>& local_xdot, const double dxdot_dx,
+        const double dx_dx, int const process_id,
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         const LocalCoupledSolutions& local_coupled_solutions)
@@ -315,14 +316,15 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
     if (process_id == _process_data.heat_conduction_process_id)
     {
         assembleWithJacobianForHeatConductionEquations(
-            t, dt, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
-            local_b_data, local_Jac_data, local_coupled_solutions);
+            t, dt, local_x, local_xdot, dxdot_dx, dx_dx, local_M_data,
+            local_K_data, local_b_data, local_Jac_data,
+            local_coupled_solutions);
         return;
     }
 
     // For the equations with deformation
     assembleWithJacobianForDeformationEquations(
-        t, dt, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
+        t, dt, local_x, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
         local_b_data, local_Jac_data, local_coupled_solutions);
 }
 
@@ -331,7 +333,7 @@ template <typename ShapeFunction, typename IntegrationMethod,
 void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
                                    DisplacementDim>::
     assembleWithJacobianForDeformationEquations(
-        const double t, double const dt,
+        const double t, double const dt, Eigen::VectorXd const& local_x,
         const std::vector<double>& /*local_xdot*/, const double /*dxdot_dx*/,
         const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
         std::vector<double>& /*local_K_data*/,
@@ -339,10 +341,7 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
         const LocalCoupledSolutions& local_coupled_solutions)
 {
     auto const local_T =
-        Eigen::Map<typename ShapeMatricesType::template VectorType<
-            temperature_size> const>(
-            &local_coupled_solutions.local_coupled_xs[temperature_index],
-            temperature_size);
+        local_x.template segment<temperature_size>(temperature_index);
 
     auto const local_T0 =
         Eigen::Map<typename ShapeMatricesType::template VectorType<
@@ -350,10 +349,8 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
             &local_coupled_solutions.local_coupled_xs0[temperature_index],
             temperature_size);
 
-    auto const u = Eigen::Map<typename ShapeMatricesType::template VectorType<
-        displacement_size> const>(
-        &local_coupled_solutions.local_coupled_xs[displacement_index],
-        displacement_size);
+    auto const u =
+        local_x.template segment<displacement_size>(displacement_index);
 
     auto local_Jac = MathLib::createZeroedMatrix<
         typename ShapeMatricesType::template MatrixType<displacement_size,
@@ -465,7 +462,7 @@ template <typename ShapeFunction, typename IntegrationMethod,
 void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
                                    DisplacementDim>::
     assembleWithJacobianForHeatConductionEquations(
-        const double t, double const dt,
+        const double t, double const dt, Eigen::VectorXd const& local_x,
         const std::vector<double>& /*local_xdot*/, const double /*dxdot_dx*/,
         const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
         std::vector<double>& /*local_K_data*/,
@@ -473,10 +470,7 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
         const LocalCoupledSolutions& local_coupled_solutions)
 {
     auto const local_T =
-        Eigen::Map<typename ShapeMatricesType::template VectorType<
-            temperature_size> const>(
-            &local_coupled_solutions.local_coupled_xs[temperature_index],
-            temperature_size);
+        local_x.template segment<temperature_size>(temperature_index);
 
     auto const local_T0 =
         Eigen::Map<typename ShapeMatricesType::template VectorType<
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
index 2197ea7c1b6972a606840390cb34815b82b172b7..830b78d71bc40bb6ac0f8f13371f7e8f9fe71892 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
@@ -134,8 +134,9 @@ public:
     }
 
     void assembleWithJacobianForStaggeredScheme(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx, int const process_id,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, int const process_id,
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions) override;
@@ -229,10 +230,11 @@ private:
      */
 
     void assembleWithJacobianForDeformationEquations(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, std::vector<double>& local_M_data,
+        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
+        std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions);
 
     /**
@@ -262,10 +264,11 @@ private:
      *                                displacement of an element.
      */
     void assembleWithJacobianForHeatConductionEquations(
-        double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
-        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
+        double const t, double const dt, Eigen::VectorXd const& local_x,
+        std::vector<double> const& local_xdot, const double dxdot_dx,
+        const double dx_dx, std::vector<double>& local_M_data,
+        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
+        std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions);
 
     std::size_t setSigma(double const* values);
diff --git a/ProcessLib/VectorMatrixAssembler.cpp b/ProcessLib/VectorMatrixAssembler.cpp
index b5db72ee845073e08115d8f69458ee401805d810..eed34a643709933c789d3b68655fb50d310c053e 100644
--- a/ProcessLib/VectorMatrixAssembler.cpp
+++ b/ProcessLib/VectorMatrixAssembler.cpp
@@ -70,15 +70,18 @@ void VectorMatrixAssembler::assemble(
     {
         auto local_coupled_xs0 = getCoupledLocalSolutions(cpl_xs->coupled_xs_t0,
                                                           indices_of_processes);
+
         auto local_coupled_xs =
             getCoupledLocalSolutions(x, indices_of_processes);
 
+        auto const local_x = MathLib::toVector(local_coupled_xs);
+
         ProcessLib::LocalCoupledSolutions local_coupled_solutions(
-            std::move(local_coupled_xs0), std::move(local_coupled_xs));
+            std::move(local_coupled_xs0));
 
         local_assembler.assembleForStaggeredScheme(
-            t, dt, process_id, _local_M_data, _local_K_data, _local_b_data,
-            local_coupled_solutions);
+            t, dt, local_x, process_id, _local_M_data, _local_K_data,
+            _local_b_data, local_coupled_solutions);
     }
 
     auto const num_r_c = indices.size();
@@ -138,16 +141,19 @@ void VectorMatrixAssembler::assembleWithJacobian(
     {
         auto local_coupled_xs0 = getCoupledLocalSolutions(cpl_xs->coupled_xs_t0,
                                                           indices_of_processes);
+
         auto local_coupled_xs =
             getCoupledLocalSolutions(x, indices_of_processes);
 
+        auto const local_x = MathLib::toVector(local_coupled_xs);
+
         ProcessLib::LocalCoupledSolutions local_coupled_solutions(
-            std::move(local_coupled_xs0), std::move(local_coupled_xs));
+            std::move(local_coupled_xs0));
 
         _jacobian_assembler->assembleWithJacobianForStaggeredScheme(
-            local_assembler, t, dt, local_xdot, dxdot_dx, dx_dx, process_id,
-            _local_M_data, _local_K_data, _local_b_data, _local_Jac_data,
-            local_coupled_solutions);
+            local_assembler, t, dt, local_x, local_xdot, dxdot_dx, dx_dx,
+            process_id, _local_M_data, _local_K_data, _local_b_data,
+            _local_Jac_data, local_coupled_solutions);
     }
 
     auto const num_r_c = indices.size();