Merge pull request #2005 from endJunction/PreAssemble

Add ode system preAssemble() call.

NumLib/ODESolver/ODESystem.h

@@ -44,6 +44,9 @@ public:
     static const ODESystemTag ODETag =
         ODESystemTag::FirstOrderImplicitQuasilinear;
 
+    //! Calls process' pre-assembly with the provided state (\c t, \c x).
+    virtual void preAssemble(const double t, GlobalVector const& x) = 0;
+
     //! Assemble \c M, \c K and \c b at the provided state (\c t, \c x).
     virtual void assemble(
         const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
@@ -73,6 +76,10 @@ class ODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
                        NonlinearSolverTag::Picard>
 {
 public:
+
+    //! Calls process' pre-assembly with the provided state (\c t, \c x).
+    virtual void preAssemble(const double t, GlobalVector const& x) = 0;
+
     /*! Assemble \c M, \c K, \c b and the Jacobian
      * \f$ \mathtt{Jac} := \partial r/\partial x_N \f$
      * at the provided state (\c t, \c x).

NumLib/ODESolver/TimeDiscretizedODESystem.cpp

@@ -84,6 +84,7 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
     _b->setZero();
     _Jac->setZero();
 
+    _ode.preAssemble(t, x_curr);
     _ode.assembleWithJacobian(t, x_curr, xdot, dxdot_dx, dx_dx, *_M, *_K, *_b,
                               *_Jac);
 
@@ -186,6 +187,7 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
     _K->setZero();
     _b->setZero();
 
+    _ode.preAssemble(t, x_curr);
     _ode.assemble(t, x_curr, *_M, *_K, *_b);
 
     LinAlg::finalizeAssembly(*_M);

ProcessLib/LocalAssemblerInterface.h

@@ -35,6 +35,9 @@ class LocalAssemblerInterface
 public:
     virtual ~LocalAssemblerInterface() = default;
 
+    virtual void preAssemble(double const /*t*/,
+                             std::vector<double> const& /*local_x*/){};
+
     virtual void assemble(double const t, std::vector<double> const& local_x,
                           std::vector<double>& local_M_data,
                           std::vector<double>& local_K_data,

ProcessLib/Process.cpp

@@ -132,6 +132,11 @@ MathLib::MatrixSpecifications Process::getMatrixSpecifications() const
             &l.getGhostIndices(), &_sparsity_pattern};
 }
 
+void Process::preAssemble(const double t, GlobalVector const& x)
+{
+    preAssembleConcreteProcess(t, x);
+}
+
 void Process::assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
                        GlobalMatrix& K, GlobalVector& b)
 {

ProcessLib/Process.h

@@ -77,6 +77,8 @@ public:
         _coupled_solutions = coupled_solutions;
     }
 
+    void preAssemble(const double t, GlobalVector const& x) override final;
+
     virtual void setCoupledSolutionsForStaggeredSchemeToLocalAssemblers() {}
     void assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
                   GlobalMatrix& K, GlobalVector& b) final;
@@ -142,6 +144,11 @@ private:
         MeshLib::Mesh const& mesh,
         unsigned const integration_order) = 0;
 
+    virtual void preAssembleConcreteProcess(const double /*t*/,
+                                            GlobalVector const& /*x*/)
+    {
+    }
+
     virtual void assembleConcreteProcess(const double t, GlobalVector const& x,
                                          GlobalMatrix& M, GlobalMatrix& K,
                                          GlobalVector& b) = 0;

ProcessLib/VectorMatrixAssembler.cpp

@@ -55,6 +55,17 @@ VectorMatrixAssembler::VectorMatrixAssembler(
 {
 }
 
+void VectorMatrixAssembler::preAssemble(
+    const std::size_t mesh_item_id, LocalAssemblerInterface& local_assembler,
+    const NumLib::LocalToGlobalIndexMap& dof_table, const double t,
+    const GlobalVector& x)
+{
+    auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
+    auto const local_x = x.get(indices);
+
+    local_assembler.preAssemble(t, local_x);
+}
+
 void VectorMatrixAssembler::assemble(
     const std::size_t mesh_item_id, LocalAssemblerInterface& local_assembler,
     const NumLib::LocalToGlobalIndexMap& dof_table, const double t,

ProcessLib/VectorMatrixAssembler.h

@@ -33,6 +33,11 @@ public:
     explicit VectorMatrixAssembler(
         std::unique_ptr<AbstractJacobianAssembler>&& jacobian_assembler);
 
+    void preAssemble(const std::size_t mesh_item_id,
+                     LocalAssemblerInterface& local_assembler,
+                     const NumLib::LocalToGlobalIndexMap& dof_table,
+                     const double t, const GlobalVector& x);
+
     //! Assembles\c M, \c K, and \c b.
     //! \remark Jacobian is not assembled here, see assembleWithJacobian().
     void assemble(std::size_t const mesh_item_id,

Tests/NumLib/ODEs.h

@@ -26,6 +26,8 @@ class ODE1 final : public NumLib::ODESystem<
                        NumLib::NonlinearSolverTag::Newton>
 {
 public:
+    void preAssemble(const double t, GlobalVector const& x) override {}
+
     void assemble(const double /*t*/, GlobalVector const& /*x*/,
                   GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override
     {
@@ -102,6 +104,8 @@ class ODE2 final : public NumLib::ODESystem<
                        NumLib::NonlinearSolverTag::Newton>
 {
 public:
+    void preAssemble(const double t, GlobalVector const& x) override {}
+
     void assemble(const double /*t*/, GlobalVector const& x, GlobalMatrix& M,
                   GlobalMatrix& K, GlobalVector& b) override
     {
@@ -187,6 +191,8 @@ class ODE3 final : public NumLib::ODESystem<
                        NumLib::NonlinearSolverTag::Newton>
 {
 public:
+    void preAssemble(const double t, GlobalVector const& x) override {}
+
     void assemble(const double /*t*/, GlobalVector const& x_curr, GlobalMatrix& M,
                   GlobalMatrix& K, GlobalVector& b) override
     {