[PL] made Process Jacobian-assembly-ready

ProcessLib/LocalAssemblerInterface.cpp

@@ -8,43 +8,21 @@
  */
 
 #include "LocalAssemblerInterface.h"
+#include <cassert>
 #include "NumLib/DOF/DOFTableUtil.h"
 
 namespace ProcessLib
 {
-void LocalAssemblerInterface::assemble(
-    const std::size_t mesh_item_id,
-    const NumLib::LocalToGlobalIndexMap& dof_table, const double t,
-    const GlobalVector& x, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
-{
-    auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
-    auto const local_x = x.get(indices);
-    auto const r_c_indices =
-        NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices);
-
-    assembleConcrete(t, local_x, r_c_indices, M, K, b);
-}
-
-void LocalAssemblerInterface::assembleJacobian(
-    const std::size_t mesh_item_id,
-    const NumLib::LocalToGlobalIndexMap& dof_table, const double t,
-    const GlobalVector& x, GlobalMatrix& Jac)
-{
-    auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
-    auto const local_x = x.get(indices);
-    auto const r_c_indices =
-        NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices);
-
-    assembleJacobianConcrete(t, local_x, r_c_indices, Jac);
-}
-
-void LocalAssemblerInterface::assembleJacobianConcrete(
-        double const /*t*/, std::vector<double> const& /*local_x*/,
-        NumLib::LocalToGlobalIndexMap::RowColumnIndices const& /*indices*/,
-        GlobalMatrix& /*Jac*/)
+void LocalAssemblerInterface::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*/)
 {
     OGS_FATAL(
-        "The assembleJacobian() function is not implemented in the local "
+        "The assembleWithJacobian() function is not implemented in the local "
         "assembler.");
 }
 

ProcessLib/LocalAssemblerInterface.h

@@ -10,9 +10,13 @@
 #ifndef PROCESSLIB_LOCALASSEMBLERINTERFACE_H
 #define PROCESSLIB_LOCALASSEMBLERINTERFACE_H
 
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
 #include "NumLib/NumericsConfig.h"
 
+namespace NumLib
+{
+class LocalToGlobalIndexMap;
+}  // NumLib
+
 namespace ProcessLib
 {
 
@@ -26,15 +30,19 @@ class LocalAssemblerInterface
 public:
     virtual ~LocalAssemblerInterface() = default;
 
-    void assemble(std::size_t const mesh_item_id,
-                  NumLib::LocalToGlobalIndexMap const& dof_table,
-                  double const t, GlobalVector const& x, GlobalMatrix& M,
-                  GlobalMatrix& K, GlobalVector& b);
+    virtual void 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) = 0;
 
-    void assembleJacobian(std::size_t const mesh_item_id,
-                          NumLib::LocalToGlobalIndexMap const& dof_table,
-                          double const t, GlobalVector const& x,
-                          GlobalMatrix& Jac);
+    virtual void 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);
 
     virtual void preTimestep(std::size_t const mesh_item_id,
                              NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -45,17 +53,7 @@ public:
                               NumLib::LocalToGlobalIndexMap const& dof_table,
                               GlobalVector const& x);
 
-protected:
-    virtual void assembleConcrete(
-            double const t, std::vector<double> const& local_x,
-            NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
-            GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) = 0;
-
-    virtual void assembleJacobianConcrete(
-        double const t, std::vector<double> const& local_x,
-        NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
-        GlobalMatrix& Jac);
-
+private:
     virtual void preTimestepConcrete(std::vector<double> const& /*local_x*/,
                                      double const /*t*/, double const /*dt*/)
     {

ProcessLib/Process.cpp

@@ -20,6 +20,7 @@ namespace ProcessLib
 {
 Process::Process(
     MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
     SecondaryVariableCollection&& secondary_variables,
@@ -27,6 +28,7 @@ Process::Process(
     : _mesh(mesh),
       _secondary_variables(std::move(secondary_variables)),
       _named_function_caller(std::move(named_function_caller)),
+      _global_assembler(std::move(jacobian_assembler)),
       _process_variables(std::move(process_variables)),
       _boundary_conditions(parameters)
 {
@@ -115,43 +117,16 @@ void Process::assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
     _boundary_conditions.applyNaturalBC(t, x, K, b);
 }
 
-void Process::assembleJacobian(const double t, GlobalVector const& x,
-                               GlobalVector const& xdot, const double dxdot_dx,
-                               GlobalMatrix const& M, const double dx_dx,
-                               GlobalMatrix const& K, GlobalMatrix& Jac)
+void Process::assembleWithJacobian(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)
 {
-    assembleJacobianConcreteProcess(t, x, xdot, dxdot_dx, M, dx_dx, K, Jac);
-
-    // TODO In this method one could check if the user wants to use an
-    //      analytical or a numerical Jacobian. Then the right
-    //      assembleJacobianConcreteProcess() method will be chosen.
-    //      Additionally in the default implementation of said method one
-    //      could provide a fallback to a numerical Jacobian. However, that
-    //      would be in a sense implicit behaviour and it might be better to
-    //      just abort, as is currently the case.
-    //      In order to implement the Jacobian assembly entirely, in addition
-    //      to the operator() in VectorMatrixAssembler there has to be a method
-    //      that dispatches the Jacobian assembly.
-    //      Similarly, then the NeumannBC specialization of VectorMatrixAssembler
-    //      probably can be merged into the main class s.t. one has only one
-    //      type of VectorMatrixAssembler (for each equation type) with the
-    //      three methods assemble(), assembleJacobian() and assembleNeumannBC().
-    //      That list can be extended, e.g. by methods for the assembly of
-    //      source terms.
-    //      UPDATE: Probably it is better to keep a separate NeumannBC version of the
-    //      VectoMatrixAssembler since that will work for all kinds of processes.
-}
+    assembleWithJacobianConcreteProcess(t, x, xdot, dxdot_dx, dx_dx, M, K, b, Jac);
 
-void Process::assembleJacobianConcreteProcess(
-    const double /*t*/, GlobalVector const& /*x*/, GlobalVector const& /*xdot*/,
-    const double /*dxdot_dx*/, GlobalMatrix const& /*M*/,
-    const double /*dx_dx*/, GlobalMatrix const& /*K*/, GlobalMatrix& /*Jac*/)
-{
-    OGS_FATAL(
-        "The concrete implementation of this Process did not override the"
-        " assembleJacobianConcreteProcess() method."
-        " Hence, no analytical Jacobian is provided for this process"
-        " and the Newton-Raphson method cannot be used to solve it.");
+    // TODO apply BCs to Jacobian.
+    _boundary_conditions.applyNaturalBC(t, x, K, b);
 }
 
 void Process::constructDofTable()

ProcessLib/Process.h

@@ -40,6 +40,8 @@ public:
 
     Process(
         MeshLib::Mesh& mesh,
+        std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
+            jacobian_assembler,
         std::vector<std::unique_ptr<ParameterBase>> const& parameters,
         std::vector<std::reference_wrapper<ProcessVariable>>&&
             process_variables,
@@ -68,14 +70,14 @@ public:
     void assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
                   GlobalMatrix& K, GlobalVector& b) override final;
 
-    void assembleJacobian(const double t, GlobalVector const& x,
-                          GlobalVector const& xdot, const double dxdot_dx,
-                          GlobalMatrix const& M, const double dx_dx,
-                          GlobalMatrix const& K,
-                          GlobalMatrix& Jac) override final;
+    void assembleWithJacobian(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 final;
+
     std::vector<NumLib::IndexValueVector<GlobalIndexType>> const*
     getKnownSolutions(
-
         double const t) const override final
     {
         return _boundary_conditions.getKnownSolutions(t);
@@ -121,6 +123,12 @@ private:
                                          GlobalMatrix& M, GlobalMatrix& K,
                                          GlobalVector& b) = 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;
+
     virtual void preTimestepConcreteProcess(GlobalVector const& /*x*/,
                                             const double /*t*/,
                                             const double /*delta_t*/)
@@ -169,6 +177,8 @@ protected:
         _cached_secondary_variables;
     SecondaryVariableContext _secondary_variable_context;
 
+    VectorMatrixAssembler _global_assembler;
+
 private:
     unsigned const _integration_order = 2;
     GlobalSparsityPattern _sparsity_pattern;