From 42c19d6feb2e77b3c8bde13e698471159f23d6dc Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Wed, 21 Dec 2016 15:46:25 +0100
Subject: [PATCH] [Coupling] Added a coupling argument to global assembly
 functions of all process classes

coupled terms.
---
 .../GroundwaterFlowProcess.cpp                |  9 +-
 .../GroundwaterFlow/GroundwaterFlowProcess.h  |  7 +-
 ProcessLib/HT/HTProcess.cpp                   |  9 +-
 ProcessLib/HT/HTProcess.h                     |  7 +-
 .../HeatConduction/HeatConductionProcess.cpp  |  9 +-
 .../HeatConduction/HeatConductionProcess.h    |  7 +-
 ProcessLib/HydroMechanics/HydroMechanicsFEM.h |  4 +-
 .../HydroMechanics/HydroMechanicsProcess.h    |  9 +-
 .../HydroMechanics/HydroMechanicsProcess.h    |  4 +-
 .../SmallDeformationProcess.h                 |  9 +-
 ProcessLib/LiquidFlow/LiquidFlowProcess.cpp   |  9 +-
 ProcessLib/LiquidFlow/LiquidFlowProcess.h     | 12 +--
 ProcessLib/Process.cpp                        |  6 +-
 ProcessLib/Process.h                          |  7 +-
 .../RichardsFlow/RichardsFlowProcess.cpp      |  9 +-
 ProcessLib/RichardsFlow/RichardsFlowProcess.h |  7 +-
 .../SmallDeformationProcess.h                 |  9 +-
 ProcessLib/StaggeredCouplingTerm.cpp          |  7 ++
 ProcessLib/StaggeredCouplingTerm.h            | 15 +++
 ProcessLib/TES/TESProcess.cpp                 |  9 +-
 ProcessLib/TES/TESProcess.h                   |  7 +-
 .../TwoPhaseFlowWithPPProcess.cpp             |  9 +-
 .../TwoPhaseFlowWithPPProcess.h               |  7 +-
 .../TwoPhaseFlowWithPrhoProcess.cpp           | 11 ++-
 .../TwoPhaseFlowWithPrhoProcess.h             |  6 +-
 ProcessLib/VectorMatrixAssembler.cpp          | 93 ++++++++++++++++---
 ProcessLib/VectorMatrixAssembler.h            |  7 +-
 27 files changed, 227 insertions(+), 77 deletions(-)

diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
index 42367c1b05e..93589087db5 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
@@ -76,7 +76,9 @@ void GroundwaterFlowProcess::assembleConcreteProcess(const double t,
                                                      GlobalVector const& x,
                                                      GlobalMatrix& M,
                                                      GlobalMatrix& K,
-                                                     GlobalVector& b)
+                                                     GlobalVector& b,
+                                                     StaggeredCouplingTerm
+                                                     const& /*coupled_term*/)
 {
     DBUG("Assemble GroundwaterFlowProcess.");
 
@@ -89,7 +91,8 @@ void GroundwaterFlowProcess::assembleConcreteProcess(const double t,
 void GroundwaterFlowProcess::assembleWithJacobianConcreteProcess(
     const double t, GlobalVector const& x, GlobalVector const& xdot,
     const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b, GlobalMatrix& Jac)
+    GlobalVector& b, GlobalMatrix& Jac,
+    StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("AssembleWithJacobian GroundwaterFlowProcess.");
 
@@ -97,7 +100,7 @@ void GroundwaterFlowProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac);
+        dx_dx, M, K, b, Jac, coupled_term);
 }
 
 
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
index eaebcab80d6..69bb1a685a9 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
@@ -101,12 +101,15 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override;
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const& coupled_term
+                                ) override;
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     GroundwaterFlowProcessData _process_data;
 
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index 2ce02f4b28c..d8f4eddb63a 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -71,7 +71,9 @@ void HTProcess::assembleConcreteProcess(const double t,
                                         GlobalVector const& x,
                                         GlobalMatrix& M,
                                         GlobalMatrix& K,
-                                        GlobalVector& b)
+                                        GlobalVector& b,
+                                        StaggeredCouplingTerm const&
+                                        /*coupled_term*/)
 {
     DBUG("Assemble HTProcess.");
 
@@ -84,7 +86,8 @@ void HTProcess::assembleConcreteProcess(const double t,
 void HTProcess::assembleWithJacobianConcreteProcess(
     const double t, GlobalVector const& x, GlobalVector const& xdot,
     const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b, GlobalMatrix& Jac)
+    GlobalVector& b, GlobalMatrix& Jac,
+    StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("AssembleWithJacobian HTProcess.");
 
@@ -92,7 +95,7 @@ void HTProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac);
+        dx_dx, M, K, b, Jac, coupled_term);
 }
 
 }  // namespace HT
diff --git a/ProcessLib/HT/HTProcess.h b/ProcessLib/HT/HTProcess.h
index 63a9f058dd0..1da30247709 100644
--- a/ProcessLib/HT/HTProcess.h
+++ b/ProcessLib/HT/HTProcess.h
@@ -68,12 +68,15 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override;
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const& coupled_term
+                                ) override;
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     HTProcessData _process_data;
 
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.cpp b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
index b82273f8b6f..3ceed6b6453 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.cpp
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
@@ -72,7 +72,9 @@ void HeatConductionProcess::assembleConcreteProcess(const double t,
                                                     GlobalVector const& x,
                                                     GlobalMatrix& M,
                                                     GlobalMatrix& K,
-                                                    GlobalVector& b)
+                                                    GlobalVector& b,
+                                                    StaggeredCouplingTerm
+                                                    const& /*coupled_term*/)
 {
     DBUG("Assemble HeatConductionProcess.");
 
@@ -85,7 +87,8 @@ void HeatConductionProcess::assembleConcreteProcess(const double t,
 void HeatConductionProcess::assembleWithJacobianConcreteProcess(
     const double t, GlobalVector const& x, GlobalVector const& xdot,
     const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b, GlobalMatrix& Jac)
+    GlobalVector& b, GlobalMatrix& Jac,
+    StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("AssembleWithJacobian HeatConductionProcess.");
 
@@ -93,7 +96,7 @@ void HeatConductionProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac);
+        dx_dx, M, K, b, Jac, coupled_term);
 }
 
 void HeatConductionProcess::computeSecondaryVariableConcrete(const double t,
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.h b/ProcessLib/HeatConduction/HeatConductionProcess.h
index 4f59addd53d..5c94f4f95b6 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.h
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.h
@@ -52,12 +52,15 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override;
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const& coupled_term
+                                ) override;
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     HeatConductionProcessData _process_data;
 
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index 97e3f0d0eab..72540162dfb 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -200,7 +200,9 @@ public:
                               std::vector<double>& /*local_M_data*/,
                               std::vector<double>& /*local_K_data*/,
                               std::vector<double>& local_rhs_data,
-                              std::vector<double>& local_Jac_data) override
+                              std::vector<double>& local_Jac_data,
+                              LocalCouplingTerm const& /*coupled_term*/
+                             ) override
     {
         assert(local_x.size() == pressure_size + displacement_size);
 
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
index 8f0fd00cc45..0a6dc8c707a 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
@@ -109,7 +109,9 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const&
+                                 /*coupled_term*/) override
     {
         DBUG("Assemble HydroMechanicsProcess.");
 
@@ -122,7 +124,8 @@ private:
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override
     {
         DBUG("AssembleJacobian HydroMechanicsProcess.");
 
@@ -130,7 +133,7 @@ private:
         GlobalExecutor::executeMemberDereferenced(
             _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
             _local_assemblers, *_local_to_global_index_map, t, x, xdot,
-            dxdot_dx, dx_dx, M, K, b, Jac);
+            dxdot_dx, dx_dx, M, K, b, Jac, coupled_term);
     }
 
     void preTimestepConcreteProcess(GlobalVector const& x, double const t,
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
index 93dd3463446..18f260b49bb 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
@@ -70,7 +70,9 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const&
+                                 /*coupled_term*/) override
     {
         DBUG("Assemble HydroMechanicsProcess.");
 
diff --git a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
index 03adda654a9..efd4da63603 100644
--- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
+++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
@@ -68,7 +68,9 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const&
+                                 /*coupled_term*/) override
     {
         DBUG("Assemble SmallDeformationProcess.");
 
@@ -81,7 +83,8 @@ private:
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override
     {
         DBUG("AssembleWithJacobian SmallDeformationProcess.");
 
@@ -89,7 +92,7 @@ private:
         GlobalExecutor::executeMemberDereferenced(
             _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
             _local_assemblers, *_local_to_global_index_map, t, x, xdot,
-            dxdot_dx, dx_dx, M, K, b, Jac);
+            dxdot_dx, dx_dx, M, K, b, Jac, coupled_term);
     }
 
     void preTimestepConcreteProcess(GlobalVector const& x, double const t,
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index 1d9d9389dbc..b5f351e009a 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -92,7 +92,9 @@ void LiquidFlowProcess::assembleConcreteProcess(const double t,
                                                 GlobalVector const& x,
                                                 GlobalMatrix& M,
                                                 GlobalMatrix& K,
-                                                GlobalVector& b)
+                                                GlobalVector& b,
+                                                StaggeredCouplingTerm const&
+                                                /*coupled_term*/)
 {
     DBUG("Assemble LiquidFlowProcess.");
     // Call global assembler for each local assembly item.
@@ -104,7 +106,8 @@ void LiquidFlowProcess::assembleConcreteProcess(const double t,
 void LiquidFlowProcess::assembleWithJacobianConcreteProcess(
     const double t, GlobalVector const& x, GlobalVector const& xdot,
     const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b, GlobalMatrix& Jac)
+    GlobalVector& b, GlobalMatrix& Jac,
+    StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("AssembleWithJacobian LiquidFlowProcess.");
 
@@ -112,7 +115,7 @@ void LiquidFlowProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac);
+        dx_dx, M, K, b, Jac, coupled_term);
 }
 
 void LiquidFlowProcess::computeSecondaryVariableConcrete(const double t,
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.h b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
index ec69f532a3c..d7431f62589 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
@@ -80,11 +80,6 @@ public:
     ProcessType getProcessType() const override
                      {return ProcessLib::ProcessType::LiquidFlowProcess;}
 
-    MaterialLib::Fluid::FluidProperties* getFluidProperties() const
-    {
-        return _material_properties->getFluidProperties();
-    }
-
 private:
     void initializeConcreteProcess(
         NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -92,12 +87,15 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override;
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const& coupled_term
+                                ) override;
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     const int _gravitational_axis_id;
     const double _gravitational_acceleration;
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 4e69b98a00c..f41e79a61af 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -122,7 +122,8 @@ void Process::assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
                        GlobalMatrix& K, GlobalVector& b,
                        StaggeredCouplingTerm const& coupled_term)
 {
-    assembleConcreteProcess(t, x, M, K, b);
+    assembleConcreteProcess(t, x, M, K, b, coupled_term);
+
     _boundary_conditions.applyNaturalBC(t, x, K, b);
 }
 
@@ -133,7 +134,8 @@ void Process::assembleWithJacobian(const double t, GlobalVector const& x,
                                    GlobalVector& b, GlobalMatrix& Jac,
                                    StaggeredCouplingTerm const& coupled_term)
 {
-    assembleWithJacobianConcreteProcess(t, x, xdot, dxdot_dx, dx_dx, M, K, b, Jac);
+    assembleWithJacobianConcreteProcess(t, x, xdot, dxdot_dx, dx_dx, M, K, b,
+                                        Jac, coupled_term);
 
     // TODO apply BCs to Jacobian.
     _boundary_conditions.applyNaturalBC(t, x, K, b);
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 44df54f18ca..1f55b9e9b64 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -139,13 +139,16 @@ private:
 
     virtual void assembleConcreteProcess(const double t, GlobalVector const& x,
                                          GlobalMatrix& M, GlobalMatrix& K,
-                                         GlobalVector& b) = 0;
+                                         GlobalVector& b,
+                                         StaggeredCouplingTerm const&
+                                         coupled_term) = 0;
 
     virtual void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x,
         GlobalVector const& xdot, const double dxdot_dx,
         const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-        GlobalVector& b, GlobalMatrix& Jac) = 0;
+        GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) = 0;
 
     virtual void preTimestepConcreteProcess(GlobalVector const& /*x*/,
                                             const double /*t*/,
diff --git a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
index 0838ccf7d6e..f51c2d734e8 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
+++ b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
@@ -55,7 +55,9 @@ void RichardsFlowProcess::assembleConcreteProcess(const double t,
                                                   GlobalVector const& x,
                                                   GlobalMatrix& M,
                                                   GlobalMatrix& K,
-                                                  GlobalVector& b)
+                                                  GlobalVector& b,
+                                                  StaggeredCouplingTerm const&
+                                                  /*coupled_term*/)
 {
     DBUG("Assemble RichardsFlowProcess.");
 
@@ -68,7 +70,8 @@ void RichardsFlowProcess::assembleConcreteProcess(const double t,
 void RichardsFlowProcess::assembleWithJacobianConcreteProcess(
     const double t, GlobalVector const& x, GlobalVector const& xdot,
     const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b, GlobalMatrix& Jac)
+    GlobalVector& b, GlobalMatrix& Jac,
+    StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("AssembleWithJacobian RichardsFlowProcess.");
 
@@ -76,7 +79,7 @@ void RichardsFlowProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac);
+        dx_dx, M, K, b, Jac, coupled_term);
 }
 
 }  // namespace RichardsFlow
diff --git a/ProcessLib/RichardsFlow/RichardsFlowProcess.h b/ProcessLib/RichardsFlow/RichardsFlowProcess.h
index a0ec5904731..bcab9778a86 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowProcess.h
+++ b/ProcessLib/RichardsFlow/RichardsFlowProcess.h
@@ -52,12 +52,15 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override;
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const& coupled_term
+                                ) override;
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     RichardsFlowProcessData _process_data;
 
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.h b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
index 875f892c8fb..e0e4bf81528 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
@@ -146,7 +146,9 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const&
+                                 /*coupled_term*/) override
     {
         DBUG("Assemble SmallDeformationProcess.");
 
@@ -159,7 +161,8 @@ private:
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override
     {
         DBUG("AssembleWithJacobian SmallDeformationProcess.");
 
@@ -167,7 +170,7 @@ private:
         GlobalExecutor::executeMemberDereferenced(
             _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
             _local_assemblers, *_local_to_global_index_map, t, x, xdot,
-            dxdot_dx, dx_dx, M, K, b, Jac);
+            dxdot_dx, dx_dx, M, K, b, Jac, coupled_term);
     }
 
     void preTimestepConcreteProcess(GlobalVector const& x, double const t,
diff --git a/ProcessLib/StaggeredCouplingTerm.cpp b/ProcessLib/StaggeredCouplingTerm.cpp
index b588710308b..99e207857d5 100644
--- a/ProcessLib/StaggeredCouplingTerm.cpp
+++ b/ProcessLib/StaggeredCouplingTerm.cpp
@@ -22,5 +22,12 @@ const StaggeredCouplingTerm createVoidStaggeredCouplingTerm()
     return StaggeredCouplingTerm(coupled_pcsss, coupled_xs);
 }
 
+const LocalCouplingTerm createVoidLocalCouplingTerm()
+{
+    std::map<ProcessType, Process const&> coupled_pcsss;
+    std::map<ProcessType, const std::vector<double>> local_coupled_xs;
+    return LocalCouplingTerm(coupled_pcsss, std::move(local_coupled_xs));
+}
+
 } // end of ProcessLib
 
diff --git a/ProcessLib/StaggeredCouplingTerm.h b/ProcessLib/StaggeredCouplingTerm.h
index c7a36166d8e..76d056f14cb 100644
--- a/ProcessLib/StaggeredCouplingTerm.h
+++ b/ProcessLib/StaggeredCouplingTerm.h
@@ -35,7 +35,22 @@ struct StaggeredCouplingTerm
     std::map<ProcessType, GlobalVector const*> const& coupled_xs;
 };
 
+struct LocalCouplingTerm
+{
+    LocalCouplingTerm(
+        std::map<ProcessType, Process const&> const& coupled_processes_,
+        std::map<ProcessType, const std::vector<double>>&& local_coupled_xs_)
+    : coupled_processes(coupled_processes_),
+      local_coupled_xs(std::move(local_coupled_xs_))
+    {
+    }
+
+    std::map<ProcessType, Process const&> const& coupled_processes;
+    std::map<ProcessType, const std::vector<double>> const local_coupled_xs;
+};
+
 const StaggeredCouplingTerm createVoidStaggeredCouplingTerm();
+const LocalCouplingTerm createVoidLocalCouplingTerm();
 
 } // end of ProcessLib
 #endif /* OGS_STAGGERED_COUPLING_TERM_H */
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
index bad86ba6555..c5f2ccc466e 100644
--- a/ProcessLib/TES/TESProcess.cpp
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -226,7 +226,9 @@ void TESProcess::assembleConcreteProcess(const double t,
                                          GlobalVector const& x,
                                          GlobalMatrix& M,
                                          GlobalMatrix& K,
-                                         GlobalVector& b)
+                                         GlobalVector& b,
+                                         StaggeredCouplingTerm
+                                         const& /*coupled_term*/)
 {
     DBUG("Assemble TESProcess.");
 
@@ -239,12 +241,13 @@ void TESProcess::assembleConcreteProcess(const double t,
 void TESProcess::assembleWithJacobianConcreteProcess(
     const double t, GlobalVector const& x, GlobalVector const& xdot,
     const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b, GlobalMatrix& Jac)
+    GlobalVector& b, GlobalMatrix& Jac,
+    StaggeredCouplingTerm const& coupled_term)
 {
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac);
+        dx_dx, M, K, b, Jac, coupled_term);
 }
 
 void TESProcess::preTimestepConcreteProcess(GlobalVector const& x,
diff --git a/ProcessLib/TES/TESProcess.h b/ProcessLib/TES/TESProcess.h
index bed0ee6f6ac..5a340fbb6c7 100644
--- a/ProcessLib/TES/TESProcess.h
+++ b/ProcessLib/TES/TESProcess.h
@@ -59,14 +59,17 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override;
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const& coupled_term
+                                ) override;
 
     void initializeSecondaryVariables();
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     GlobalVector const& computeVapourPartialPressure(
         GlobalVector const& x,
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
index 669fa71c3f6..61e4c559d55 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
@@ -72,7 +72,9 @@ void TwoPhaseFlowWithPPProcess::assembleConcreteProcess(const double t,
                                                         GlobalVector const& x,
                                                         GlobalMatrix& M,
                                                         GlobalMatrix& K,
-                                                        GlobalVector& b)
+                                                        GlobalVector& b,
+                                                        StaggeredCouplingTerm
+                                                        const& /*coupled_term*/)
 {
     DBUG("Assemble TwoPhaseFlowWithPPProcess.");
     // Call global assembler for each local assembly item.
@@ -84,7 +86,8 @@ void TwoPhaseFlowWithPPProcess::assembleConcreteProcess(const double t,
 void TwoPhaseFlowWithPPProcess::assembleWithJacobianConcreteProcess(
     const double t, GlobalVector const& x, GlobalVector const& xdot,
     const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b, GlobalMatrix& Jac)
+    GlobalVector& b, GlobalMatrix& Jac,
+    StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("AssembleWithJacobian TwoPhaseFlowWithPPProcess.");
 
@@ -92,7 +95,7 @@ void TwoPhaseFlowWithPPProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac);
+        dx_dx, M, K, b, Jac, coupled_term);
 }
 
 }  // end of namespace
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h
index b3030be3987..155deaef901 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h
@@ -63,12 +63,15 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override;
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const& coupled_term
+                                ) override;
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     TwoPhaseFlowWithPPProcessData _process_data;
 
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
index 2125e18045d..9af42f4c48f 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
@@ -72,19 +72,22 @@ void TwoPhaseFlowWithPrhoProcess::assembleConcreteProcess(const double t,
                                                           GlobalVector const& x,
                                                           GlobalMatrix& M,
                                                           GlobalMatrix& K,
-                                                          GlobalVector& b)
+                                                          GlobalVector& b,
+                                                          StaggeredCouplingTerm
+                                                          const& coupled_term)
 {
     DBUG("Assemble TwoPhaseFlowWithPrhoProcess.");
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b);
+        *_local_to_global_index_map, t, x, M, K, b, coupled_term);
 }
 
 void TwoPhaseFlowWithPrhoProcess::assembleWithJacobianConcreteProcess(
     const double t, GlobalVector const& x, GlobalVector const& xdot,
     const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b, GlobalMatrix& Jac)
+    GlobalVector& b, GlobalMatrix& Jac,
+    StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("AssembleWithJacobian TwoPhaseFlowWithPrhoProcess.");
 
@@ -92,7 +95,7 @@ void TwoPhaseFlowWithPrhoProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac);
+        dx_dx, M, K, b, Jac, coupled_term);
 }
 void TwoPhaseFlowWithPrhoProcess::preTimestepConcreteProcess(
     GlobalVector const& x, double const t, double const dt)
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h
index ccf2163602f..decfa1103e6 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h
@@ -56,12 +56,14 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override;
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const& coupled_term) override;
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     void preTimestepConcreteProcess(GlobalVector const& x, const double t,
                                     const double delta_t) override;
diff --git a/ProcessLib/VectorMatrixAssembler.cpp b/ProcessLib/VectorMatrixAssembler.cpp
index dda159cf6ab..00628c642e7 100644
--- a/ProcessLib/VectorMatrixAssembler.cpp
+++ b/ProcessLib/VectorMatrixAssembler.cpp
@@ -15,6 +15,8 @@
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "LocalAssemblerInterface.h"
 
+#include "Process.h"
+
 namespace ProcessLib
 {
 VectorMatrixAssembler::VectorMatrixAssembler(
@@ -26,7 +28,8 @@ VectorMatrixAssembler::VectorMatrixAssembler(
 void VectorMatrixAssembler::assemble(
     const std::size_t mesh_item_id, LocalAssemblerInterface& local_assembler,
     const NumLib::LocalToGlobalIndexMap& dof_table, const double t,
-    const GlobalVector& x, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
+    const GlobalVector& x, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+    const StaggeredCouplingTerm& coupled_term)
 {
     auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
     auto const local_x = x.get(indices);
@@ -35,21 +38,50 @@ void VectorMatrixAssembler::assemble(
     _local_K_data.clear();
     _local_b_data.clear();
 
-    local_assembler.assemble(t, local_x, _local_M_data, _local_K_data, _local_b_data);
+    if (coupled_term.coupled_processes.empty())
+    {
+        local_assembler.assemble(t, local_x, _local_M_data, _local_K_data,
+                                 _local_b_data,
+                                 ProcessLib::createVoidLocalCouplingTerm());
+    }
+    else
+    {
+        std::map<ProcessLib::ProcessType, const std::vector<double>>
+            local_coupled_xs;
+        auto it = coupled_term.coupled_xs.begin();
+        while (it != coupled_term.coupled_xs.end())
+        {
+            GlobalVector const* coupled_x = it->second;
+            auto const local_coupled_x = coupled_x->get(indices);
+            BaseLib::insertMapIfKeyUniqueElseError(local_coupled_xs, it->first,
+                                                   local_coupled_x,
+                                                   "local_coupled_x");
+            it++;
+        }
+
+        ProcessLib::LocalCouplingTerm local_coupling_term(
+            coupled_term.coupled_processes, std::move(local_coupled_xs));
+
+        local_assembler.assemble(t, local_x, _local_M_data, _local_K_data,
+                                 _local_b_data, local_coupling_term);
+    }
 
     auto const num_r_c = indices.size();
     auto const r_c_indices =
         NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices);
 
-    if (!_local_M_data.empty()) {
+    if (!_local_M_data.empty())
+    {
         auto const local_M = MathLib::toMatrix(_local_M_data, num_r_c, num_r_c);
         M.add(r_c_indices, local_M);
     }
-    if (!_local_K_data.empty()) {
+    if (!_local_K_data.empty())
+    {
         auto const local_K = MathLib::toMatrix(_local_K_data, num_r_c, num_r_c);
         K.add(r_c_indices, local_K);
     }
-    if (!_local_b_data.empty()) {
+    if (!_local_b_data.empty())
+    {
         assert(_local_b_data.size() == num_r_c);
         b.add(indices, _local_b_data);
     }
@@ -60,7 +92,7 @@ void VectorMatrixAssembler::assembleWithJacobian(
     NumLib::LocalToGlobalIndexMap const& dof_table, const double t,
     GlobalVector const& x, GlobalVector const& xdot, const double dxdot_dx,
     const double dx_dx, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
-    GlobalMatrix& Jac)
+    GlobalMatrix& Jac, const StaggeredCouplingTerm& coupled_term)
 {
     auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
     auto const local_x = x.get(indices);
@@ -71,31 +103,64 @@ void VectorMatrixAssembler::assembleWithJacobian(
     _local_b_data.clear();
     _local_Jac_data.clear();
 
-    _jacobian_assembler->assembleWithJacobian(
-        local_assembler, t, local_x, local_xdot, dxdot_dx, dx_dx, _local_M_data,
-        _local_K_data, _local_b_data, _local_Jac_data);
+    if (coupled_term.coupled_processes.empty())
+    {
+        _jacobian_assembler->assembleWithJacobian(
+            local_assembler, t, local_x, local_xdot, dxdot_dx, dx_dx,
+            _local_M_data, _local_K_data, _local_b_data, _local_Jac_data,
+            ProcessLib::createVoidLocalCouplingTerm());
+    }
+    else
+    {
+        std::map<ProcessLib::ProcessType, const std::vector<double>>
+            local_coupled_xs;
+        auto it = coupled_term.coupled_xs.begin();
+        while (it != coupled_term.coupled_xs.end())
+        {
+            GlobalVector const* coupled_x = it->second;
+            auto const local_coupled_x = coupled_x->get(indices);
+            BaseLib::insertMapIfKeyUniqueElseError(local_coupled_xs, it->first,
+                                                   local_coupled_x,
+                                                   "local_coupled_x");
+            it++;
+        }
+
+        ProcessLib::LocalCouplingTerm local_coupling_term(
+            coupled_term.coupled_processes, std::move(local_coupled_xs));
+
+        _jacobian_assembler->assembleWithJacobian(
+            local_assembler, t, local_x, local_xdot, dxdot_dx, dx_dx,
+            _local_M_data, _local_K_data, _local_b_data, _local_Jac_data,
+            local_coupling_term);
+    }
 
     auto const num_r_c = indices.size();
     auto const r_c_indices =
         NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices);
 
-    if (!_local_M_data.empty()) {
+    if (!_local_M_data.empty())
+    {
         auto const local_M = MathLib::toMatrix(_local_M_data, num_r_c, num_r_c);
         M.add(r_c_indices, local_M);
     }
-    if (!_local_K_data.empty()) {
+    if (!_local_K_data.empty())
+    {
         auto const local_K = MathLib::toMatrix(_local_K_data, num_r_c, num_r_c);
         K.add(r_c_indices, local_K);
     }
-    if (!_local_b_data.empty()) {
+    if (!_local_b_data.empty())
+    {
         assert(_local_b_data.size() == num_r_c);
         b.add(indices, _local_b_data);
     }
-    if (!_local_Jac_data.empty()) {
+    if (!_local_Jac_data.empty())
+    {
         auto const local_Jac =
             MathLib::toMatrix(_local_Jac_data, num_r_c, num_r_c);
         Jac.add(r_c_indices, local_Jac);
-    } else {
+    }
+    else
+    {
         OGS_FATAL(
             "No Jacobian has been assembled! This might be due to programming "
             "errors in the local assembler of the current process.");
diff --git a/ProcessLib/VectorMatrixAssembler.h b/ProcessLib/VectorMatrixAssembler.h
index a017a00ab00..5466d6d71c0 100644
--- a/ProcessLib/VectorMatrixAssembler.h
+++ b/ProcessLib/VectorMatrixAssembler.h
@@ -12,6 +12,7 @@
 #include <vector>
 #include "NumLib/NumericsConfig.h"
 #include "AbstractJacobianAssembler.h"
+#include "StaggeredCouplingTerm.h"
 
 namespace NumLib
 {
@@ -38,7 +39,8 @@ public:
                   LocalAssemblerInterface& local_assembler,
                   NumLib::LocalToGlobalIndexMap const& dof_table,
                   double const t, GlobalVector const& x, GlobalMatrix& M,
-                  GlobalMatrix& K, GlobalVector& b);
+                  GlobalMatrix& K, GlobalVector& b,
+                  const StaggeredCouplingTerm& coupled_term);
 
     //! Assembles \c M, \c K, \c b, and the Jacobian \c Jac of the residual.
     //! \note The Jacobian must be assembled.
@@ -49,7 +51,8 @@ public:
                               GlobalVector const& xdot, const double dxdot_dx,
                               const double dx_dx, GlobalMatrix& M,
                               GlobalMatrix& K, GlobalVector& b,
-                              GlobalMatrix& Jac);
+                              GlobalMatrix& Jac,
+                              const StaggeredCouplingTerm& coupled_term);
 
 private:
     // temporary data only stored here in order to avoid frequent memory
-- 
GitLab