From 0d6378357732b9d9120b5d2ef75101727e9f30f1 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Wed, 21 Dec 2016 17:29:02 +0100
Subject: [PATCH] [Coupling] Added local assembler of the coupled terms

 for the staggered scheme.
---
 ProcessLib/AbstractJacobianAssembler.h        | 13 +++++++++
 .../GroundwaterFlowProcess.cpp                |  4 +--
 ProcessLib/HT/HTProcess.cpp                   |  4 +--
 .../HeatConduction/HeatConductionProcess.cpp  |  4 +--
 ProcessLib/HydroMechanics/HydroMechanicsFEM.h |  4 +--
 .../HydroMechanics/HydroMechanicsProcess.h    |  5 ++--
 .../HydroMechanics/HydroMechanicsProcess.h    | 10 ++++---
 .../SmallDeformationProcess.h                 |  5 ++--
 ProcessLib/LiquidFlow/LiquidFlowProcess.cpp   |  4 +--
 ProcessLib/LocalAssemblerInterface.cpp        | 27 +++++++++++++++++++
 ProcessLib/LocalAssemblerInterface.h          | 17 ++++++++++++
 .../RichardsFlow/RichardsFlowProcess.cpp      |  4 +--
 .../SmallDeformationProcess.h                 |  5 ++--
 ProcessLib/TES/TESProcess.cpp                 |  4 +--
 .../TwoPhaseFlowWithPPProcess.cpp             |  4 +--
 ProcessLib/VectorMatrixAssembler.cpp          | 13 +++++----
 16 files changed, 93 insertions(+), 34 deletions(-)

diff --git a/ProcessLib/AbstractJacobianAssembler.h b/ProcessLib/AbstractJacobianAssembler.h
index 7989d383b16..b4267fadb7a 100644
--- a/ProcessLib/AbstractJacobianAssembler.h
+++ b/ProcessLib/AbstractJacobianAssembler.h
@@ -14,6 +14,7 @@
 namespace ProcessLib
 {
 class LocalAssemblerInterface;
+struct LocalCouplingTerm;
 
 //! Base class for Jacobian assemblers.
 class AbstractJacobianAssembler
@@ -29,6 +30,18 @@ public:
         std::vector<double>& local_K_data, std::vector<double>& local_b_data,
         std::vector<double>& local_Jac_data) = 0;
 
+    //! Assembles the Jacobian, the matrices \f$M\f$ and \f$K\f$, and the vector
+    //! \f$b\f$ with coupling.
+    virtual void coupling_assembleWithJacobian(
+        LocalAssemblerInterface& /*local_assembler*/, double const /*t*/,
+        std::vector<double> 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*/,
+        LocalCouplingTerm const& /*coupled_term*/) {}
+
     virtual ~AbstractJacobianAssembler() = default;
 };
 
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
index 93589087db5..11370da33e0 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
@@ -78,14 +78,14 @@ void GroundwaterFlowProcess::assembleConcreteProcess(const double t,
                                                      GlobalMatrix& K,
                                                      GlobalVector& b,
                                                      StaggeredCouplingTerm
-                                                     const& /*coupled_term*/)
+                                                     const& coupled_term)
 {
     DBUG("Assemble GroundwaterFlowProcess.");
 
     // 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 GroundwaterFlowProcess::assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index d8f4eddb63a..a9d0b6d19e1 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -73,14 +73,14 @@ void HTProcess::assembleConcreteProcess(const double t,
                                         GlobalMatrix& K,
                                         GlobalVector& b,
                                         StaggeredCouplingTerm const&
-                                        /*coupled_term*/)
+                                        coupled_term)
 {
     DBUG("Assemble HTProcess.");
 
     // 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 HTProcess::assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.cpp b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
index 3ceed6b6453..d181c782f10 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.cpp
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
@@ -74,14 +74,14 @@ void HeatConductionProcess::assembleConcreteProcess(const double t,
                                                     GlobalMatrix& K,
                                                     GlobalVector& b,
                                                     StaggeredCouplingTerm
-                                                    const& /*coupled_term*/)
+                                                    const& coupled_term)
 {
     DBUG("Assemble HeatConductionProcess.");
 
     // 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 HeatConductionProcess::assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index 72540162dfb..97e3f0d0eab 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -200,9 +200,7 @@ 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,
-                              LocalCouplingTerm const& /*coupled_term*/
-                             ) override
+                              std::vector<double>& local_Jac_data) override
     {
         assert(local_x.size() == pressure_size + displacement_size);
 
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
index 0a6dc8c707a..88d59391b9a 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
@@ -111,14 +111,15 @@ private:
                                  GlobalMatrix& M, GlobalMatrix& K,
                                  GlobalVector& b,
                                  StaggeredCouplingTerm const&
-                                 /*coupled_term*/) override
+                                 coupled_term) override
     {
         DBUG("Assemble HydroMechanicsProcess.");
 
         // 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_assemblers, *_local_to_global_index_map, t, x, M, K, b,
+            coupled_term);
     }
 
     void assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
index 18f260b49bb..8638c0fbf6e 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
@@ -72,20 +72,22 @@ private:
                                  GlobalMatrix& M, GlobalMatrix& K,
                                  GlobalVector& b,
                                  StaggeredCouplingTerm const&
-                                 /*coupled_term*/) override
+                                 coupled_term) override
     {
         DBUG("Assemble HydroMechanicsProcess.");
 
         // 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_assemblers, *_local_to_global_index_map, t, x, M, K, b,
+            coupled_term);
     }
 
     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 HydroMechanicsProcess.");
 
@@ -93,7 +95,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/SmallDeformation/SmallDeformationProcess.h b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
index efd4da63603..52fe3451935 100644
--- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
+++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
@@ -70,14 +70,15 @@ private:
                                  GlobalMatrix& M, GlobalMatrix& K,
                                  GlobalVector& b,
                                  StaggeredCouplingTerm const&
-                                 /*coupled_term*/) override
+                                 coupled_term) override
     {
         DBUG("Assemble SmallDeformationProcess.");
 
         // 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_assemblers, *_local_to_global_index_map, t, x, M, K, b,
+            coupled_term);
     }
 
     void assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index b5f351e009a..e9598f66f35 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -94,13 +94,13 @@ void LiquidFlowProcess::assembleConcreteProcess(const double t,
                                                 GlobalMatrix& K,
                                                 GlobalVector& b,
                                                 StaggeredCouplingTerm const&
-                                                /*coupled_term*/)
+                                                coupled_term)
 {
     DBUG("Assemble LiquidFlowProcess.");
     // 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 LiquidFlowProcess::assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/LocalAssemblerInterface.cpp b/ProcessLib/LocalAssemblerInterface.cpp
index 508f73b6bfa..ab83fdaef16 100644
--- a/ProcessLib/LocalAssemblerInterface.cpp
+++ b/ProcessLib/LocalAssemblerInterface.cpp
@@ -13,6 +13,19 @@
 
 namespace ProcessLib
 {
+
+void LocalAssemblerInterface::coupling_assemble(
+    double const /*t*/, std::vector<double> const& /*local_x*/,
+    std::vector<double>& /*local_M_data*/,
+    std::vector<double>& /*local_K_data*/,
+    std::vector<double>& /*local_b_data*/,
+    LocalCouplingTerm const& /*coupled_term*/)
+{
+    OGS_FATAL(
+        "The coupling_assemble() function is not implemented in the local "
+        "assembler.");
+}
+
 void LocalAssemblerInterface::assembleWithJacobian(
     double const /*t*/, std::vector<double> const& /*local_x*/,
     std::vector<double> const& /*local_xdot*/, const double /*dxdot_dx*/,
@@ -26,6 +39,20 @@ void LocalAssemblerInterface::assembleWithJacobian(
         "assembler.");
 }
 
+void LocalAssemblerInterface::coupling_assembleWithJacobian(
+    double const /*t*/, std::vector<double> 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*/,
+    LocalCouplingTerm const& /*coupled_term*/)
+{
+    OGS_FATAL(
+        "The coupling_assembleWithJacobian() function is not implemented in"
+        " the local.");
+}
+
 void LocalAssemblerInterface::computeSecondaryVariable(
                               std::size_t const mesh_item_id,
                               NumLib::LocalToGlobalIndexMap const& dof_table,
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index 8176485f679..aeb7329343d 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -11,6 +11,7 @@
 
 #include "NumLib/NumericsConfig.h"
 #include "MathLib/Point3d.h"
+#include "StaggeredCouplingTerm.h"
 
 namespace NumLib
 {
@@ -35,6 +36,12 @@ public:
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data) = 0;
 
+    virtual void coupling_assemble(
+        double const t, std::vector<double> const& local_x,
+        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
+        std::vector<double>& local_b_data,
+        LocalCouplingTerm const& coupled_term);
+
     virtual void assembleWithJacobian(double const t,
                                       std::vector<double> const& local_x,
                                       std::vector<double> const& local_xdot,
@@ -44,6 +51,16 @@ public:
                                       std::vector<double>& local_b_data,
                                       std::vector<double>& local_Jac_data);
 
+    virtual void coupling_assembleWithJacobian(double const t,
+                                      std::vector<double> 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,
+                                      LocalCouplingTerm const& coupled_term);
+
     virtual void computeSecondaryVariable(std::size_t const mesh_item_id,
                               NumLib::LocalToGlobalIndexMap const& dof_table,
                               const double t, GlobalVector const& x);
diff --git a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
index f51c2d734e8..7e5858a6dc4 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
+++ b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
@@ -57,14 +57,14 @@ void RichardsFlowProcess::assembleConcreteProcess(const double t,
                                                   GlobalMatrix& K,
                                                   GlobalVector& b,
                                                   StaggeredCouplingTerm const&
-                                                  /*coupled_term*/)
+                                                  coupled_term)
 {
     DBUG("Assemble RichardsFlowProcess.");
 
     // 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 RichardsFlowProcess::assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.h b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
index e0e4bf81528..c94a7b001a4 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
@@ -148,14 +148,15 @@ private:
                                  GlobalMatrix& M, GlobalMatrix& K,
                                  GlobalVector& b,
                                  StaggeredCouplingTerm const&
-                                 /*coupled_term*/) override
+                                 coupled_term) override
     {
         DBUG("Assemble SmallDeformationProcess.");
 
         // 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_assemblers, *_local_to_global_index_map, t, x, M, K, b,
+            coupled_term);
     }
 
     void assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
index c5f2ccc466e..1b9c77a9b9b 100644
--- a/ProcessLib/TES/TESProcess.cpp
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -228,14 +228,14 @@ void TESProcess::assembleConcreteProcess(const double t,
                                          GlobalMatrix& K,
                                          GlobalVector& b,
                                          StaggeredCouplingTerm
-                                         const& /*coupled_term*/)
+                                         const& coupled_term)
 {
     DBUG("Assemble TESProcess.");
 
     // 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 TESProcess::assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
index 61e4c559d55..79b9619105f 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
@@ -74,13 +74,13 @@ void TwoPhaseFlowWithPPProcess::assembleConcreteProcess(const double t,
                                                         GlobalMatrix& K,
                                                         GlobalVector& b,
                                                         StaggeredCouplingTerm
-                                                        const& /*coupled_term*/)
+                                                        const& coupled_term)
 {
     DBUG("Assemble TwoPhaseFlowWithPPProcess.");
     // 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 TwoPhaseFlowWithPPProcess::assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/VectorMatrixAssembler.cpp b/ProcessLib/VectorMatrixAssembler.cpp
index 00628c642e7..313ac0e0301 100644
--- a/ProcessLib/VectorMatrixAssembler.cpp
+++ b/ProcessLib/VectorMatrixAssembler.cpp
@@ -41,8 +41,7 @@ void VectorMatrixAssembler::assemble(
     if (coupled_term.coupled_processes.empty())
     {
         local_assembler.assemble(t, local_x, _local_M_data, _local_K_data,
-                                 _local_b_data,
-                                 ProcessLib::createVoidLocalCouplingTerm());
+                                 _local_b_data);
     }
     else
     {
@@ -62,8 +61,9 @@ void VectorMatrixAssembler::assemble(
         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);
+        local_assembler.coupling_assemble(t, local_x, _local_M_data,
+                                          _local_K_data, _local_b_data,
+                                          local_coupling_term);
     }
 
     auto const num_r_c = indices.size();
@@ -107,8 +107,7 @@ void VectorMatrixAssembler::assembleWithJacobian(
     {
         _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());
+            _local_M_data, _local_K_data, _local_b_data, _local_Jac_data);
     }
     else
     {
@@ -128,7 +127,7 @@ void VectorMatrixAssembler::assembleWithJacobian(
         ProcessLib::LocalCouplingTerm local_coupling_term(
             coupled_term.coupled_processes, std::move(local_coupled_xs));
 
-        _jacobian_assembler->assembleWithJacobian(
+        _jacobian_assembler->coupling_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);
-- 
GitLab