From 0be2ee6810a7024f8e411cf65a5ccf805b8a6111 Mon Sep 17 00:00:00 2001
From: Francesco Parisio <francesco.parisio@protonmail.com>
Date: Fri, 10 Feb 2017 12:11:23 +0100
Subject: [PATCH] Add ode system preAssemble() call.
In some processes (using non-local methods, for example)
the assembly is split in two parts: computation of local
quantities and actual (non-local) assembly (integration)
into the local matrices.
---
NumLib/ODESolver/ODESystem.h | 7 +++++++
NumLib/ODESolver/TimeDiscretizedODESystem.cpp | 2 ++
ProcessLib/LocalAssemblerInterface.h | 3 +++
ProcessLib/Process.cpp | 5 +++++
ProcessLib/Process.h | 7 +++++++
ProcessLib/VectorMatrixAssembler.cpp | 11 +++++++++++
ProcessLib/VectorMatrixAssembler.h | 5 +++++
Tests/NumLib/ODEs.h | 6 ++++++
8 files changed, 46 insertions(+)
diff --git a/NumLib/ODESolver/ODESystem.h b/NumLib/ODESolver/ODESystem.h
index b362403f9b5..7dd004d058d 100644
--- a/NumLib/ODESolver/ODESystem.h
+++ b/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).
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index a1b6f29685c..064c3d40c56 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/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);
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index bfb983f87d4..e3fc5cb7983 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/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,
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 808bb909c59..119e9063dd9 100644
--- a/ProcessLib/Process.cpp
+++ b/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)
{
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index d067d0ff747..ef2ef4d1b09 100644
--- a/ProcessLib/Process.h
+++ 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;
diff --git a/ProcessLib/VectorMatrixAssembler.cpp b/ProcessLib/VectorMatrixAssembler.cpp
index 398c7a5d45b..690c9bcc23c 100644
--- a/ProcessLib/VectorMatrixAssembler.cpp
+++ b/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,
diff --git a/ProcessLib/VectorMatrixAssembler.h b/ProcessLib/VectorMatrixAssembler.h
index 5f305cc482b..a1203361b4d 100644
--- a/ProcessLib/VectorMatrixAssembler.h
+++ b/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,
diff --git a/Tests/NumLib/ODEs.h b/Tests/NumLib/ODEs.h
index f5326f712ff..fbe3f271563 100644
--- a/Tests/NumLib/ODEs.h
+++ b/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
{
--
GitLab