diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.h b/ProcessLib/BoundaryCondition/BoundaryCondition.h
index 460a0bd06146143fdf8f421d2a61940f0397de9e..67468d969420f9d9082a5097ea9376c8d4e00eea 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.h
@@ -37,9 +37,10 @@ class BoundaryCondition
 public:
     //! Applies natural BCs (i.e. non-Dirichlet BCs) to the stiffness matrix
     //! \c K and the vector \c b.
-    virtual void applyNaturalBC(const double /*t*/, GlobalVector const& /*x*/,
-                                GlobalMatrix& /*K*/, GlobalVector& /*b*/,
-                                GlobalMatrix* /*Jac*/)
+    virtual void applyNaturalBC(const double /*t*/,
+                                std::vector<GlobalVector*> const& /*x*/,
+                                int const /*process_id*/, GlobalMatrix& /*K*/,
+                                GlobalVector& /*b*/, GlobalMatrix* /*Jac*/)
     {
         // By default it is assumed that the BC is not a natural BC. Therefore
         // there is nothing to do here.
diff --git a/ProcessLib/BoundaryCondition/BoundaryConditionCollection.cpp b/ProcessLib/BoundaryCondition/BoundaryConditionCollection.cpp
index 7f88f1f8a5d08434c0229c06a528660d6bc7a47e..103cd70a805b317589af53fc392dd43e4e5f288c 100644
--- a/ProcessLib/BoundaryCondition/BoundaryConditionCollection.cpp
+++ b/ProcessLib/BoundaryCondition/BoundaryConditionCollection.cpp
@@ -12,15 +12,13 @@
 
 namespace ProcessLib
 {
-void BoundaryConditionCollection::applyNaturalBC(const double t,
-                                                 GlobalVector const& x,
-                                                 GlobalMatrix& K,
-                                                 GlobalVector& b,
-                                                 GlobalMatrix* Jac)
+void BoundaryConditionCollection::applyNaturalBC(
+    const double t, std::vector<GlobalVector*> const& x, int const process_id,
+    GlobalMatrix& K, GlobalVector& b, GlobalMatrix* Jac)
 {
     for (auto const& bc : _boundary_conditions)
     {
-        bc->applyNaturalBC(t, x, K, b, Jac);
+        bc->applyNaturalBC(t, x, process_id, K, b, Jac);
     }
 }
 
diff --git a/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h b/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h
index f9844d7558c653373e623c521b0aebd9cad5f41d..632be72894e3778b00dedb248a8a81bf6d1038c2 100644
--- a/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h
+++ b/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h
@@ -27,8 +27,9 @@ public:
     {
     }
 
-    void applyNaturalBC(const double t, GlobalVector const& x, GlobalMatrix& K,
-                        GlobalVector& b, GlobalMatrix* Jac);
+    void applyNaturalBC(const double t, std::vector<GlobalVector*> const& x,
+                        int const process_id, GlobalMatrix& K, GlobalVector& b,
+                        GlobalMatrix* Jac);
 
     std::vector<NumLib::IndexValueVector<GlobalIndexType>> const*
     getKnownSolutions(double const t, GlobalVector const& x) const
diff --git a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
index 587170d86728deb18c199d781fed61bb6650fe23..981a415f94d30872970d52135aaf5f18dfb02ff5 100644
--- a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
@@ -28,7 +28,8 @@ ConstraintDirichletBoundaryCondition::ConstraintDirichletBoundaryCondition(
     unsigned const integration_order, MeshLib::Mesh const& bulk_mesh,
     double const constraint_threshold, bool const lower,
     std::function<Eigen::Vector3d(std::size_t const, MathLib::Point3d const&,
-                                  double const, GlobalVector const&)>
+                                  double const,
+                                  std::vector<GlobalVector*> const&)>
         getFlux)
     : _parameter(parameter),
       _variable_id(variable_id),
@@ -107,20 +108,20 @@ ConstraintDirichletBoundaryCondition::ConstraintDirichletBoundaryCondition(
 }
 
 void ConstraintDirichletBoundaryCondition::preTimestep(
-    double t, std::vector<GlobalVector*> const& x, int const process_id)
+    double const t, std::vector<GlobalVector*> const& x,
+    int const /*process_id*/)
 {
     DBUG(
         "ConstraintDirichletBoundaryCondition::preTimestep: computing flux "
         "constraints");
-    auto const& solution = *x[process_id];
     for (auto const* boundary_element : _bc_mesh.getElements())
     {
         _flux_values[boundary_element->getID()] =
             _local_assemblers[boundary_element->getID()]->integrate(
-                solution, t,
+                x, t,
                 [this](std::size_t const element_id,
                        MathLib::Point3d const& pnt, double const t,
-                       GlobalVector const& x) {
+                       std::vector<GlobalVector*> const& x) {
                     return _getFlux(element_id, pnt, t, x);
                 });
     }
@@ -323,7 +324,7 @@ createConstraintDirichletBoundaryCondition(
         lower,
         [&constraining_process](std::size_t const element_id,
                                 MathLib::Point3d const& pnt, double const t,
-                                GlobalVector const& x) {
+                                std::vector<GlobalVector*> const& x) {
             return constraining_process.getFlux(element_id, pnt, t, x);
         });
 }
diff --git a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.h
index 186264607d394acd34e8e4be91e6aa84d789a89c..e119ec9e9d107b9b203e50534d10888ff0b32fad 100644
--- a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.h
@@ -58,10 +58,10 @@ public:
         bool const lower,
         std::function<Eigen::Vector3d(std::size_t const,
                                       MathLib::Point3d const&, double const,
-                                      GlobalVector const&)>
+                                      std::vector<GlobalVector*> const&)>
             getFlux);
 
-    void preTimestep(double t, std::vector<GlobalVector*> const& x,
+    void preTimestep(double const t, std::vector<GlobalVector*> const& x,
                      int const process_id) override;
 
     void getEssentialBCValues(
@@ -113,9 +113,10 @@ private:
     MeshLib::Mesh const& _bulk_mesh;
 
     /// The function _getFlux calculates the flux through the boundary element.
-    std::function<Eigen::Vector3d(
-            std::size_t const, MathLib::Point3d const&, double const,
-            GlobalVector const&)> _getFlux;
+    std::function<Eigen::Vector3d(std::size_t const, MathLib::Point3d const&,
+                                  double const,
+                                  std::vector<GlobalVector*> const&)>
+        _getFlux;
 };
 
 /// The function parses the config tree and creates a
diff --git a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
index 125b9c8f53226c4a6117b7a73d635a85c4298bb6..ef97c86120cb4a5f09b225a1c5ce25301e3a4540 100644
--- a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
@@ -50,10 +50,10 @@ public:
         default;
 
     virtual double integrate(
-        GlobalVector const& x, double const t,
-        std::function<Eigen::Vector3d(std::size_t const,
-                                      MathLib::Point3d const&, double const,
-                                      GlobalVector const&)> const& getFlux) = 0;
+        std::vector<GlobalVector*> const& x, double const t,
+        std::function<Eigen::Vector3d(
+            std::size_t const, MathLib::Point3d const&, double const,
+            std::vector<GlobalVector*> const&)> const& getFlux) = 0;
 };
 
 template <typename ShapeFunction, typename IntegrationMethod,
@@ -124,10 +124,10 @@ public:
     /// @param getFlux The function of the constraining process used to
     /// calculate the flux.
     double integrate(
-        GlobalVector const& x, double const t,
+        std::vector<GlobalVector*> const& x, double const t,
         std::function<Eigen::Vector3d(
             std::size_t const, MathLib::Point3d const&, double const,
-            GlobalVector const&)> const& getFlux) override
+            std::vector<GlobalVector*> const&)> const& getFlux) override
     {
         auto const n_integration_points =
             _integration_method.getNumberOfPoints();
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
index bc431a6fcd5ee226cf0fe70225e161b1d4190eda..bc716873285ea35bca49501ba9018b63587c491a 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
@@ -80,17 +80,18 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
 template <typename BoundaryConditionData,
           template <typename, typename, unsigned>
           class LocalAssemblerImplementation>
-void GenericNaturalBoundaryCondition<
-    BoundaryConditionData,
-    LocalAssemblerImplementation>::applyNaturalBC(const double t,
-                                                  const GlobalVector& x,
-                                                  GlobalMatrix& K,
-                                                  GlobalVector& b,
-                                                  GlobalMatrix* Jac)
+void GenericNaturalBoundaryCondition<BoundaryConditionData,
+                                     LocalAssemblerImplementation>::
+    applyNaturalBC(const double t,
+                   std::vector<GlobalVector*> const& x,
+                   int const process_id,
+                   GlobalMatrix& K,
+                   GlobalVector& b,
+                   GlobalMatrix* Jac)
 {
     GlobalExecutor::executeMemberOnDereferenced(
         &GenericNaturalBoundaryConditionLocalAssemblerInterface::assemble,
-        _local_assemblers, *_dof_table_boundary, t, x, K, b, Jac);
+        _local_assemblers, *_dof_table_boundary, t, x, process_id, K, b, Jac);
 }
 
 }  // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
index 07504106cf969afae31a250424c69bf87fece109..9d35e1dac29531e398b5b06ed2b0a4e8f42d2f82 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
@@ -35,8 +35,9 @@ public:
 
     /// Calls local assemblers which calculate their contributions to the global
     /// matrix and the right-hand-side.
-    void applyNaturalBC(const double t, GlobalVector const& x, GlobalMatrix& K,
-                        GlobalVector& b, GlobalMatrix* Jac) override;
+    void applyNaturalBC(const double t, std::vector<GlobalVector*> const& x,
+                        int const process_id, GlobalMatrix& K, GlobalVector& b,
+                        GlobalMatrix* Jac) override;
 
 private:
     /// Data used in the assembly of the specific boundary condition.
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h
index 0971216b1da448980eb3cd9276f2db08248f0643..aed6474a82f02071d37e9883165245f89720e694 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h
@@ -25,8 +25,8 @@ public:
     virtual void assemble(
         std::size_t const id,
         NumLib::LocalToGlobalIndexMap const& dof_table_boundary, double const t,
-        const GlobalVector& x, GlobalMatrix& K, GlobalVector& b,
-        GlobalMatrix* Jac) = 0;
+        std::vector<GlobalVector*> const& x, int const process_id,
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix* Jac) = 0;
 };
 
 template <typename ShapeFunction, typename IntegrationMethod,
diff --git a/ProcessLib/BoundaryCondition/HCNonAdvectiveFreeComponentFlowBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/HCNonAdvectiveFreeComponentFlowBoundaryConditionLocalAssembler.h
index cd2c8bf4a4201ada21f56f0adc2573b6233cda32..46c1b2747e21d59f6fbd0d3aff77738d22a4628a 100644
--- a/ProcessLib/BoundaryCondition/HCNonAdvectiveFreeComponentFlowBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/HCNonAdvectiveFreeComponentFlowBoundaryConditionLocalAssembler.h
@@ -69,8 +69,9 @@ public:
 
     void assemble(std::size_t const mesh_item_id,
                   NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
-                  double const t, const GlobalVector& x, GlobalMatrix& /*K*/,
-                  GlobalVector& b, GlobalMatrix* /*Jac*/) override
+                  double const t, std::vector<GlobalVector*> const& x,
+                  int const process_id, GlobalMatrix& /*K*/, GlobalVector& b,
+                  GlobalMatrix* /*Jac*/) override
     {
         NodalVectorType _local_rhs = NodalVectorType::Zero(_local_matrix_size);
         // Get element nodes for the interpolation from nodes to
@@ -83,7 +84,7 @@ public:
 
         auto const indices =
             NumLib::getIndices(mesh_item_id, dof_table_boundary);
-        std::vector<double> const local_values = x.get(indices);
+        std::vector<double> const local_values = x[process_id]->get(indices);
         std::size_t const bulk_element_id =
             _data.bulk_element_ids[Base::_element.getID()];
         std::size_t const bulk_face_id =
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
index 87072689c5618465721d88371f56c815720af3c0..11939130fa9d3a520edce452c4de0c75adaacdd4 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
@@ -46,9 +46,9 @@ public:
 
     void assemble(std::size_t const id,
                   NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
-                  double const t, const GlobalVector& /*x*/,
-                  GlobalMatrix& /*K*/, GlobalVector& b,
-                  GlobalMatrix* /*Jac*/) override
+                  double const t, std::vector<GlobalVector*> const& /*x*/,
+                  int const /*process_id*/, GlobalMatrix& /*K*/,
+                  GlobalVector& b, GlobalMatrix* /*Jac*/) override
     {
         _local_rhs.setZero();
 
diff --git a/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h
index 85bb366d0ecb05070b9113172eefaf018914590f..4895ee49e6f64689d0e0487ca597127a37acfcc2 100644
--- a/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h
@@ -61,12 +61,10 @@ NormalTractionBoundaryCondition<LocalAssemblerImplementation>::
 
 template <template <typename, typename, unsigned>
           class LocalAssemblerImplementation>
-void NormalTractionBoundaryCondition<
-    LocalAssemblerImplementation>::applyNaturalBC(const double t,
-                                                  const GlobalVector& x,
-                                                  GlobalMatrix& K,
-                                                  GlobalVector& b,
-                                                  GlobalMatrix* Jac)
+void NormalTractionBoundaryCondition<LocalAssemblerImplementation>::
+    applyNaturalBC(const double t, std::vector<GlobalVector*> const& x,
+                   int const /*process_id*/, GlobalMatrix& K, GlobalVector& b,
+                   GlobalMatrix* Jac)
 {
     GlobalExecutor::executeMemberOnDereferenced(
         &NormalTractionBoundaryConditionLocalAssemblerInterface::assemble,
diff --git a/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition.h b/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition.h
index b13faa3dbffc2ea795c71a8c8d96970b93588791..103aff697f54562e289dc3a1dcc2a163954e0a0d 100644
--- a/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition.h
@@ -44,8 +44,9 @@ public:
 
     /// Calls local assemblers which calculate their contributions to the global
     /// matrix and the right-hand-side.
-    void applyNaturalBC(const double t, GlobalVector const& x, GlobalMatrix& K,
-                        GlobalVector& b, GlobalMatrix* Jac) override;
+    void applyNaturalBC(const double t, std::vector<GlobalVector*> const& x,
+                        int const process_id, GlobalMatrix& K, GlobalVector& b,
+                        GlobalMatrix* Jac) override;
 
 private:
     MeshLib::Mesh const& _bc_mesh;
diff --git a/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h
index b6365f249751d04dbdb0f6488e0edc0166379bfa..53d6897a062b59aa26ef7a6b52fbabe4f64c8930 100644
--- a/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h
@@ -43,8 +43,8 @@ public:
     virtual void assemble(
         std::size_t const id,
         NumLib::LocalToGlobalIndexMap const& dof_table_boundary, double const t,
-        const GlobalVector& /*x*/, GlobalMatrix& /*K*/, GlobalVector& b,
-        GlobalMatrix* /*Jac*/) = 0;
+        std::vector<GlobalVector*> const& /*x*/, GlobalMatrix& /*K*/,
+        GlobalVector& b, GlobalMatrix* /*Jac*/) = 0;
     virtual ~NormalTractionBoundaryConditionLocalAssemblerInterface() = default;
 };
 
@@ -115,7 +115,7 @@ public:
 
     void assemble(std::size_t const id,
                   NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
-                  double const t, const GlobalVector& /*x*/,
+                  double const t, std::vector<GlobalVector*> const& /*x*/,
                   GlobalMatrix& /*K*/, GlobalVector& local_rhs,
                   GlobalMatrix* /*Jac*/) override
     {
diff --git a/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.cpp
index 97691e447279137434c1cfaf983d850bf949c9dc..695c61a5cc19f3567a55e9adb163c71d06660af6 100644
--- a/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.cpp
@@ -179,10 +179,9 @@ void PythonBoundaryCondition::getEssentialBCValues(
     }
 }
 
-void PythonBoundaryCondition::applyNaturalBC(const double t,
-                                             const GlobalVector& x,
-                                             GlobalMatrix& K, GlobalVector& b,
-                                             GlobalMatrix* Jac)
+void PythonBoundaryCondition::applyNaturalBC(
+    const double t, std::vector<GlobalVector*> const& x, int const process_id,
+    GlobalMatrix& K, GlobalVector& b, GlobalMatrix* Jac)
 {
     FlushStdoutGuard guard(_flush_stdout);
 
@@ -190,7 +189,8 @@ void PythonBoundaryCondition::applyNaturalBC(const double t,
     {
         GlobalExecutor::executeMemberOnDereferenced(
             &GenericNaturalBoundaryConditionLocalAssemblerInterface::assemble,
-            _local_assemblers, *_dof_table_boundary, t, x, K, b, Jac);
+            _local_assemblers, *_dof_table_boundary, t, x, process_id, K, b,
+            Jac);
     }
     catch (MethodNotOverriddenInDerivedClassException const& /*e*/)
     {
diff --git a/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.h b/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.h
index ef847a684b59f690f465a0c3fd57dfabba90bd50..be380f8649bbc56a1a05549ef22b31910bb6d9a3 100644
--- a/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.h
@@ -54,8 +54,9 @@ public:
         const double t, const GlobalVector& x,
         NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
 
-    void applyNaturalBC(const double t, const GlobalVector& x, GlobalMatrix& K,
-                        GlobalVector& b, GlobalMatrix* Jac) override;
+    void applyNaturalBC(const double t, std::vector<GlobalVector*> const& x,
+                        int const process_id, GlobalMatrix& K, GlobalVector& b,
+                        GlobalMatrix* Jac) override;
 
 private:
     //! Auxiliary data.
diff --git a/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h
index 81b1edf9d215a967d20aeb222e9ad9b3bac6f397..e7fac3247b59cfe6ad295710912114d4cb077a9e 100644
--- a/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h
@@ -48,8 +48,9 @@ public:
 
     void assemble(std::size_t const boundary_element_id,
                   NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
-                  double const t, const GlobalVector& x, GlobalMatrix& /*K*/,
-                  GlobalVector& b, GlobalMatrix* Jac) override
+                  double const t, std::vector<GlobalVector*> const& x,
+                  int const process_id, GlobalMatrix& /*K*/, GlobalVector& b,
+                  GlobalMatrix* Jac) override
     {
         using ShapeMatricesType =
             ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
@@ -104,7 +105,7 @@ public:
                             bulk_node_id, var, comp);
                     }
                     primary_variables_mat(element_node_id, global_component) =
-                        x[dof_idx];
+                        (*x[process_id])[dof_idx];
                 }
             }
         }
diff --git a/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
index 3e9a0be6c7a4f90d6af6d05ec93a979ae0d3a8de..64d9190bc2c520bae87f7ad6c137685d8e360d5d 100644
--- a/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
@@ -47,8 +47,9 @@ public:
     // TODO also implement derivative for Jacobian in Newton scheme.
     void assemble(std::size_t const id,
                   NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
-                  double const t, const GlobalVector& /*x*/, GlobalMatrix& K,
-                  GlobalVector& b, GlobalMatrix* /*Jac*/) override
+                  double const t, std::vector<GlobalVector*> const& /*x*/,
+                  int const /*process_id*/, GlobalMatrix& K, GlobalVector& b,
+                  GlobalMatrix* /*Jac*/) override
     {
         _local_K.setZero();
         _local_rhs.setZero();
diff --git a/ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryConditionLocalAssembler.h
index 9aa29b5536c8a180a73cd2536f31e34e2388384b..afcd4f9b264bd168f4b70f4d133c6586f6616ecb 100644
--- a/ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryConditionLocalAssembler.h
@@ -57,8 +57,9 @@ public:
 
     void assemble(std::size_t const mesh_item_id,
                   NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
-                  double const t, const GlobalVector& x, GlobalMatrix& /*K*/,
-                  GlobalVector& b, GlobalMatrix* /*Jac*/) override
+                  double const t, std::vector<GlobalVector*> const& x,
+                  int const process_id, GlobalMatrix& /*K*/, GlobalVector& b,
+                  GlobalMatrix* /*Jac*/) override
     {
         NodalVectorType _local_rhs(_local_matrix_size);
         _local_rhs.setZero();
@@ -87,9 +88,9 @@ public:
         auto const indices_other_variable = NumLib::getIndices(
             mesh_item_id, _data.dof_table_boundary_other_variable);
         std::vector<double> const local_current_variable =
-            x.get(indices_current_variable);
+            x[process_id]->get(indices_current_variable);
         std::vector<double> const local_other_variable =
-            x.get(indices_other_variable);
+            x[process_id]->get(indices_other_variable);
 
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index 0611892de5a41005e3a0d29fc39ac97105cacd3b..9faec0fd6276265e23bb298db143502465ef5a03 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -148,29 +148,22 @@ void ComponentTransportProcess::assembleWithJacobianConcreteProcess(
         dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 }
 
-Eigen::Vector3d ComponentTransportProcess::getFlux(std::size_t const element_id,
-                                                   MathLib::Point3d const& p,
-                                                   double const t,
-                                                   GlobalVector const& x) const
+Eigen::Vector3d ComponentTransportProcess::getFlux(
+    std::size_t const element_id,
+    MathLib::Point3d const& p,
+    double const t,
+    std::vector<GlobalVector*> const& x) const
 {
     std::vector<GlobalIndexType> indices_cache;
     auto const r_c_indices = NumLib::getRowColumnIndices(
         element_id, *_local_to_global_index_map, indices_cache);
 
-    if (_use_monolithic_scheme)
-    {
-        std::vector<double> local_x(x.get(r_c_indices.rows));
-
-        return _local_assemblers[element_id]->getFlux(p, t, local_x);
-    }
-
-        std::vector<std::vector<GlobalIndexType>>
-            indices_of_all_coupled_processes{
-                _coupled_solutions->coupled_xs.size(), r_c_indices.rows};
-        auto const local_xs = getCoupledLocalSolutions(
-            _coupled_solutions->coupled_xs, indices_of_all_coupled_processes);
+    std::vector<std::vector<GlobalIndexType>> indices_of_all_coupled_processes{
+        x.size(), r_c_indices.rows};
+    auto const local_xs =
+        getCoupledLocalSolutions(x, indices_of_all_coupled_processes);
 
-        return _local_assemblers[element_id]->getFlux(p, t, local_xs);
+    return _local_assemblers[element_id]->getFlux(p, t, local_xs);
 }
 
 void ComponentTransportProcess::
@@ -237,8 +230,7 @@ void ComponentTransportProcess::postTimestepConcreteProcess(
 
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
-    _surfaceflux->integrate(*x[process_id], t, *this, process_id,
-                            _integration_order, _mesh,
+    _surfaceflux->integrate(x, t, *this, process_id, _integration_order, _mesh,
                             pv.getActiveElementIDs());
     _surfaceflux->save(t);
 }
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.h b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
index 2486a1a4dadd5c94cbbdda95df01be1e1979da7d..1e23711b29a096066fbb5b98c036d1b742f7daa0 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
@@ -115,7 +115,7 @@ public:
 
     Eigen::Vector3d getFlux(std::size_t const element_id,
                             MathLib::Point3d const& p, double const t,
-                            GlobalVector const& x) const override;
+                            std::vector<GlobalVector*> const& x) const override;
 
     std::vector<std::pair<int, std::string>> const&
     getProcessIDToComponentNameMap() const
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
index 7e2544aa93b53db4539f543318ff336a728556d2..15fe117171531cfd65ce010c63cd0f2dd39d8074 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
@@ -50,13 +50,13 @@ public:
     Eigen::Vector3d getFlux(std::size_t element_id,
                             MathLib::Point3d const& p,
                             double const t,
-                            GlobalVector const& x) const override
+                            std::vector<GlobalVector*> const& x) const override
     {
         // fetch local_x from primary variable
         std::vector<GlobalIndexType> indices_cache;
         auto const r_c_indices = NumLib::getRowColumnIndices(
             element_id, *_local_to_global_index_map, indices_cache);
-        std::vector<double> local_x(x.get(r_c_indices.rows));
+        std::vector<double> local_x(x[0]->get(r_c_indices.rows));
 
         return _local_assemblers[element_id]->getFlux(p, t, local_x);
     }
@@ -79,9 +79,8 @@ public:
 
         ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
-        _surfaceflux->integrate(*x[process_id], t, *this, process_id,
-                                _integration_order, _mesh,
-                                pv.getActiveElementIDs());
+        _surfaceflux->integrate(x, t, *this, process_id, _integration_order,
+                                _mesh, pv.getActiveElementIDs());
         _surfaceflux->save(t);
     }
 
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index c4637fa191f2e34a448ba81471519d7a554af25d..b61a7efa25bfafca38631f5ca355a746d7c0bf2b 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -243,13 +243,16 @@ void HTProcess::setCoupledSolutionsOfPreviousTimeStep()
 Eigen::Vector3d HTProcess::getFlux(std::size_t element_id,
                                    MathLib::Point3d const& p,
                                    double const t,
-                                   GlobalVector const& x) const
+                                   std::vector<GlobalVector*> const& x) const
 {
     // fetch local_x from primary variable
     std::vector<GlobalIndexType> indices_cache;
     auto const r_c_indices = NumLib::getRowColumnIndices(
         element_id, *_local_to_global_index_map, indices_cache);
-    std::vector<double> local_x(x.get(r_c_indices.rows));
+    std::vector<std::vector<GlobalIndexType>> indices_of_all_coupled_processes{
+        x.size(), r_c_indices.rows};
+    auto const local_x =
+        getCoupledLocalSolutions(x, indices_of_all_coupled_processes);
 
     return _local_assemblers[element_id]->getFlux(p, t, local_x);
 }
@@ -279,8 +282,7 @@ void HTProcess::postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
 
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
-    _surfaceflux->integrate(*x[process_id], t, *this, process_id,
-                            _integration_order, _mesh,
+    _surfaceflux->integrate(x, t, *this, process_id, _integration_order, _mesh,
                             pv.getActiveElementIDs());
     _surfaceflux->save(t);
 }
diff --git a/ProcessLib/HT/HTProcess.h b/ProcessLib/HT/HTProcess.h
index 17ef8de52d1c78aeabb3f419fee3c65338088d1d..1d6a7e313c0e432733d995dcb45d2ebe4f46a702 100644
--- a/ProcessLib/HT/HTProcess.h
+++ b/ProcessLib/HT/HTProcess.h
@@ -76,7 +76,7 @@ public:
     Eigen::Vector3d getFlux(std::size_t element_id,
                             MathLib::Point3d const& p,
                             double const t,
-                            GlobalVector const& x) const override;
+                            std::vector<GlobalVector*> const& x) const override;
 
     void setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
         int const process_id) override;
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 48d21c8d52ccc2534146bd6c4a96e9919c0c3d83..f5582c6a6628a3d3757b305da0171ee71c6b1a68 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -205,7 +205,7 @@ void Process::assemble(const double t, double const dt,
     assembleConcreteProcess(t, dt, x, process_id, M, K, b);
 
     // the last argument is for the jacobian, nullptr is for a unused jacobian
-    _boundary_conditions[process_id].applyNaturalBC(t, *x[process_id], K, b,
+    _boundary_conditions[process_id].applyNaturalBC(t, x, process_id, K, b,
                                                     nullptr);
 
     // the last argument is for the jacobian, nullptr is for a unused jacobian
@@ -228,7 +228,7 @@ void Process::assembleWithJacobian(const double t, double const dt,
                                         process_id, M, K, b, Jac);
 
     // TODO: apply BCs to Jacobian.
-    _boundary_conditions[process_id].applyNaturalBC(t, *x[process_id], K, b,
+    _boundary_conditions[process_id].applyNaturalBC(t, x, process_id, K, b,
                                                     &Jac);
 
     // the last argument is for the jacobian, nullptr is for a unused jacobian
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index b901f5bfc14db811c8b6f4eab31fd6c5fff59fdb..2c10d9c9ad50288b6740a924898f17254104ca81 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -146,10 +146,11 @@ public:
 
     // Used as a call back for CalculateSurfaceFlux process.
 
-    virtual Eigen::Vector3d getFlux(std::size_t /*element_id*/,
-                                    MathLib::Point3d const& /*p*/,
-                                    double const /*t*/,
-                                    GlobalVector const& /*x*/) const
+    virtual Eigen::Vector3d getFlux(
+        std::size_t /*element_id*/,
+        MathLib::Point3d const& /*p*/,
+        double const /*t*/,
+        std::vector<GlobalVector*> const& /*x*/) const
     {
         return Eigen::Vector3d{};
     }
diff --git a/ProcessLib/SurfaceFlux/SurfaceFlux.cpp b/ProcessLib/SurfaceFlux/SurfaceFlux.cpp
index 90a9a8f2f76c5f134ea0902047973628b053f84c..4f3911284077eaa525ae693c116101bee6ac505a 100644
--- a/ProcessLib/SurfaceFlux/SurfaceFlux.cpp
+++ b/ProcessLib/SurfaceFlux/SurfaceFlux.cpp
@@ -54,14 +54,14 @@ SurfaceFlux::SurfaceFlux(
 }
 
 void SurfaceFlux::integrate(
-    GlobalVector const& x,
+    std::vector<GlobalVector*> const& x,
     MeshLib::PropertyVector<double>& balance,
     double const t,
     MeshLib::Mesh const& bulk_mesh,
     std::vector<std::size_t> const& active_element_ids,
-    std::function<Eigen::Vector3d(std::size_t const, MathLib::Point3d const&,
-                                  double const, GlobalVector const&)> const&
-        getFlux)
+    std::function<Eigen::Vector3d(
+        std::size_t const, MathLib::Point3d const&, double const,
+        std::vector<GlobalVector*> const&)> const& getFlux)
 {
     DBUG("Integrate SurfaceFlux.");
 
diff --git a/ProcessLib/SurfaceFlux/SurfaceFlux.h b/ProcessLib/SurfaceFlux/SurfaceFlux.h
index 30a32cffe5abd1ebeb4c7ff0e34f2f36ec979aa2..7580ba35ea43505d4551887ea5c927d8c1587204 100644
--- a/ProcessLib/SurfaceFlux/SurfaceFlux.h
+++ b/ProcessLib/SurfaceFlux/SurfaceFlux.h
@@ -40,14 +40,14 @@ public:
     /// @param getFlux function that calculates the flux in the integration
     /// points of the face elements of the bulk element that belongs to the
     /// surface.
-    void integrate(GlobalVector const& x,
+    void integrate(std::vector<GlobalVector*> const& x,
                    MeshLib::PropertyVector<double>& balance,
                    double const t,
                    MeshLib::Mesh const& bulk_mesh,
                    std::vector<std::size_t> const& active_element_ids,
                    std::function<Eigen::Vector3d(
                        std::size_t const, MathLib::Point3d const&, double const,
-                       GlobalVector const&)> const& getFlux);
+                       std::vector<GlobalVector*> const&)> const& getFlux);
 
 private:
     // the local assemblers for each element of the mesh
diff --git a/ProcessLib/SurfaceFlux/SurfaceFluxData.h b/ProcessLib/SurfaceFlux/SurfaceFluxData.h
index aa79b3188478e81d902ea64619b5ad4de13c57c3..b03e379998a7a2f9d27e19c9207a49c2691ab13e 100644
--- a/ProcessLib/SurfaceFlux/SurfaceFluxData.h
+++ b/ProcessLib/SurfaceFlux/SurfaceFluxData.h
@@ -76,9 +76,9 @@ struct SurfaceFluxData
             std::move(surfaceflux_out_fname));
     }
 
-    void integrate(GlobalVector const& x, double const t, Process const& p,
-                   int const process_id, int const integration_order,
-                   MeshLib::Mesh const& bulk_mesh,
+    void integrate(std::vector<GlobalVector*> const& x, double const t,
+                   Process const& p, int const process_id,
+                   int const integration_order, MeshLib::Mesh const& bulk_mesh,
                    std::vector<std::size_t> const& active_element_ids)
     {
         auto* const surfaceflux_pv = MeshLib::getOrCreateMeshProperty<double>(
@@ -93,7 +93,7 @@ struct SurfaceFluxData
         surfaceflux_process.integrate(
             x, *surfaceflux_pv, t, bulk_mesh, active_element_ids,
             [&p](std::size_t const element_id, MathLib::Point3d const& pnt,
-                 double const t, GlobalVector const& x) {
+                 double const t, std::vector<GlobalVector*> const& x) {
                 return p.getFlux(element_id, pnt, t, x);
             });
     }
diff --git a/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h b/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h
index c4df9da0baad7c170ba34e1989ae8d0eb0f6d05b..8d499b2621fa9aa5066e97e0b976260c65c83cf2 100644
--- a/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h
+++ b/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h
@@ -29,14 +29,15 @@ class SurfaceFluxLocalAssemblerInterface
 public:
     virtual ~SurfaceFluxLocalAssemblerInterface() = default;
 
-    virtual void integrate(std::size_t const element_id,
-                           GlobalVector const& x,
-                           MeshLib::PropertyVector<double>& specific_flux,
-                           double const t,
-                           MeshLib::Mesh const& bulk_mesh,
-                           std::function<Eigen::Vector3d(
-                               std::size_t const, MathLib::Point3d const&,
-                               double const, GlobalVector const&)>) = 0;
+    virtual void integrate(
+        std::size_t const element_id,
+        std::vector<GlobalVector*> const& x,
+        MeshLib::PropertyVector<double>& specific_flux,
+        double const t,
+        MeshLib::Mesh const& bulk_mesh,
+        std::function<Eigen::Vector3d(std::size_t const,
+                                      MathLib::Point3d const&, double const,
+                                      std::vector<GlobalVector*> const&)>) = 0;
 };
 
 template <typename ShapeFunction, typename IntegrationMethod,
@@ -110,13 +111,13 @@ public:
     /// points of the face elements of the bulk element that belongs to the
     /// surface.
     void integrate(std::size_t const element_id,
-                   GlobalVector const& x,
+                   std::vector<GlobalVector*> const& x,
                    MeshLib::PropertyVector<double>& specific_flux,
                    double const t,
                    MeshLib::Mesh const& bulk_mesh,
                    std::function<Eigen::Vector3d(
                        std::size_t const, MathLib::Point3d const&, double const,
-                       GlobalVector const&)>
+                       std::vector<GlobalVector*> const&)>
                        getFlux) override
     {
         auto surface_element_normal =