diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 58c6ff0169443a8a9eceece6b25d7dee7db5017d..ad78dd94513b33841626e7f1ac5ee906af3dd2d2 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -330,8 +330,7 @@ void ProjectData::parseTimeStepping(BaseLib::ConfigTree const& timestepping_conf
{
DBUG("Reading time loop configuration.");
- _time_loop = ApplicationsLib::createUncoupledProcessesTimeLoop<
- GlobalMatrix, GlobalVector>(timestepping_config);
+ _time_loop = ApplicationsLib::createUncoupledProcessesTimeLoop(timestepping_config);
if (!_time_loop)
{
@@ -372,7 +371,7 @@ void ProjectData::parseNonlinearSolvers(BaseLib::ConfigTree const& config)
auto const name = conf.getConfigParameter<std::string>("name");
BaseLib::insertIfKeyUniqueElseError(_nonlinear_solvers,
name,
- NumLib::createNonlinearSolver<GlobalMatrix, GlobalVector>(
+ NumLib::createNonlinearSolver(
*linear_solver, conf).first,
"The nonlinear solver name is not unique");
}
diff --git a/Applications/ApplicationsLib/ProjectData.h b/Applications/ApplicationsLib/ProjectData.h
index 9d5abda1704d66154d445d199698f6fbe3ab79b9..9397f5a75e6d918fa4fad2eeb22af6e0006b35ce 100644
--- a/Applications/ApplicationsLib/ProjectData.h
+++ b/Applications/ApplicationsLib/ProjectData.h
@@ -44,7 +44,6 @@ class NonlinearSolverBase;
namespace ApplicationsLib
{
-template<typename Matrix, typename Vector>
class UncoupledProcessesTimeLoop;
}
@@ -57,8 +56,7 @@ class ProjectData final
{
public:
/// The time loop type used to solve this project's processes.
- using TimeLoop = ApplicationsLib::UncoupledProcessesTimeLoop<
- GlobalMatrix, GlobalVector>;
+ using TimeLoop = ApplicationsLib::UncoupledProcessesTimeLoop;
/// The empty constructor used in the gui, for example, when the project's
/// configuration is not loaded yet.
diff --git a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ff891e6503b96805515dbeb94f5a5d6884f703fd
--- /dev/null
+++ b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp
@@ -0,0 +1,229 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "UncoupledProcessesTimeLoop.h"
+
+namespace ApplicationsLib
+{
+
+std::unique_ptr<UncoupledProcessesTimeLoop>
+createUncoupledProcessesTimeLoop(BaseLib::ConfigTree const& conf)
+{
+ //! \ogs_file_param{prj__time_stepping__type}
+ auto const type = conf.peekConfigParameter<std::string>("type");
+
+ std::unique_ptr<NumLib::ITimeStepAlgorithm> timestepper;
+
+ if (type == "SingleStep") {
+ conf.ignoreConfigParameter("type");
+ timestepper.reset(new NumLib::FixedTimeStepping(0.0, 1.0, 1.0));
+ } else if (type == "FixedTimeStepping") {
+ timestepper = NumLib::FixedTimeStepping::newInstance(conf);
+ } else {
+ OGS_FATAL("Unknown timestepper type: `%s'.", type.c_str());
+ }
+
+ using TimeLoop = UncoupledProcessesTimeLoop;
+ return std::unique_ptr<TimeLoop>{new TimeLoop{std::move(timestepper)}};
+}
+
+
+std::vector<typename UncoupledProcessesTimeLoop::SingleProcessData>
+UncoupledProcessesTimeLoop::
+initInternalData(ProjectData& project)
+{
+ auto const num_processes = std::distance(project.processesBegin(),
+ project.processesEnd());
+
+ std::vector<SingleProcessData> per_process_data;
+ per_process_data.reserve(num_processes);
+
+ // create a time discretized ODE system for each process
+ for (auto p = project.processesBegin(); p != project.processesEnd(); ++p)
+ {
+ auto& pcs = **p;
+ auto& nonlinear_solver = pcs.getNonlinearSolver();
+ auto& time_disc = pcs.getTimeDiscretization();
+
+ per_process_data.emplace_back(
+ makeSingleProcessData(nonlinear_solver, **p, time_disc));
+ }
+
+ return per_process_data;
+}
+
+
+void
+UncoupledProcessesTimeLoop::
+setInitialConditions(ProjectData& project,
+ double const t0,
+ std::vector<SingleProcessData>& per_process_data)
+{
+ auto const num_processes = std::distance(project.processesBegin(),
+ project.processesEnd());
+
+ _process_solutions.reserve(num_processes);
+
+ unsigned pcs_idx = 0;
+ for (auto p = project.processesBegin(); p != project.processesEnd();
+ ++p, ++pcs_idx)
+ {
+ auto& pcs = **p;
+ auto& time_disc = pcs.getTimeDiscretization();
+
+ auto& ppd = per_process_data[pcs_idx];
+ auto& ode_sys = *ppd.tdisc_ode_sys;
+ auto const nl_tag = ppd.nonlinear_solver_tag;
+
+ // append a solution vector of suitable size
+ _process_solutions.emplace_back(
+ &MathLib::GlobalVectorProvider<GlobalVector>::provider.getVector(
+ ode_sys.getMatrixSpecifications()));
+
+ auto& x0 = *_process_solutions[pcs_idx];
+ pcs.setInitialConditions(x0);
+ MathLib::BLAS::finalizeAssembly(x0);
+
+ time_disc.setInitialState(t0, x0); // push IC
+
+ if (time_disc.needsPreload())
+ {
+ auto& nonlinear_solver = ppd.nonlinear_solver;
+ auto& mat_strg = ppd.mat_strg;
+
+ setEquationSystem(nonlinear_solver, ode_sys, nl_tag);
+ nonlinear_solver.assemble(x0);
+ time_disc.pushState(t0, x0, mat_strg); // TODO: that might do duplicate work
+ }
+ }
+}
+
+
+bool
+UncoupledProcessesTimeLoop::
+solveOneTimeStepOneProcess(
+ GlobalVector& x, double const t, double const delta_t,
+ SingleProcessData& process_data,
+ UncoupledProcessesTimeLoop::Process& process)
+{
+ auto& time_disc = process.getTimeDiscretization();
+ auto& ode_sys = *process_data.tdisc_ode_sys;
+ auto& nonlinear_solver = process_data.nonlinear_solver;
+ auto const nl_tag = process_data.nonlinear_solver_tag;
+
+ setEquationSystem(nonlinear_solver, ode_sys, nl_tag);
+
+ // Note: Order matters!
+ // First advance to the next timestep, then set known solutions at that
+ // time, afterwards pass the right solution vector and time to the
+ // preTimestep() hook.
+
+ time_disc.nextTimestep(t, delta_t);
+
+ applyKnownSolutions(ode_sys, nl_tag, x);
+
+ process.preTimestep(x, t, delta_t);
+
+ bool nonlinear_solver_succeeded = nonlinear_solver.solve(x);
+
+ auto& mat_strg = process_data.mat_strg;
+ time_disc.pushState(t, x, mat_strg);
+
+ process.postTimestep(x);
+
+ return nonlinear_solver_succeeded;
+}
+
+
+bool UncoupledProcessesTimeLoop::loop(ProjectData& project)
+{
+ auto per_process_data = initInternalData(project);
+
+ auto& out_ctrl = project.getOutputControl();
+ out_ctrl.initialize(project.processesBegin(), project.processesEnd());
+
+ auto const t0 = _timestepper->getTimeStep().current(); // time of the IC
+
+ // init solution storage
+ setInitialConditions(project, t0, per_process_data);
+
+ // output initial conditions
+ {
+ unsigned pcs_idx = 0;
+ for (auto p = project.processesBegin(); p != project.processesEnd();
+ ++p, ++pcs_idx)
+ {
+ auto const& x0 = *_process_solutions[pcs_idx];
+ out_ctrl.doOutput(**p, 0, t0, x0);
+ }
+ }
+
+ double t = t0;
+ std::size_t timestep = 1; // the first timestep really is number one
+ bool nonlinear_solver_succeeded = true;
+
+ while (_timestepper->next())
+ {
+ auto const ts = _timestepper->getTimeStep();
+ auto const delta_t = ts.dt();
+ t = ts.current();
+ timestep = ts.steps();
+
+ INFO("=== timestep #%u (t=%gs, dt=%gs) ==============================",
+ timestep, t, delta_t);
+
+ // TODO use process name
+ unsigned pcs_idx = 0;
+ for (auto p = project.processesBegin(); p != project.processesEnd();
+ ++p, ++pcs_idx)
+ {
+ auto& x = *_process_solutions[pcs_idx];
+
+ nonlinear_solver_succeeded = solveOneTimeStepOneProcess(
+ x, t, delta_t, per_process_data[pcs_idx], **p);
+
+ if (!nonlinear_solver_succeeded) {
+ ERR("The nonlinear solver failed in timestep #%u at t = %g s"
+ " for process #%u.", timestep, t, pcs_idx);
+
+ // save unsuccessful solution
+ out_ctrl.doOutputAlways(**p, timestep, t, x);
+
+ break;
+ } else {
+ out_ctrl.doOutput(**p, timestep, t, x);
+ }
+ }
+
+ if (!nonlinear_solver_succeeded) break;
+ }
+
+ // output last timestep
+ if (nonlinear_solver_succeeded)
+ {
+ unsigned pcs_idx = 0;
+ for (auto p = project.processesBegin(); p != project.processesEnd();
+ ++p, ++pcs_idx)
+ {
+ auto const& x = *_process_solutions[pcs_idx];
+ out_ctrl.doOutputLastTimestep(**p, timestep, t, x);
+ }
+ }
+
+ return nonlinear_solver_succeeded;
+}
+
+
+UncoupledProcessesTimeLoop::~UncoupledProcessesTimeLoop()
+{
+ for (auto * x : _process_solutions)
+ MathLib::GlobalVectorProvider<GlobalVector>::provider.releaseVector(*x);
+}
+
+} // namespace ApplicationsLib
diff --git a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h
index 2ce2b8a2c37dbc30ad6293408b8bf1ddd0fa5db3..a89d8d15826e6f1b28df181a68ef95cb01387660 100644
--- a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h
+++ b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h
@@ -26,7 +26,6 @@ namespace ApplicationsLib
{
//! Time loop capable of time-integrating several uncoupled processes at once.
-template<typename Matrix, typename Vector>
class UncoupledProcessesTimeLoop
{
public:
@@ -41,28 +40,28 @@ public:
private:
//! An abstract nonlinear solver
- using AbstractNLSolver = NumLib::NonlinearSolverBase<Matrix, Vector>;
+ using AbstractNLSolver = NumLib::NonlinearSolverBase<GlobalMatrix, GlobalVector>;
//! An abstract equations system
- using EquationSystem = NumLib::EquationSystem<Vector>;
+ using EquationSystem = NumLib::EquationSystem<GlobalVector>;
//! An abstract process
using Process = ProcessLib::Process;
//! An abstract time discretization
- using TimeDisc = NumLib::TimeDiscretization<Vector>;
+ using TimeDisc = NumLib::TimeDiscretization<GlobalVector>;
- std::vector<Vector*> _process_solutions;
+ std::vector<GlobalVector*> _process_solutions;
std::unique_ptr<NumLib::ITimeStepAlgorithm> _timestepper;
struct SingleProcessData
{
template<NumLib::ODESystemTag ODETag, NumLib::NonlinearSolverTag NLTag>
SingleProcessData(
- NumLib::NonlinearSolver<Matrix, Vector, NLTag>& nonlinear_solver,
+ NumLib::NonlinearSolver<GlobalMatrix, GlobalVector, NLTag>& nonlinear_solver,
TimeDisc& time_disc,
- NumLib::ODESystem<Matrix, Vector, ODETag, NLTag>& ode_sys)
+ NumLib::ODESystem<GlobalMatrix, GlobalVector, ODETag, NLTag>& ode_sys)
: nonlinear_solver_tag(NLTag)
, nonlinear_solver(nonlinear_solver)
, tdisc_ode_sys(
- new NumLib::TimeDiscretizedODESystem<Matrix, Vector, ODETag, NLTag>(
+ new NumLib::TimeDiscretizedODESystem<GlobalMatrix, GlobalVector, ODETag, NLTag>(
ode_sys, time_disc))
, mat_strg(dynamic_cast<NumLib::InternalMatrixStorage&>(*tdisc_ode_sys))
{}
@@ -94,7 +93,7 @@ private:
*
* \note Currently the \c ODETag is automatically inferred from the given
* \c ody_sys. This works as long as \c Process derives from
- * \c ODESystem<Matrix, Vector, ODETag>, i.e. as long we only deal with
+ * \c ODESystem<GlobalMatrix, GlobalVector, ODETag>, i.e. as long we only deal with
* one type of ODE. When we introduce more types, this method will have
* to be extended slightly.
*/
@@ -102,17 +101,17 @@ private:
SingleProcessData
makeSingleProcessData(
AbstractNLSolver& nonlinear_solver,
- NumLib::ODESystem<Matrix, Vector, ODETag,
+ NumLib::ODESystem<GlobalMatrix, GlobalVector, ODETag,
NumLib::NonlinearSolverTag::Picard>& ode_sys,
TimeDisc& time_disc)
{
using Tag = NumLib::NonlinearSolverTag;
// A concrete Picard solver
using NonlinearSolverPicard =
- NumLib::NonlinearSolver<Matrix, Vector, Tag::Picard>;
+ NumLib::NonlinearSolver<GlobalMatrix, GlobalVector, Tag::Picard>;
// A concrete Newton solver
using NonlinearSolverNewton =
- NumLib::NonlinearSolver<Matrix, Vector, Tag::Newton>;
+ NumLib::NonlinearSolver<GlobalMatrix, GlobalVector, Tag::Newton>;
if (auto* nonlinear_solver_picard =
dynamic_cast<NonlinearSolverPicard*>(&nonlinear_solver))
@@ -129,7 +128,7 @@ private:
{
// The Newton-Raphson method needs a Newton-ready ODE.
- using ODENewton = NumLib::ODESystem<Matrix, Vector, ODETag, Tag::Newton>;
+ using ODENewton = NumLib::ODESystem<GlobalMatrix, GlobalVector, ODETag, Tag::Newton>;
if (auto* ode_newton = dynamic_cast<ODENewton*>(&ode_sys))
{
return SingleProcessData{
@@ -156,7 +155,7 @@ private:
//! Solves one timestep for the given \c process.
bool solveOneTimeStepOneProcess(
- Vector& x, double const t, double const delta_t,
+ GlobalVector& x, double const t, double const delta_t,
SingleProcessData& process_data,
Process& process);
@@ -166,8 +165,8 @@ private:
static void setEquationSystem(AbstractNLSolver& nonlinear_solver,
EquationSystem& eq_sys)
{
- using Solver = NumLib::NonlinearSolver<Matrix, Vector, NLTag>;
- using EqSys = NumLib::NonlinearSystem<Matrix, Vector, NLTag>;
+ using Solver = NumLib::NonlinearSolver<GlobalMatrix, GlobalVector, NLTag>;
+ using EqSys = NumLib::NonlinearSystem<GlobalMatrix, GlobalVector, NLTag>;
assert(dynamic_cast<Solver*>(&nonlinear_solver) != nullptr);
assert(dynamic_cast<EqSys*> (&eq_sys) != nullptr);
@@ -200,9 +199,9 @@ private:
//! for equation systems linearized with either the Picard or Newton method.
template<NumLib::NonlinearSolverTag NLTag>
static void applyKnownSolutions(
- EquationSystem const& eq_sys, Vector& x)
+ EquationSystem const& eq_sys, GlobalVector& x)
{
- using EqSys = NumLib::NonlinearSystem<Matrix, Vector, NLTag>;
+ using EqSys = NumLib::NonlinearSystem<GlobalMatrix, GlobalVector, NLTag>;
assert(dynamic_cast<EqSys const*> (&eq_sys) != nullptr);
auto& eq_sys_ = static_cast<EqSys const&> (eq_sys);
@@ -213,7 +212,7 @@ private:
//! for equation systems linearized with either the Picard or Newton method.
static void applyKnownSolutions(
EquationSystem const& eq_sys,
- NumLib::NonlinearSolverTag const nl_tag, Vector& x)
+ NumLib::NonlinearSolverTag const nl_tag, GlobalVector& x)
{
using Tag = NumLib::NonlinearSolverTag;
switch (nl_tag)
@@ -225,228 +224,8 @@ private:
};
//! Builds an UncoupledProcessesTimeLoop from the given configuration.
-template<typename Matrix, typename Vector>
-std::unique_ptr<UncoupledProcessesTimeLoop<Matrix, Vector> >
-createUncoupledProcessesTimeLoop(BaseLib::ConfigTree const& conf)
-{
- //! \ogs_file_param{prj__time_stepping__type}
- auto const type = conf.peekConfigParameter<std::string>("type");
-
- std::unique_ptr<NumLib::ITimeStepAlgorithm> timestepper;
-
- if (type == "SingleStep") {
- conf.ignoreConfigParameter("type");
- timestepper.reset(new NumLib::FixedTimeStepping(0.0, 1.0, 1.0));
- } else if (type == "FixedTimeStepping") {
- timestepper = NumLib::FixedTimeStepping::newInstance(conf);
- } else {
- OGS_FATAL("Unknown timestepper type: `%s'.", type.c_str());
- }
-
- using TimeLoop = UncoupledProcessesTimeLoop<Matrix, Vector>;
- return std::unique_ptr<TimeLoop>{new TimeLoop{std::move(timestepper)}};
-}
-
-
-template<typename Matrix, typename Vector>
-std::vector<typename UncoupledProcessesTimeLoop<Matrix, Vector>::SingleProcessData>
-UncoupledProcessesTimeLoop<Matrix, Vector>::
-initInternalData(ProjectData& project)
-{
- auto const num_processes = std::distance(project.processesBegin(),
- project.processesEnd());
-
- std::vector<SingleProcessData> per_process_data;
- per_process_data.reserve(num_processes);
-
- // create a time discretized ODE system for each process
- for (auto p = project.processesBegin(); p != project.processesEnd(); ++p)
- {
- auto& pcs = **p;
- auto& nonlinear_solver = pcs.getNonlinearSolver();
- auto& time_disc = pcs.getTimeDiscretization();
-
- per_process_data.emplace_back(
- makeSingleProcessData(nonlinear_solver, **p, time_disc));
- }
-
- return per_process_data;
-}
-
-
-template<typename Matrix, typename Vector>
-void
-UncoupledProcessesTimeLoop<Matrix, Vector>::
-setInitialConditions(ProjectData& project,
- double const t0,
- std::vector<SingleProcessData>& per_process_data)
-{
- auto const num_processes = std::distance(project.processesBegin(),
- project.processesEnd());
-
- _process_solutions.reserve(num_processes);
-
- unsigned pcs_idx = 0;
- for (auto p = project.processesBegin(); p != project.processesEnd();
- ++p, ++pcs_idx)
- {
- auto& pcs = **p;
- auto& time_disc = pcs.getTimeDiscretization();
-
- auto& ppd = per_process_data[pcs_idx];
- auto& ode_sys = *ppd.tdisc_ode_sys;
- auto const nl_tag = ppd.nonlinear_solver_tag;
-
- // append a solution vector of suitable size
- _process_solutions.emplace_back(
- &MathLib::GlobalVectorProvider<Vector>::provider.getVector(
- ode_sys.getMatrixSpecifications()));
-
- auto& x0 = *_process_solutions[pcs_idx];
- pcs.setInitialConditions(x0);
- MathLib::BLAS::finalizeAssembly(x0);
-
- time_disc.setInitialState(t0, x0); // push IC
-
- if (time_disc.needsPreload())
- {
- auto& nonlinear_solver = ppd.nonlinear_solver;
- auto& mat_strg = ppd.mat_strg;
-
- setEquationSystem(nonlinear_solver, ode_sys, nl_tag);
- nonlinear_solver.assemble(x0);
- time_disc.pushState(t0, x0, mat_strg); // TODO: that might do duplicate work
- }
- }
-}
-
-
-template<typename Matrix, typename Vector>
-bool
-UncoupledProcessesTimeLoop<Matrix, Vector>::
-solveOneTimeStepOneProcess(
- Vector& x, double const t, double const delta_t,
- SingleProcessData& process_data,
- typename UncoupledProcessesTimeLoop<Matrix, Vector>::Process& process)
-{
- auto& time_disc = process.getTimeDiscretization();
- auto& ode_sys = *process_data.tdisc_ode_sys;
- auto& nonlinear_solver = process_data.nonlinear_solver;
- auto const nl_tag = process_data.nonlinear_solver_tag;
-
- setEquationSystem(nonlinear_solver, ode_sys, nl_tag);
-
- // Note: Order matters!
- // First advance to the next timestep, then set known solutions at that
- // time, afterwards pass the right solution vector and time to the
- // preTimestep() hook.
-
- time_disc.nextTimestep(t, delta_t);
-
- applyKnownSolutions(ode_sys, nl_tag, x);
-
- process.preTimestep(x, t, delta_t);
-
- bool nonlinear_solver_succeeded = nonlinear_solver.solve(x);
-
- auto& mat_strg = process_data.mat_strg;
- time_disc.pushState(t, x, mat_strg);
-
- process.postTimestep(x);
-
- return nonlinear_solver_succeeded;
-}
-
-
-template<typename Matrix, typename Vector>
-bool
-UncoupledProcessesTimeLoop<Matrix, Vector>::
-loop(ProjectData& project)
-{
- auto per_process_data = initInternalData(project);
-
- auto& out_ctrl = project.getOutputControl();
- out_ctrl.initialize(project.processesBegin(), project.processesEnd());
-
- auto const t0 = _timestepper->getTimeStep().current(); // time of the IC
-
- // init solution storage
- setInitialConditions(project, t0, per_process_data);
-
- // output initial conditions
- {
- unsigned pcs_idx = 0;
- for (auto p = project.processesBegin(); p != project.processesEnd();
- ++p, ++pcs_idx)
- {
- auto const& x0 = *_process_solutions[pcs_idx];
- out_ctrl.doOutput(**p, 0, t0, x0);
- }
- }
-
- double t = t0;
- std::size_t timestep = 1; // the first timestep really is number one
- bool nonlinear_solver_succeeded = true;
-
- while (_timestepper->next())
- {
- auto const ts = _timestepper->getTimeStep();
- auto const delta_t = ts.dt();
- t = ts.current();
- timestep = ts.steps();
-
- INFO("=== timestep #%u (t=%gs, dt=%gs) ==============================",
- timestep, t, delta_t);
-
- // TODO use process name
- unsigned pcs_idx = 0;
- for (auto p = project.processesBegin(); p != project.processesEnd();
- ++p, ++pcs_idx)
- {
- auto& x = *_process_solutions[pcs_idx];
-
- nonlinear_solver_succeeded = solveOneTimeStepOneProcess(
- x, t, delta_t, per_process_data[pcs_idx], **p);
-
- if (!nonlinear_solver_succeeded) {
- ERR("The nonlinear solver failed in timestep #%u at t = %g s"
- " for process #%u.", timestep, t, pcs_idx);
-
- // save unsuccessful solution
- out_ctrl.doOutputAlways(**p, timestep, t, x);
-
- break;
- } else {
- out_ctrl.doOutput(**p, timestep, t, x);
- }
- }
-
- if (!nonlinear_solver_succeeded) break;
- }
-
- // output last timestep
- if (nonlinear_solver_succeeded)
- {
- unsigned pcs_idx = 0;
- for (auto p = project.processesBegin(); p != project.processesEnd();
- ++p, ++pcs_idx)
- {
- auto const& x = *_process_solutions[pcs_idx];
- out_ctrl.doOutputLastTimestep(**p, timestep, t, x);
- }
- }
-
- return nonlinear_solver_succeeded;
-}
-
-
-template<typename Matrix, typename Vector>
-UncoupledProcessesTimeLoop<Matrix, Vector>::
-~UncoupledProcessesTimeLoop()
-{
- for (auto * x : _process_solutions)
- MathLib::GlobalVectorProvider<Vector>::provider.releaseVector(*x);
-}
+std::unique_ptr<UncoupledProcessesTimeLoop>
+createUncoupledProcessesTimeLoop(BaseLib::ConfigTree const& conf);
} // namespace ApplicationsLib
diff --git a/NumLib/Assembler/VectorMatrixAssembler.h b/NumLib/Assembler/VectorMatrixAssembler.h
index 1a0433a8985f674a5250f6b0645844a97883453b..99e682e32c1b47f02cb0287c172b13a02fd839f0 100644
--- a/NumLib/Assembler/VectorMatrixAssembler.h
+++ b/NumLib/Assembler/VectorMatrixAssembler.h
@@ -37,8 +37,7 @@ getIndices(std::size_t const id,
return indices;
}
-template <typename GlobalVector>
-std::vector<double>
+inline std::vector<double>
getLocalNodalDOFs(GlobalVector const& x,
std::vector<GlobalIndexType> const& dof_indices)
{
@@ -64,15 +63,14 @@ getLocalNodalDOFs(GlobalVector const& x,
* Each type of equation as flagged by the \c ODETag will have a different
* VectorMatrixAssembler type.
*/
-template<typename GlobalMatrix, typename GlobalVector,
- typename LocalAssembler,
+template<typename LocalAssembler,
NumLib::ODESystemTag ODETag>
class VectorMatrixAssembler;
//! Specialization for first-order implicit quasi-linear systems.
-template<typename GlobalMatrix, typename GlobalVector, typename LocalAssembler>
-class VectorMatrixAssembler<GlobalMatrix, GlobalVector, LocalAssembler,
+template<typename LocalAssembler>
+class VectorMatrixAssembler<LocalAssembler,
NumLib::ODESystemTag::FirstOrderImplicitQuasilinear> final
{
public:
@@ -167,8 +165,8 @@ private:
//! Specialization used to add Neumann boundary conditions.
-template<typename GlobalMatrix, typename GlobalVector, typename LocalAssembler>
-class VectorMatrixAssembler<GlobalMatrix, GlobalVector, LocalAssembler,
+template<typename LocalAssembler>
+class VectorMatrixAssembler<LocalAssembler,
NumLib::ODESystemTag::NeumannBC> final
{
public:
diff --git a/NumLib/Extrapolation/Extrapolator.h b/NumLib/Extrapolation/Extrapolator.h
index 797012aab204eaf33e53bc226da1eaf4d03dcd7e..0867cfa389a937efd852c24313d88d3516e5e823 100644
--- a/NumLib/Extrapolation/Extrapolator.h
+++ b/NumLib/Extrapolation/Extrapolator.h
@@ -12,6 +12,9 @@
#include <vector>
+#include <Eigen/Eigen>
+#include "NumLib/NumericsConfig.h"
+
namespace NumLib
{
@@ -21,11 +24,10 @@ namespace NumLib
* Local assemblers that want to have some integration point values extrapolated
* using Extrapolator have to implement this interface.
*
- * \tparam GlobalVector type of the global vector
* \tparam PropertyTag type of the property used to query a specific kind of
* integration point value, usually an enum.
*/
-template <typename GlobalVector, typename PropertyTag>
+template <typename PropertyTag>
class Extrapolatable
{
public:
@@ -55,13 +57,12 @@ public:
/*! Interface for classes that extrapolate integration point values to nodal
* values.
*
- * \tparam GlobalVector type of the global vector
* \tparam PropertyTag type of the property used to query a specific kind of
* integration point value, usually an enum.
* \tparam LocalAssembler type of the local assembler being queried for
* integration point values.
*/
-template<typename GlobalVector, typename PropertyTag, typename LocalAssembler>
+template<typename PropertyTag, typename LocalAssembler>
class Extrapolator
{
public:
diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h
index 9215984c16f18c4b77dea9e879bdc4287b16537b..b95f39cb29127f7ebb503293b209ea1ce25a664a 100644
--- a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h
+++ b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h
@@ -21,9 +21,9 @@
namespace NumLib
{
-template<typename GlobalVector, typename PropertyTag, typename LocalAssembler>
+template<typename PropertyTag, typename LocalAssembler>
void
-LocalLinearLeastSquaresExtrapolator<GlobalVector, PropertyTag, LocalAssembler>::
+LocalLinearLeastSquaresExtrapolator<PropertyTag, LocalAssembler>::
extrapolate(LocalAssemblers const& local_assemblers, PropertyTag const property)
{
_nodal_values.setZero();
@@ -35,7 +35,7 @@ extrapolate(LocalAssemblers const& local_assemblers, PropertyTag const property)
counts->setZero(); // TODO BLAS?
using Self = LocalLinearLeastSquaresExtrapolator<
- GlobalVector, PropertyTag, LocalAssembler>;
+ PropertyTag, LocalAssembler>;
NumLib::SerialExecutor::executeMemberDereferenced(
*this, &Self::extrapolateElement, local_assemblers, property, *counts);
@@ -43,24 +43,24 @@ extrapolate(LocalAssemblers const& local_assemblers, PropertyTag const property)
MathLib::BLAS::componentwiseDivide(_nodal_values, _nodal_values, *counts);
}
-template<typename GlobalVector, typename PropertyTag, typename LocalAssembler>
+template<typename PropertyTag, typename LocalAssembler>
void
-LocalLinearLeastSquaresExtrapolator<GlobalVector, PropertyTag, LocalAssembler>::
+LocalLinearLeastSquaresExtrapolator<PropertyTag, LocalAssembler>::
calculateResiduals(LocalAssemblers const& local_assemblers,
PropertyTag const property)
{
assert(static_cast<std::size_t>(_residuals.size()) == local_assemblers.size());
using Self = LocalLinearLeastSquaresExtrapolator<
- GlobalVector, PropertyTag, LocalAssembler>;
+ PropertyTag, LocalAssembler>;
NumLib::SerialExecutor::executeMemberDereferenced(
*this, &Self::calculateResiudalElement, local_assemblers, property);
}
-template<typename GlobalVector, typename PropertyTag, typename LocalAssembler>
+template<typename PropertyTag, typename LocalAssembler>
void
-LocalLinearLeastSquaresExtrapolator<GlobalVector, PropertyTag, LocalAssembler>::
+LocalLinearLeastSquaresExtrapolator<PropertyTag, LocalAssembler>::
extrapolateElement(std::size_t const element_index,
LocalAssembler const& loc_asm, PropertyTag const property,
GlobalVector& counts)
@@ -115,9 +115,9 @@ extrapolateElement(std::size_t const element_index,
counts.add(global_indices, std::vector<double>(global_indices.size(), 1.0));
}
-template<typename GlobalVector, typename PropertyTag, typename LocalAssembler>
+template<typename PropertyTag, typename LocalAssembler>
void
-LocalLinearLeastSquaresExtrapolator<GlobalVector, PropertyTag, LocalAssembler>::
+LocalLinearLeastSquaresExtrapolator<PropertyTag, LocalAssembler>::
calculateResiudalElement(std::size_t const element_index,
LocalAssembler const& loc_asm, PropertyTag const property)
{
diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
index 6fd3066409dfd33c1cd7d435c95da7c9ae07c739..1336a25994f71c676e5ceaf3638ccefec46c30f0 100644
--- a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
+++ b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
@@ -10,6 +10,7 @@
#ifndef NUMLIB_LOCAL_LLSQ_EXTRAPOLATOR_H
#define NUMLIB_LOCAL_LLSQ_EXTRAPOLATOR_H
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
#include "NumLib/DOF/GlobalMatrixProviders.h"
#include "Extrapolator.h"
@@ -34,13 +35,13 @@ namespace NumLib
* use of the least squares which requires an exact or overdetermined equation system.
* \endparblock
*/
-template<typename GlobalVector, typename PropertyTag, typename LocalAssembler>
+template<typename PropertyTag, typename LocalAssembler>
class LocalLinearLeastSquaresExtrapolator
- : public Extrapolator<GlobalVector, PropertyTag, LocalAssembler>
+ : public Extrapolator<PropertyTag, LocalAssembler>
{
public:
using LocalAssemblers = typename Extrapolator<
- GlobalVector, PropertyTag, LocalAssembler>::LocalAssemblers;
+ PropertyTag, LocalAssembler>::LocalAssemblers;
/*! Constructs a new instance
*
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
index aa2d57f99cf2967f93304e4a91d65b7bc125b578..4f21fdc34310f2970b1c79d8c58258f09c26c35a 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
@@ -34,19 +34,16 @@ enum class IntegrationPointValue {
const unsigned NUM_NODAL_DOF = 1;
-template <typename GlobalMatrix, typename GlobalVector>
class GroundwaterFlowLocalAssemblerInterface
- : public ProcessLib::LocalAssemblerInterface<GlobalMatrix, GlobalVector>
- , public NumLib::Extrapolatable<GlobalVector, IntegrationPointValue>
+ : public ProcessLib::LocalAssemblerInterface
+ , public NumLib::Extrapolatable<IntegrationPointValue>
{};
template <typename ShapeFunction,
typename IntegrationMethod,
- typename GlobalMatrix,
- typename GlobalVector,
unsigned GlobalDim>
class LocalAssemblerData
- : public GroundwaterFlowLocalAssemblerInterface<GlobalMatrix, GlobalVector>
+ : public GroundwaterFlowLocalAssemblerInterface
{
using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1c86cc3642a1c3e4e8bbdcfa8a0a5c89dbcd096f
--- /dev/null
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
@@ -0,0 +1,164 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "GroundwaterFlowProcess.h"
+
+#include <cassert>
+
+#include "NumLib/Assembler/VectorMatrixAssembler.h"
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+
+
+namespace ProcessLib
+{
+namespace GroundwaterFlow
+{
+
+GroundwaterFlowProcess::GroundwaterFlowProcess(
+ MeshLib::Mesh& mesh,
+ Base::NonlinearSolver& nonlinear_solver,
+ std::unique_ptr<Base::TimeDiscretization>&& time_discretization,
+ std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
+ GroundwaterFlowProcessData&& process_data,
+ SecondaryVariableCollection&& secondary_variables,
+ ProcessOutput&& process_output
+ )
+ : Process(mesh, nonlinear_solver, std::move(time_discretization),
+ std::move(process_variables),
+ std::move(secondary_variables),
+ std::move(process_output))
+ , _process_data(std::move(process_data))
+{
+ if (dynamic_cast<NumLib::ForwardEuler<GlobalVector>*>(
+ &Base::getTimeDiscretization()) != nullptr)
+ {
+ OGS_FATAL("GroundwaterFlowProcess can not be solved with the ForwardEuler"
+ " time discretization scheme. Aborting");
+ // Because the M matrix is not assembled. Thus, the linearized system
+ // would be singular. The same applies to CrankNicolson with theta = 0.0,
+ // but this case is not checked here.
+ // Anyway, the GroundwaterFlowProcess shall be transferred to a simpler
+ // ODESystemTag in the future.
+ }
+}
+
+void GroundwaterFlowProcess::initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order)
+{
+ DBUG("Create global assembler.");
+ _global_assembler.reset(new GlobalAssembler(dof_table));
+
+ ProcessLib::createLocalAssemblers<LocalAssemblerData>(
+ mesh.getDimension(), mesh.getElements(),
+ dof_table, integration_order, _local_assemblers,
+ _process_data);
+
+ // TOOD Later on the DOF table can change during the simulation!
+ _extrapolator.reset(
+ new ExtrapolatorImplementation(Base::getMatrixSpecifications(),
+ *Base::_local_to_global_index_map));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "darcy_velocity_x", 1,
+ makeExtrapolator(IntegrationPointValue::DarcyVelocityX, *_extrapolator,
+ _local_assemblers));
+
+ if (mesh.getDimension() > 1) {
+ Base::_secondary_variables.addSecondaryVariable(
+ "darcy_velocity_y", 1,
+ makeExtrapolator(IntegrationPointValue::DarcyVelocityY, *_extrapolator,
+ _local_assemblers));
+ }
+ if (mesh.getDimension() > 2) {
+ Base::_secondary_variables.addSecondaryVariable(
+ "darcy_velocity_z", 1,
+ makeExtrapolator(IntegrationPointValue::DarcyVelocityZ, *_extrapolator,
+ _local_assemblers));
+ }
+}
+
+void GroundwaterFlowProcess::assembleConcreteProcess(const double t, GlobalVector const& x,
+ GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
+{
+ DBUG("Assemble GroundwaterFlowProcess.");
+
+ // Call global assembler for each local assembly item.
+ GlobalExecutor::executeMemberDereferenced(
+ *_global_assembler, &GlobalAssembler::assemble,
+ _local_assemblers, t, x, M, K, b);
+}
+
+std::unique_ptr<GroundwaterFlowProcess>
+createGroundwaterFlowProcess(
+ MeshLib::Mesh& mesh,
+ Process::NonlinearSolver& nonlinear_solver,
+ std::unique_ptr<Process::TimeDiscretization>&& time_discretization,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{process__type}
+ config.checkConfigParameter("type", "GROUNDWATER_FLOW");
+
+ DBUG("Create GroundwaterFlowProcess.");
+
+ // Process variable.
+ auto process_variables = findProcessVariables(variables, config, {
+ //! \ogs_file_param_special{process__GROUNDWATER_FLOW__process_variables__process_variable}
+ "process_variable"
+ });
+
+ // Hydraulic conductivity parameter.
+ auto& hydraulic_conductivity =
+ findParameter<double, MeshLib::Element const&>(
+ config,
+ //! \ogs_file_param_special{process__GROUNDWATER_FLOW__hydraulic_conductivity}
+ "hydraulic_conductivity",
+ parameters);
+
+ DBUG("Use \'%s\' as hydraulic conductivity parameter.",
+ hydraulic_conductivity.name.c_str());
+
+ GroundwaterFlowProcessData process_data {
+ hydraulic_conductivity
+ };
+
+ SecondaryVariableCollection secondary_variables {
+ //! \ogs_file_param{process__secondary_variables}
+ config.getConfigSubtreeOptional("secondary_variables"),
+ {
+ //! \ogs_file_param_special{process__GROUNDWATER_FLOW__secondary_variables__darcy_velocity_x}
+ "darcy_velocity_x",
+ //! \ogs_file_param_special{process__GROUNDWATER_FLOW__secondary_variables__darcy_velocity_y}
+ "darcy_velocity_y",
+ //! \ogs_file_param_special{process__GROUNDWATER_FLOW__secondary_variables__darcy_velocity_z}
+ "darcy_velocity_z"
+ }
+ };
+
+ ProcessOutput
+ //! \ogs_file_param{process__output}
+ process_output{config.getConfigSubtree("output"),
+ process_variables, secondary_variables};
+
+ return std::unique_ptr<GroundwaterFlowProcess>{
+ new GroundwaterFlowProcess{
+ mesh, nonlinear_solver,std::move(time_discretization),
+ std::move(process_variables),
+ std::move(process_data),
+ std::move(secondary_variables),
+ std::move(process_output)
+ }
+ };
+}
+
+} // namespace GroundwaterFlow
+} // namespace ProcessLib
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
index d5d2610ba9c27fee1d990ba4478231ea5aa9c9a3..279e3a2c0525e52bcc8464f0af48f15dd4ba8156 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
@@ -10,13 +10,8 @@
#ifndef PROCESS_LIB_GROUNDWATERFLOWPROCESS_H_
#define PROCESS_LIB_GROUNDWATERFLOWPROCESS_H_
-#include <cassert>
-
-#include "NumLib/Assembler/VectorMatrixAssembler.h"
#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
#include "ProcessLib/Process.h"
-#include "ProcessLib/Utils/CreateLocalAssemblers.h"
-
#include "GroundwaterFlowFEM.h"
#include "GroundwaterFlowProcessData.h"
@@ -38,27 +33,9 @@ public:
std::unique_ptr<Base::TimeDiscretization>&& time_discretization,
std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
GroundwaterFlowProcessData&& process_data,
- SecondaryVariableCollection<GlobalVector>&& secondary_variables,
- ProcessOutput<GlobalVector>&& process_output
- )
- : Process(mesh, nonlinear_solver, std::move(time_discretization),
- std::move(process_variables),
- std::move(secondary_variables),
- std::move(process_output))
- , _process_data(std::move(process_data))
- {
- if (dynamic_cast<NumLib::ForwardEuler<GlobalVector>*>(
- &Base::getTimeDiscretization()) != nullptr)
- {
- OGS_FATAL("GroundwaterFlowProcess can not be solved with the ForwardEuler"
- " time discretization scheme. Aborting");
- // Because the M matrix is not assembled. Thus, the linearized system
- // would be singular. The same applies to CrankNicolson with theta = 0.0,
- // but this case is not checked here.
- // Anyway, the GroundwaterFlowProcess shall be transferred to a simpler
- // ODESystemTag in the future.
- }
- }
+ SecondaryVariableCollection&& secondary_variables,
+ ProcessOutput&& process_output
+ );
//! \name ODESystem interface
//! @{
@@ -71,71 +48,28 @@ public:
//! @}
private:
- using LocalAssemblerInterface =
- GroundwaterFlowLocalAssemblerInterface<GlobalMatrix, GlobalVector>;
-
using GlobalAssembler = NumLib::VectorMatrixAssembler<
- GlobalMatrix, GlobalVector, LocalAssemblerInterface,
+ GroundwaterFlowLocalAssemblerInterface,
NumLib::ODESystemTag::FirstOrderImplicitQuasilinear>;
using ExtrapolatorInterface = NumLib::Extrapolator<
- GlobalVector, IntegrationPointValue, LocalAssemblerInterface>;
+ IntegrationPointValue, GroundwaterFlowLocalAssemblerInterface>;
using ExtrapolatorImplementation = NumLib::LocalLinearLeastSquaresExtrapolator<
- GlobalVector, IntegrationPointValue, LocalAssemblerInterface>;
+ IntegrationPointValue, GroundwaterFlowLocalAssemblerInterface>;
void initializeConcreteProcess(
NumLib::LocalToGlobalIndexMap const& dof_table,
MeshLib::Mesh const& mesh,
- unsigned const integration_order) override
- {
- DBUG("Create global assembler.");
- _global_assembler.reset(new GlobalAssembler(dof_table));
-
- ProcessLib::createLocalAssemblers<LocalAssemblerData>(
- mesh.getDimension(), mesh.getElements(),
- dof_table, integration_order, _local_assemblers,
- _process_data);
-
- // TOOD Later on the DOF table can change during the simulation!
- _extrapolator.reset(
- new ExtrapolatorImplementation(Base::getMatrixSpecifications(),
- *Base::_local_to_global_index_map));
-
- Base::_secondary_variables.addSecondaryVariable(
- "darcy_velocity_x", 1,
- makeExtrapolator(IntegrationPointValue::DarcyVelocityX, *_extrapolator,
- _local_assemblers));
-
- if (mesh.getDimension() > 1) {
- Base::_secondary_variables.addSecondaryVariable(
- "darcy_velocity_y", 1,
- makeExtrapolator(IntegrationPointValue::DarcyVelocityY, *_extrapolator,
- _local_assemblers));
- }
- if (mesh.getDimension() > 2) {
- Base::_secondary_variables.addSecondaryVariable(
- "darcy_velocity_z", 1,
- makeExtrapolator(IntegrationPointValue::DarcyVelocityZ, *_extrapolator,
- _local_assemblers));
- }
- }
+ unsigned const integration_order) override;
void assembleConcreteProcess(const double t, GlobalVector const& x,
- GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override
- {
- DBUG("Assemble GroundwaterFlowProcess.");
-
- // Call global assembler for each local assembly item.
- GlobalExecutor::executeMemberDereferenced(
- *_global_assembler, &GlobalAssembler::assemble,
- _local_assemblers, t, x, M, K, b);
- }
+ GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override;
GroundwaterFlowProcessData _process_data;
std::unique_ptr<GlobalAssembler> _global_assembler;
- std::vector<std::unique_ptr<LocalAssemblerInterface>> _local_assemblers;
+ std::vector<std::unique_ptr<GroundwaterFlowLocalAssemblerInterface>> _local_assemblers;
std::unique_ptr<ExtrapolatorInterface> _extrapolator;
};
@@ -147,62 +81,7 @@ createGroundwaterFlowProcess(
std::unique_ptr<Process::TimeDiscretization>&& time_discretization,
std::vector<ProcessVariable> const& variables,
std::vector<std::unique_ptr<ParameterBase>> const& parameters,
- BaseLib::ConfigTree const& config)
-{
- //! \ogs_file_param{process__type}
- config.checkConfigParameter("type", "GROUNDWATER_FLOW");
-
- DBUG("Create GroundwaterFlowProcess.");
-
- // Process variable.
- auto process_variables = findProcessVariables(variables, config, {
- //! \ogs_file_param_special{process__GROUNDWATER_FLOW__process_variables__process_variable}
- "process_variable"
- });
-
- // Hydraulic conductivity parameter.
- auto& hydraulic_conductivity =
- findParameter<double, MeshLib::Element const&>(
- config,
- //! \ogs_file_param_special{process__GROUNDWATER_FLOW__hydraulic_conductivity}
- "hydraulic_conductivity",
- parameters);
-
- DBUG("Use \'%s\' as hydraulic conductivity parameter.",
- hydraulic_conductivity.name.c_str());
-
- GroundwaterFlowProcessData process_data {
- hydraulic_conductivity
- };
-
- SecondaryVariableCollection<GlobalVector> secondary_variables {
- //! \ogs_file_param{process__secondary_variables}
- config.getConfigSubtreeOptional("secondary_variables"),
- {
- //! \ogs_file_param_special{process__GROUNDWATER_FLOW__secondary_variables__darcy_velocity_x}
- "darcy_velocity_x",
- //! \ogs_file_param_special{process__GROUNDWATER_FLOW__secondary_variables__darcy_velocity_y}
- "darcy_velocity_y",
- //! \ogs_file_param_special{process__GROUNDWATER_FLOW__secondary_variables__darcy_velocity_z}
- "darcy_velocity_z"
- }
- };
-
- ProcessOutput<GlobalVector>
- //! \ogs_file_param{process__output}
- process_output{config.getConfigSubtree("output"),
- process_variables, secondary_variables};
-
- return std::unique_ptr<GroundwaterFlowProcess>{
- new GroundwaterFlowProcess{
- mesh, nonlinear_solver,std::move(time_discretization),
- std::move(process_variables),
- std::move(process_data),
- std::move(secondary_variables),
- std::move(process_output)
- }
- };
-}
+ BaseLib::ConfigTree const& config);
} // namespace GroundwaterFlow
} // namespace ProcessLib
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index 338faa9376d3a32c70a3d5a8bd299b05496c3c37..8c43b50380a754a2b0c2c9befe8ee103ae8494b0 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -18,7 +18,6 @@ namespace ProcessLib
*
* \todo Generalize to other NumLib::ODESystemTag's.
*/
-template <typename GlobalMatrix, typename GlobalVector>
class LocalAssemblerInterface
{
public:
diff --git a/ProcessLib/NeumannBc.cpp b/ProcessLib/NeumannBc.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7b14a1695aa372e8da85369fd7a4cf86f5ff89b3
--- /dev/null
+++ b/ProcessLib/NeumannBc.cpp
@@ -0,0 +1,92 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "NeumannBc.h"
+
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+
+#include "MeshLib/MeshSearch/NodeSearch.h"
+
+#include "Utils/CreateLocalAssemblers.h"
+
+namespace ProcessLib
+{
+
+NeumannBc::NeumannBc(
+ NeumannBcConfig const& bc,
+ unsigned const integration_order,
+ NumLib::LocalToGlobalIndexMap const& local_to_global_index_map,
+ int const variable_id,
+ int const component_id)
+ : _function(*bc.getFunction()),
+ _integration_order(integration_order)
+{
+ assert(component_id < local_to_global_index_map.getNumComponents());
+
+ // deep copy because the neumann bc config destroys the elements.
+ std::transform(bc.elementsBegin(), bc.elementsEnd(),
+ std::back_inserter(_elements),
+ std::mem_fn(&MeshLib::Element::clone));
+
+ std::vector<MeshLib::Node*> nodes = MeshLib::getUniqueNodes(_elements);
+
+ auto const& mesh_subsets =
+ local_to_global_index_map.getMeshSubsets(variable_id, component_id);
+
+ // TODO extend the node intersection to all parts of mesh_subsets, i.e.
+ // to each of the MeshSubset in the mesh_subsets.
+ _mesh_subset_all_nodes =
+ mesh_subsets.getMeshSubset(0).getIntersectionByNodes(nodes);
+ std::unique_ptr<MeshLib::MeshSubsets> all_mesh_subsets{
+ new MeshLib::MeshSubsets{_mesh_subset_all_nodes}};
+
+ // Create local DOF table from intersected mesh subsets for the given
+ // variable and component ids.
+ _local_to_global_index_map.reset(
+ local_to_global_index_map.deriveBoundaryConstrainedMap(
+ variable_id, component_id, std::move(all_mesh_subsets),
+ _elements));
+}
+
+NeumannBc::~NeumannBc()
+{
+ delete _mesh_subset_all_nodes;
+
+ for (auto e : _elements)
+ delete e;
+}
+
+void NeumannBc::integrate(const double t, GlobalVector& b)
+{
+ GlobalExecutor::executeMemberDereferenced(
+ *_global_assembler, &GlobalAssembler::assemble,
+ _local_assemblers, t, b);
+}
+
+void NeumannBc::initialize(unsigned global_dim)
+{
+ DBUG("Create global assembler.");
+ _global_assembler.reset(
+ new GlobalAssembler(*_local_to_global_index_map));
+
+ auto elementValueLookup = [this](MeshLib::Element const&)
+ {
+ return _function();
+ };
+
+ createLocalAssemblers<LocalNeumannBcAsmData>(
+ global_dim, _elements,
+ *_local_to_global_index_map, _integration_order,
+ _local_assemblers,
+ elementValueLookup
+ );
+}
+
+} // namespace ProcessLib
diff --git a/ProcessLib/NeumannBc.h b/ProcessLib/NeumannBc.h
index c666ec4787b65e53d09b3a21da5afcf56290ff5f..0808f0e1f3ec2df13ccf7f5693ac313fac3aa71f 100644
--- a/ProcessLib/NeumannBc.h
+++ b/ProcessLib/NeumannBc.h
@@ -13,15 +13,9 @@
#include <memory>
#include <vector>
-#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
-#include "NumLib/Fem/ShapeMatrixPolicy.h"
-
-#include "NumLib/Assembler/VectorMatrixAssembler.h"
#include "MeshLib/MeshSubset.h"
-#include "MeshLib/MeshSearch/NodeSearch.h"
-
-#include "Utils/CreateLocalAssemblers.h"
-
+#include "NumLib/NumericsConfig.h"
+#include "NumLib/Assembler/VectorMatrixAssembler.h"
#include "NeumannBcConfig.h"
#include "NeumannBcAssembler.h"
@@ -52,72 +46,15 @@ public:
unsigned const integration_order,
NumLib::LocalToGlobalIndexMap const& local_to_global_index_map,
int const variable_id,
- int const component_id)
- : _function(*bc.getFunction()),
- _integration_order(integration_order)
- {
- assert(component_id < local_to_global_index_map.getNumComponents());
-
- // deep copy because the neumann bc config destroys the elements.
- std::transform(bc.elementsBegin(), bc.elementsEnd(),
- std::back_inserter(_elements),
- std::mem_fn(&MeshLib::Element::clone));
-
- std::vector<MeshLib::Node*> nodes = MeshLib::getUniqueNodes(_elements);
-
- auto const& mesh_subsets =
- local_to_global_index_map.getMeshSubsets(variable_id, component_id);
-
- // TODO extend the node intersection to all parts of mesh_subsets, i.e.
- // to each of the MeshSubset in the mesh_subsets.
- _mesh_subset_all_nodes =
- mesh_subsets.getMeshSubset(0).getIntersectionByNodes(nodes);
- std::unique_ptr<MeshLib::MeshSubsets> all_mesh_subsets{
- new MeshLib::MeshSubsets{_mesh_subset_all_nodes}};
-
- // Create local DOF table from intersected mesh subsets for the given
- // variable and component ids.
- _local_to_global_index_map.reset(
- local_to_global_index_map.deriveBoundaryConstrainedMap(
- variable_id, component_id, std::move(all_mesh_subsets),
- _elements));
- }
-
- ~NeumannBc()
- {
- delete _mesh_subset_all_nodes;
-
- for (auto e : _elements)
- delete e;
- }
+ int const component_id);
+
+ ~NeumannBc();
/// Calls local assemblers which calculate their contributions to the global
/// matrix and the right-hand-side.
- void integrate(const double t, GlobalVector& b)
- {
- GlobalExecutor::executeMemberDereferenced(
- *_global_assembler, &GlobalAssembler::assemble,
- _local_assemblers, t, b);
- }
-
- void initialize(unsigned global_dim)
- {
- DBUG("Create global assembler.");
- _global_assembler.reset(
- new GlobalAssembler(*_local_to_global_index_map));
-
- auto elementValueLookup = [this](MeshLib::Element const&)
- {
- return _function();
- };
-
- createLocalAssemblers<LocalNeumannBcAsmData>(
- global_dim, _elements,
- *_local_to_global_index_map, _integration_order,
- _local_assemblers,
- elementValueLookup
- );
- }
+ void integrate(const double t, GlobalVector& b);
+
+ void initialize(unsigned global_dim);
private:
/// The right-hand-side function of the Neumann boundary condition given as
@@ -138,18 +75,15 @@ private:
/// the #_function.
unsigned const _integration_order;
- using LocalAssembler = LocalNeumannBcAsmDataInterface<
- GlobalMatrix, GlobalVector>;
-
using GlobalAssembler =
NumLib::VectorMatrixAssembler<
- GlobalMatrix, GlobalVector, LocalAssembler,
+ LocalNeumannBcAsmDataInterface,
NumLib::ODESystemTag::NeumannBC>;
std::unique_ptr<GlobalAssembler> _global_assembler;
/// Local assemblers for each element of #_elements.
- std::vector<std::unique_ptr<LocalAssembler>> _local_assemblers;
+ std::vector<std::unique_ptr<LocalNeumannBcAsmDataInterface>> _local_assemblers;
};
diff --git a/ProcessLib/NeumannBcAssembler.h b/ProcessLib/NeumannBcAssembler.h
index 9da05d9453cb0553ee3b87651ef6a405b20fb035..1e13a3996f075c5180bd3052f2637954620006f3 100644
--- a/ProcessLib/NeumannBcAssembler.h
+++ b/ProcessLib/NeumannBcAssembler.h
@@ -19,7 +19,6 @@
namespace ProcessLib
{
-template <typename GlobalMatrix, typename GlobalVector>
class LocalNeumannBcAsmDataInterface
{
public:
@@ -33,10 +32,8 @@ public:
template <typename ShapeFunction_,
typename IntegrationMethod_,
- typename GlobalMatrix,
- typename GlobalVector,
unsigned GlobalDim>
-class LocalNeumannBcAsmData : public LocalNeumannBcAsmDataInterface<GlobalMatrix, GlobalVector>
+class LocalNeumannBcAsmData : public LocalNeumannBcAsmDataInterface
{
public:
using ShapeFunction = ShapeFunction_;
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index d2d245c06a8b425321caf57f4675fe650b238901..46fdb7d66324a0a04a80af3ae0e8975a047a9a00 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -9,8 +9,242 @@
#include "Process.h"
+#include "NumLib/DOF/ComputeSparsityPattern.h"
+#include "DirichletBc.h"
+#include "NeumannBc.h"
+#include "NeumannBcAssembler.h"
+#include "ProcessVariable.h"
+#include "UniformDirichletBoundaryCondition.h"
+
namespace ProcessLib
{
+
+Process::Process(
+ MeshLib::Mesh& mesh,
+ NonlinearSolver& nonlinear_solver,
+ std::unique_ptr<TimeDiscretization>&& time_discretization,
+ std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
+ SecondaryVariableCollection&& secondary_variables,
+ ProcessOutput&& process_output
+ )
+ : _mesh(mesh)
+ , _secondary_variables(std::move(secondary_variables))
+ , _process_output(std::move(process_output))
+ , _nonlinear_solver(nonlinear_solver)
+ , _time_discretization(std::move(time_discretization))
+ , _process_variables(std::move(process_variables))
+{}
+
+void Process::output(std::string const& file_name,
+ const unsigned /*timestep*/,
+ GlobalVector const& x) const
+{
+ doProcessOutput(file_name, x, _mesh, *_local_to_global_index_map,
+ _process_variables, _secondary_variables, _process_output);
+}
+
+void Process::initialize()
+{
+ DBUG("Initialize process.");
+
+ DBUG("Construct dof mappings.");
+ constructDofTable();
+
+ DBUG("Compute sparsity pattern");
+ computeSparsityPattern();
+
+ initializeConcreteProcess(*_local_to_global_index_map, _mesh,
+ _integration_order);
+
+ DBUG("Initialize boundary conditions.");
+ for (int variable_id = 0;
+ variable_id < static_cast<int>(_process_variables.size());
+ ++variable_id)
+ {
+ ProcessVariable& pv = _process_variables[variable_id];
+ for (int component_id = 0;
+ component_id < pv.getNumberOfComponents();
+ ++component_id)
+ {
+ createDirichletBcs(pv, variable_id, component_id);
+ createNeumannBcs(pv, variable_id, component_id);
+ }
+ }
+
+ for (auto& bc : _neumann_bcs)
+ bc->initialize(_mesh.getDimension());
+}
+
+void Process::setInitialConditions(GlobalVector& x)
+{
+ DBUG("Set initial conditions.");
+ for (int variable_id = 0;
+ variable_id < static_cast<int>(_process_variables.size());
+ ++variable_id)
+ {
+ ProcessVariable& pv = _process_variables[variable_id];
+ for (int component_id = 0;
+ component_id < pv.getNumberOfComponents();
+ ++component_id)
+ {
+ setInitialConditions(pv, variable_id, component_id, x);
+ }
+ }
+}
+
+MathLib::MatrixSpecifications Process::getMatrixSpecifications() const
+{
+ auto const& l = *_local_to_global_index_map;
+ return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
+ &l.getGhostIndices(), &_sparsity_pattern};
+}
+
+void Process::assemble(const double t, GlobalVector const& x,
+ GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
+{
+ assembleConcreteProcess(t, x, M, K, b);
+
+ // Call global assembler for each Neumann boundary local assembler.
+ for (auto const& bc : _neumann_bcs)
+ bc->integrate(t, 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)
+{
+ 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.
+}
+
+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.");
+}
+
+void Process::constructDofTable()
+{
+ // Create single component dof in every of the mesh's nodes.
+ _mesh_subset_all_nodes.reset(
+ new MeshLib::MeshSubset(_mesh, &_mesh.getNodes()));
+
+ // Collect the mesh subsets in a vector.
+ std::vector<std::unique_ptr<MeshLib::MeshSubsets>> all_mesh_subsets;
+ for (ProcessVariable const& pv : _process_variables)
+ {
+ std::generate_n(
+ std::back_inserter(all_mesh_subsets),
+ pv.getNumberOfComponents(),
+ [&]()
+ {
+ return std::unique_ptr<MeshLib::MeshSubsets>{
+ new MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()}};
+ });
+ }
+
+ _local_to_global_index_map.reset(
+ new NumLib::LocalToGlobalIndexMap(
+ std::move(all_mesh_subsets),
+ NumLib::ComponentOrder::BY_LOCATION));
+}
+
+void Process::setInitialConditions(ProcessVariable const& variable,
+ int const variable_id,
+ int const component_id,
+ GlobalVector& x)
+{
+ std::size_t const n_nodes = _mesh.getNumberOfNodes();
+ for (std::size_t node_id = 0; node_id < n_nodes; ++node_id)
+ {
+ MeshLib::Location const l(_mesh.getID(),
+ MeshLib::MeshItemType::Node, node_id);
+ auto global_index =
+ std::abs(_local_to_global_index_map->getGlobalIndex(
+ l, variable_id, component_id));
+#ifdef USE_PETSC
+ // The global indices of the ghost entries of the global matrix or
+ // the global vectors need to be set as negative values for equation
+ // assembly, however the global indices start from zero. Therefore,
+ // any ghost entry with zero index is assigned an negative value of
+ // the vector size or the matrix dimension. To assign the initial
+ // value for the ghost entries, the negative indices of the ghost
+ // entries are restored to zero.
+ if (global_index == x.size())
+ global_index = 0;
+#endif
+ x.set(global_index,
+ variable.getInitialConditionValue(node_id, component_id));
+ }
+}
+
+void Process::createDirichletBcs(ProcessVariable& variable, int const variable_id,
+ int const component_id)
+{
+ MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher =
+ MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher(
+ variable.getMesh());
+
+ variable.initializeDirichletBCs(std::back_inserter(_dirichlet_bcs),
+ mesh_node_searcher,
+ *_local_to_global_index_map,
+ variable_id,
+ component_id);
+}
+
+void Process::createNeumannBcs(ProcessVariable& variable, int const variable_id,
+ int const component_id)
+{
+ // Find mesh nodes.
+ MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher =
+ MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher(
+ variable.getMesh());
+ MeshGeoToolsLib::BoundaryElementsSearcher mesh_element_searcher(
+ variable.getMesh(), mesh_node_searcher);
+
+ // Create a neumann BC for the process variable storing them in the
+ // _neumann_bcs vector.
+ variable.createNeumannBcs(
+ std::back_inserter(_neumann_bcs),
+ mesh_element_searcher,
+ _integration_order,
+ *_local_to_global_index_map,
+ variable_id,
+ component_id);
+}
+
+void Process::computeSparsityPattern()
+{
+ _sparsity_pattern = NumLib::computeSparsityPattern(
+ *_local_to_global_index_map, _mesh);
+}
+
ProcessVariable& findProcessVariable(
std::vector<ProcessVariable> const& variables,
BaseLib::ConfigTree const& pv_config, std::string const& tag)
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index a8d8e065f2dd3681bea130919c2da134cb4e2792..f7fe21fc50fccf897a1e6ea42882ff1ede992b52 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -10,19 +10,13 @@
#ifndef PROCESS_LIB_PROCESS_H_
#define PROCESS_LIB_PROCESS_H_
-#include "MeshLib/IO/VtkIO/VtuInterface.h"
-#include "NumLib/DOF/ComputeSparsityPattern.h"
#include "NumLib/ODESolver/ODESystem.h"
#include "NumLib/ODESolver/TimeDiscretization.h"
#include "NumLib/ODESolver/NonlinearSolver.h"
-#include "DirichletBc.h"
-#include "NeumannBc.h"
-#include "NeumannBcAssembler.h"
#include "Parameter.h"
#include "ProcessOutput.h"
-#include "ProcessVariable.h"
-#include "UniformDirichletBoundaryCondition.h"
+#include "SecondaryVariable.h"
namespace MeshLib
{
@@ -49,16 +43,9 @@ public:
NonlinearSolver& nonlinear_solver,
std::unique_ptr<TimeDiscretization>&& time_discretization,
std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
- SecondaryVariableCollection<GlobalVector>&& secondary_variables,
- ProcessOutput<GlobalVector>&& process_output
- )
- : _mesh(mesh)
- , _secondary_variables(std::move(secondary_variables))
- , _process_output(std::move(process_output))
- , _nonlinear_solver(nonlinear_solver)
- , _time_discretization(std::move(time_discretization))
- , _process_variables(std::move(process_variables))
- {}
+ SecondaryVariableCollection&& secondary_variables,
+ ProcessOutput&& process_output
+ );
/// Preprocessing before starting assembly for new timestep.
virtual void preTimestep(GlobalVector const& /*x*/,
@@ -71,105 +58,22 @@ public:
/// The file_name is indicating the name of possible output file.
void output(std::string const& file_name,
const unsigned /*timestep*/,
- GlobalVector const& x) const
- {
- doProcessOutput(file_name, x, _mesh, *_local_to_global_index_map,
- _process_variables, _secondary_variables, _process_output);
- }
+ GlobalVector const& x) const;
- void initialize()
- {
- DBUG("Initialize process.");
-
- DBUG("Construct dof mappings.");
- constructDofTable();
-
- DBUG("Compute sparsity pattern");
- computeSparsityPattern();
-
- initializeConcreteProcess(*_local_to_global_index_map, _mesh,
- _integration_order);
-
- DBUG("Initialize boundary conditions.");
- for (int variable_id = 0;
- variable_id < static_cast<int>(_process_variables.size());
- ++variable_id)
- {
- ProcessVariable& pv = _process_variables[variable_id];
- for (int component_id = 0;
- component_id < pv.getNumberOfComponents();
- ++component_id)
- {
- createDirichletBcs(pv, variable_id, component_id);
- createNeumannBcs(pv, variable_id, component_id);
- }
- }
-
- for (auto& bc : _neumann_bcs)
- bc->initialize(_mesh.getDimension());
- }
+ void initialize();
- void setInitialConditions(GlobalVector& x)
- {
- DBUG("Set initial conditions.");
- for (int variable_id = 0;
- variable_id < static_cast<int>(_process_variables.size());
- ++variable_id)
- {
- ProcessVariable& pv = _process_variables[variable_id];
- for (int component_id = 0;
- component_id < pv.getNumberOfComponents();
- ++component_id)
- {
- setInitialConditions(pv, variable_id, component_id, x);
- }
- }
- }
+ void setInitialConditions(GlobalVector& x);
- MathLib::MatrixSpecifications getMatrixSpecifications() const override final
- {
- auto const& l = *_local_to_global_index_map;
- return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
- &l.getGhostIndices(), &_sparsity_pattern};
- }
+ MathLib::MatrixSpecifications getMatrixSpecifications() const override final;
void assemble(const double t, GlobalVector const& x,
- GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override final
- {
- assembleConcreteProcess(t, x, M, K, b);
-
- // Call global assembler for each Neumann boundary local assembler.
- for (auto const& bc : _neumann_bcs)
- bc->integrate(t, b);
- }
+ 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
- {
- 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.
- }
+ GlobalMatrix& Jac) override final;
std::vector<DirichletBc<Index>> const* getKnownSolutions(
double const /*t*/) const override final
@@ -202,113 +106,26 @@ private:
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.");
- }
+ GlobalMatrix& /*Jac*/);
- void constructDofTable()
- {
- // Create single component dof in every of the mesh's nodes.
- _mesh_subset_all_nodes.reset(
- new MeshLib::MeshSubset(_mesh, &_mesh.getNodes()));
-
- // Collect the mesh subsets in a vector.
- std::vector<std::unique_ptr<MeshLib::MeshSubsets>> all_mesh_subsets;
- for (ProcessVariable const& pv : _process_variables)
- {
- std::generate_n(
- std::back_inserter(all_mesh_subsets),
- pv.getNumberOfComponents(),
- [&]()
- {
- return std::unique_ptr<MeshLib::MeshSubsets>{
- new MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()}};
- });
- }
-
- _local_to_global_index_map.reset(
- new NumLib::LocalToGlobalIndexMap(
- std::move(all_mesh_subsets),
- NumLib::ComponentOrder::BY_LOCATION));
- }
+ void constructDofTable();
/// Sets the initial condition values in the solution vector x for a given
/// process variable and component.
void setInitialConditions(ProcessVariable const& variable,
int const variable_id,
int const component_id,
- GlobalVector& x)
- {
- std::size_t const n_nodes = _mesh.getNumberOfNodes();
- for (std::size_t node_id = 0; node_id < n_nodes; ++node_id)
- {
- MeshLib::Location const l(_mesh.getID(),
- MeshLib::MeshItemType::Node, node_id);
- auto global_index =
- std::abs(_local_to_global_index_map->getGlobalIndex(
- l, variable_id, component_id));
-#ifdef USE_PETSC
- // The global indices of the ghost entries of the global matrix or
- // the global vectors need to be set as negative values for equation
- // assembly, however the global indices start from zero. Therefore,
- // any ghost entry with zero index is assigned an negative value of
- // the vector size or the matrix dimension. To assign the initial
- // value for the ghost entries, the negative indices of the ghost
- // entries are restored to zero.
- if (global_index == x.size())
- global_index = 0;
-#endif
- x.set(global_index,
- variable.getInitialConditionValue(node_id, component_id));
- }
- }
+ GlobalVector& x);
void createDirichletBcs(ProcessVariable& variable, int const variable_id,
- int const component_id)
- {
- MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher =
- MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher(
- variable.getMesh());
-
- variable.initializeDirichletBCs(std::back_inserter(_dirichlet_bcs),
- mesh_node_searcher,
- *_local_to_global_index_map,
- variable_id,
- component_id);
- }
+ int const component_id);
void createNeumannBcs(ProcessVariable& variable, int const variable_id,
- int const component_id)
- {
- // Find mesh nodes.
- MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher =
- MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher(
- variable.getMesh());
- MeshGeoToolsLib::BoundaryElementsSearcher mesh_element_searcher(
- variable.getMesh(), mesh_node_searcher);
-
- // Create a neumann BC for the process variable storing them in the
- // _neumann_bcs vector.
- variable.createNeumannBcs(
- std::back_inserter(_neumann_bcs),
- mesh_element_searcher,
- _integration_order,
- *_local_to_global_index_map,
- variable_id,
- component_id);
- }
+ int const component_id);
/// Computes and stores global matrix' sparsity pattern from given
/// DOF-table.
- void computeSparsityPattern()
- {
- _sparsity_pattern =
- NumLib::computeSparsityPattern(*_local_to_global_index_map, _mesh);
- }
+ void computeSparsityPattern();
protected:
MeshLib::Mesh& _mesh;
@@ -317,8 +134,8 @@ protected:
std::unique_ptr<NumLib::LocalToGlobalIndexMap>
_local_to_global_index_map;
- SecondaryVariableCollection<GlobalVector> _secondary_variables;
- ProcessOutput<GlobalVector> _process_output;
+ SecondaryVariableCollection _secondary_variables;
+ ProcessOutput _process_output;
private:
unsigned const _integration_order = 2;
diff --git a/ProcessLib/ProcessOutput.cpp b/ProcessLib/ProcessOutput.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6579726fba3f6eac79cafcd1852b27950a23fdc9
--- /dev/null
+++ b/ProcessLib/ProcessOutput.cpp
@@ -0,0 +1,236 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "ProcessOutput.h"
+
+#include "MeshLib/IO/VtkIO/VtuInterface.h"
+
+namespace ProcessLib
+{
+
+ProcessOutput::ProcessOutput(BaseLib::ConfigTree const& output_config,
+ std::vector<std::reference_wrapper<ProcessVariable>> const&
+ process_variables,
+ SecondaryVariableCollection const& secondary_variables)
+{
+ //! \ogs_file_param{process__output__variables}
+ auto const out_vars = output_config.getConfigSubtree("variables");
+
+ //! \ogs_file_param{process__output__variables__variable}
+ for (auto out_var : out_vars.getConfigParameterList<std::string>("variable"))
+ {
+ if (output_variables.find(out_var) != output_variables.cend())
+ {
+ OGS_FATAL("output variable `%s' specified more than once.", out_var.c_str());
+ }
+
+ auto pred = [&out_var](ProcessVariable const& pv) {
+ return pv.getName() == out_var;
+ };
+
+ // check if out_var is a process variable
+ auto const& pcs_var = std::find_if(
+ process_variables.cbegin(), process_variables.cend(), pred);
+
+ if (pcs_var == process_variables.cend()
+ && !secondary_variables.variableExists(out_var))
+ {
+ OGS_FATAL("Output variable `%s' is neither a process variable nor a"
+ " secondary variable", out_var.c_str());
+ }
+
+ DBUG("adding output variable `%s'", out_var.c_str());
+ output_variables.insert(out_var);
+ }
+
+ if (auto out_resid = output_config.getConfigParameterOptional<bool>(
+ "output_extrapolation_residuals")) {
+ output_residuals = *out_resid;
+ }
+
+ // debug output
+ if (auto const param =
+ //! \ogs_file_param{process__output__output_iteration_results}
+ output_config.getConfigParameterOptional<bool>("output_iteration_results"))
+ {
+ DBUG("output_iteration_results: %s", (*param) ? "true" : "false");
+
+ output_iteration_results = *param;
+ }
+}
+
+
+void doProcessOutput(
+ std::string const& file_name,
+ GlobalVector const& x,
+ MeshLib::Mesh& mesh,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<std::reference_wrapper<ProcessVariable>> const&
+ process_variables,
+ SecondaryVariableCollection secondary_variables,
+ ProcessOutput const& process_output)
+{
+ DBUG("Process output.");
+
+ // Copy result
+#ifdef USE_PETSC
+ // TODO It is also possible directly to copy the data for single process
+ // variable to a mesh property. It needs a vector of global indices and
+ // some PETSc magic to do so.
+ std::vector<double> x_copy(x.getLocalSize() + x.getGhostSize());
+#else
+ std::vector<double> x_copy(x.size());
+#endif
+ x.copyValues(x_copy);
+
+ auto const& output_variables = process_output.output_variables;
+
+ std::size_t const n_nodes = mesh.getNumberOfNodes();
+ int global_component_offset = 0;
+ int global_component_offset_next = 0;
+
+ // primary variables
+ for (ProcessVariable& pv : process_variables)
+ {
+ int const n_components = pv.getNumberOfComponents();
+ global_component_offset = global_component_offset_next;
+ global_component_offset_next += n_components;
+
+ if (output_variables.find(pv.getName()) == output_variables.cend())
+ continue;
+
+ DBUG(" process variable %s", pv.getName().c_str());
+
+ auto& output_data = pv.getOrCreateMeshProperty();
+
+ for (std::size_t node_id = 0; node_id < n_nodes; ++node_id)
+ {
+ MeshLib::Location const l(mesh.getID(),
+ MeshLib::MeshItemType::Node, node_id);
+ // TODO extend component ids to multiple process variables.
+ for (int component_id = 0; component_id < n_components;
+ ++component_id)
+ {
+ auto const global_component_id = global_component_offset + component_id;
+ auto const index =
+ dof_table.getLocalIndex(
+ l, global_component_id, x.getRangeBegin(),
+ x.getRangeEnd());
+
+ output_data[node_id * n_components + component_id] =
+ x_copy[index];
+ }
+ }
+ }
+
+#ifndef USE_PETSC
+ // the following section is for the output of secondary variables
+
+ auto count_mesh_items = [](
+ MeshLib::Mesh const& mesh, MeshLib::MeshItemType type) -> std::size_t
+ {
+ switch (type) {
+ case MeshLib::MeshItemType::Cell: return mesh.getNumberOfElements();
+ case MeshLib::MeshItemType::Node: return mesh.getNumberOfNodes();
+ default: break; // avoid compiler warning
+ }
+ return 0;
+ };
+
+ auto get_or_create_mesh_property = [&mesh, &count_mesh_items](
+ std::string const& property_name, MeshLib::MeshItemType type)
+ {
+ // Get or create a property vector for results.
+ boost::optional<MeshLib::PropertyVector<double>&> result;
+
+ auto const N = count_mesh_items(mesh, type);
+
+ if (mesh.getProperties().hasPropertyVector(property_name))
+ {
+ result = mesh.getProperties().template
+ getPropertyVector<double>(property_name);
+ }
+ else
+ {
+ result = mesh.getProperties().template
+ createNewPropertyVector<double>(property_name, type);
+ result->resize(N);
+ }
+ assert(result && result->size() == N);
+
+ return result;
+ };
+
+ auto add_secondary_var = [&](SecondaryVariable const& var,
+ std::string const& output_name)
+ {
+ assert(var.n_components == 1); // TODO implement other cases
+
+ {
+ DBUG(" secondary variable %s", output_name.c_str());
+
+ auto result = get_or_create_mesh_property(
+ output_name, MeshLib::MeshItemType::Node);
+ assert(result->size() == mesh.getNumberOfNodes());
+
+ std::unique_ptr<GlobalVector> result_cache;
+ auto const& nodal_values =
+ var.fcts.eval_field(x, dof_table, result_cache);
+
+ // Copy result
+ for (std::size_t i = 0; i < mesh.getNumberOfNodes(); ++i)
+ {
+ assert(!std::isnan(nodal_values[i]));
+ (*result)[i] = nodal_values[i];
+ }
+ }
+
+ if (process_output.output_residuals && var.fcts.eval_residuals)
+ {
+ DBUG(" secondary variable %s residual", output_name.c_str());
+ auto const& property_name_res = output_name + "_residual";
+
+ auto result = get_or_create_mesh_property(
+ property_name_res, MeshLib::MeshItemType::Cell);
+ assert(result->size() == mesh.getNumberOfElements());
+
+ std::unique_ptr<GlobalVector> result_cache;
+ auto const& residuals =
+ var.fcts.eval_residuals(x, dof_table, result_cache);
+
+ // Copy result
+ for (std::size_t i = 0; i < mesh.getNumberOfElements(); ++i)
+ {
+ assert(!std::isnan(residuals[i]));
+ (*result)[i] = residuals[i];
+ }
+ }
+ };
+
+ for (auto const& elem : secondary_variables)
+ {
+ auto const& var_name = elem.first;
+ if (output_variables.find(var_name) != output_variables.cend())
+ {
+ add_secondary_var(elem.second, var_name);
+ }
+ }
+
+ // secondary variables output end
+#else
+ (void) secondary_variables;
+#endif // USE_PETSC
+
+ // Write output file
+ DBUG("Writing output to \'%s\'.", file_name.c_str());
+ MeshLib::IO::VtuInterface vtu_interface(&mesh, vtkXMLWriter::Binary, true);
+ vtu_interface.writeToFile(file_name);
+}
+
+} // ProcessLib
diff --git a/ProcessLib/ProcessOutput.h b/ProcessLib/ProcessOutput.h
index 96017e01b722c35149c564d9d52d6d87787fa5e6..d791b0ac17c5f561ceceb0d1b80da9e2835ac8cf 100644
--- a/ProcessLib/ProcessOutput.h
+++ b/ProcessLib/ProcessOutput.h
@@ -10,7 +10,6 @@
#ifndef PROCESSLIB_PROCESSOUTPUT_H
#define PROCESSLIB_PROCESSOUTPUT_H
-#include "MeshLib/IO/VtkIO/VtuInterface.h"
#include "ProcessVariable.h"
#include "SecondaryVariable.h"
@@ -18,60 +17,13 @@ namespace ProcessLib
{
//! Holds information about which variables to write to output files.
-template <typename GlobalVector>
struct ProcessOutput final
{
//! Constructs a new instance.
ProcessOutput(BaseLib::ConfigTree const& output_config,
std::vector<std::reference_wrapper<ProcessVariable>> const&
process_variables,
- SecondaryVariableCollection<GlobalVector> const& secondary_variables)
- {
- //! \ogs_file_param{process__output__variables}
- auto const out_vars = output_config.getConfigSubtree("variables");
-
- //! \ogs_file_param{process__output__variables__variable}
- for (auto out_var : out_vars.getConfigParameterList<std::string>("variable"))
- {
- if (output_variables.find(out_var) != output_variables.cend())
- {
- OGS_FATAL("output variable `%s' specified more than once.", out_var.c_str());
- }
-
- auto pred = [&out_var](ProcessVariable const& pv) {
- return pv.getName() == out_var;
- };
-
- // check if out_var is a process variable
- auto const& pcs_var = std::find_if(
- process_variables.cbegin(), process_variables.cend(), pred);
-
- if (pcs_var == process_variables.cend()
- && !secondary_variables.variableExists(out_var))
- {
- OGS_FATAL("Output variable `%s' is neither a process variable nor a"
- " secondary variable", out_var.c_str());
- }
-
- DBUG("adding output variable `%s'", out_var.c_str());
- output_variables.insert(out_var);
- }
-
- if (auto out_resid = output_config.getConfigParameterOptional<bool>(
- "output_extrapolation_residuals")) {
- output_residuals = *out_resid;
- }
-
- // debug output
- if (auto const param =
- //! \ogs_file_param{process__output__output_iteration_results}
- output_config.getConfigParameterOptional<bool>("output_iteration_results"))
- {
- DBUG("output_iteration_results: %s", (*param) ? "true" : "false");
-
- output_iteration_results = *param;
- }
- }
+ SecondaryVariableCollection const& secondary_variables);
//! All variables that shall be output.
std::set<std::string> output_variables;
@@ -84,7 +36,6 @@ struct ProcessOutput final
//! Writes output to the given \c file_name using the VTU file format.
-template <typename GlobalVector>
void doProcessOutput(
std::string const& file_name,
GlobalVector const& x,
@@ -92,165 +43,8 @@ void doProcessOutput(
NumLib::LocalToGlobalIndexMap const& dof_table,
std::vector<std::reference_wrapper<ProcessVariable>> const&
process_variables,
- SecondaryVariableCollection<GlobalVector> secondary_variables,
- ProcessOutput<GlobalVector> const& process_output)
-{
- DBUG("Process output.");
-
- // Copy result
-#ifdef USE_PETSC
- // TODO It is also possible directly to copy the data for single process
- // variable to a mesh property. It needs a vector of global indices and
- // some PETSc magic to do so.
- std::vector<double> x_copy(x.getLocalSize() + x.getGhostSize());
-#else
- std::vector<double> x_copy(x.size());
-#endif
- x.copyValues(x_copy);
-
- auto const& output_variables = process_output.output_variables;
-
- std::size_t const n_nodes = mesh.getNumberOfNodes();
- int global_component_offset = 0;
- int global_component_offset_next = 0;
-
- // primary variables
- for (ProcessVariable& pv : process_variables)
- {
- int const n_components = pv.getNumberOfComponents();
- global_component_offset = global_component_offset_next;
- global_component_offset_next += n_components;
-
- if (output_variables.find(pv.getName()) == output_variables.cend())
- continue;
-
- DBUG(" process variable %s", pv.getName().c_str());
-
- auto& output_data = pv.getOrCreateMeshProperty();
-
- for (std::size_t node_id = 0; node_id < n_nodes; ++node_id)
- {
- MeshLib::Location const l(mesh.getID(),
- MeshLib::MeshItemType::Node, node_id);
- // TODO extend component ids to multiple process variables.
- for (int component_id = 0; component_id < n_components;
- ++component_id)
- {
- auto const global_component_id = global_component_offset + component_id;
- auto const index =
- dof_table.getLocalIndex(
- l, global_component_id, x.getRangeBegin(),
- x.getRangeEnd());
-
- output_data[node_id * n_components + component_id] =
- x_copy[index];
- }
- }
- }
-
-#ifndef USE_PETSC
- // the following section is for the output of secondary variables
-
- auto count_mesh_items = [](
- MeshLib::Mesh const& mesh, MeshLib::MeshItemType type) -> std::size_t
- {
- switch (type) {
- case MeshLib::MeshItemType::Cell: return mesh.getNumberOfElements();
- case MeshLib::MeshItemType::Node: return mesh.getNumberOfNodes();
- default: break; // avoid compiler warning
- }
- return 0;
- };
-
- auto get_or_create_mesh_property = [&mesh, &count_mesh_items](
- std::string const& property_name, MeshLib::MeshItemType type)
- {
- // Get or create a property vector for results.
- boost::optional<MeshLib::PropertyVector<double>&> result;
-
- auto const N = count_mesh_items(mesh, type);
-
- if (mesh.getProperties().hasPropertyVector(property_name))
- {
- result = mesh.getProperties().template
- getPropertyVector<double>(property_name);
- }
- else
- {
- result = mesh.getProperties().template
- createNewPropertyVector<double>(property_name, type);
- result->resize(N);
- }
- assert(result && result->size() == N);
-
- return result;
- };
-
- auto add_secondary_var = [&](SecondaryVariable<GlobalVector> const& var,
- std::string const& output_name)
- {
- assert(var.n_components == 1); // TODO implement other cases
-
- {
- DBUG(" secondary variable %s", output_name.c_str());
-
- auto result = get_or_create_mesh_property(
- output_name, MeshLib::MeshItemType::Node);
- assert(result->size() == mesh.getNumberOfNodes());
-
- std::unique_ptr<GlobalVector> result_cache;
- auto const& nodal_values =
- var.fcts.eval_field(x, dof_table, result_cache);
-
- // Copy result
- for (std::size_t i = 0; i < mesh.getNumberOfNodes(); ++i)
- {
- assert(!std::isnan(nodal_values[i]));
- (*result)[i] = nodal_values[i];
- }
- }
-
- if (process_output.output_residuals && var.fcts.eval_residuals)
- {
- DBUG(" secondary variable %s residual", output_name.c_str());
- auto const& property_name_res = output_name + "_residual";
-
- auto result = get_or_create_mesh_property(
- property_name_res, MeshLib::MeshItemType::Cell);
- assert(result->size() == mesh.getNumberOfElements());
-
- std::unique_ptr<GlobalVector> result_cache;
- auto const& residuals =
- var.fcts.eval_residuals(x, dof_table, result_cache);
-
- // Copy result
- for (std::size_t i = 0; i < mesh.getNumberOfElements(); ++i)
- {
- assert(!std::isnan(residuals[i]));
- (*result)[i] = residuals[i];
- }
- }
- };
-
- for (auto const& elem : secondary_variables)
- {
- auto const& var_name = elem.first;
- if (output_variables.find(var_name) != output_variables.cend())
- {
- add_secondary_var(elem.second, var_name);
- }
- }
-
- // secondary variables output end
-#else
- (void) secondary_variables;
-#endif // USE_PETSC
-
- // Write output file
- DBUG("Writing output to \'%s\'.", file_name.c_str());
- MeshLib::IO::VtuInterface vtu_interface(&mesh, vtkXMLWriter::Binary, true);
- vtu_interface.writeToFile(file_name);
-}
+ SecondaryVariableCollection secondary_variables,
+ ProcessOutput const& process_output);
} // ProcessLib
diff --git a/ProcessLib/SecondaryVariable.cpp b/ProcessLib/SecondaryVariable.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3ceff2c9e083721bc4f45eb4a17585a772aa96d9
--- /dev/null
+++ b/ProcessLib/SecondaryVariable.cpp
@@ -0,0 +1,84 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "SecondaryVariable.h"
+
+namespace ProcessLib
+{
+
+SecondaryVariableCollection::SecondaryVariableCollection(
+ boost::optional<BaseLib::ConfigTree> const& config,
+ std::initializer_list<std::string> tag_names)
+{
+ if (!config) return;
+
+ // read which variables are defined in the config
+ for (auto const& tag_name : tag_names) {
+ if (!_all_secondary_variables.insert(tag_name).second) {
+ OGS_FATAL("Tag name <%s> has been specified twice as a secondary variable.");
+ }
+
+ //! \ogs_file_special
+ if (auto var_name = config->getConfigParameterOptional<std::string>(tag_name))
+ {
+ // TODO check primary vars, too
+ BaseLib::insertIfKeyValueUniqueElseError(
+ _map_tagname_to_varname, tag_name, *var_name,
+ "Secondary variable names must be unique.");
+ }
+ }
+}
+
+bool SecondaryVariableCollection::variableExists(std::string const& variable_name) const
+{
+ auto pred = [&variable_name](std::pair<std::string, std::string> const& p) {
+ return p.second == variable_name;
+ };
+
+ // check if out_var is a secondary variable
+ auto const& var = std::find_if(
+ _map_tagname_to_varname.cbegin(), _map_tagname_to_varname.cend(), pred);
+
+ return var != _map_tagname_to_varname.cend();
+}
+
+void SecondaryVariableCollection::addSecondaryVariable(std::string const& tag_name,
+ const unsigned num_components,
+ SecondaryVariableFunctions&& fcts)
+{
+ auto it = _map_tagname_to_varname.find(tag_name);
+
+ // get user-supplied var_name for the given tag_name
+ if (it != _map_tagname_to_varname.end())
+ {
+ auto const& var_name = it->first;
+
+ if (!_configured_secondary_variables
+ .emplace(std::make_pair(
+ var_name,
+ SecondaryVariable{
+ var_name, num_components, std::move(fcts)}))
+ .second)
+ {
+ OGS_FATAL("The secondary variable with name `%s' has already been "
+ "set up.",
+ var_name.c_str());
+ }
+ }
+ else if (_all_secondary_variables.find(tag_name) ==
+ _all_secondary_variables.end())
+ {
+ OGS_FATAL("The tag <%s> has not been registered to mark a secondary "
+ "variable.",
+ tag_name.c_str());
+ }
+}
+
+
+} // namespace ProcessLib
diff --git a/ProcessLib/SecondaryVariable.h b/ProcessLib/SecondaryVariable.h
index 480a9feedd2e780d6f7081e98ab279502daa989f..bf99e2e4c3fb4dfef1072c3be335b9caac102063 100644
--- a/ProcessLib/SecondaryVariable.h
+++ b/ProcessLib/SecondaryVariable.h
@@ -21,7 +21,6 @@ namespace ProcessLib
//! Holder for function objects that compute secondary variables,
//! and (optionally) also the residuals (e.g., in case of extrapolation)
-template<typename GlobalVector>
struct SecondaryVariableFunctions final
{
/*! Type of functions used.
@@ -83,18 +82,16 @@ struct SecondaryVariableFunctions final
};
//! Stores information about a specific secondary variable
-template<typename GlobalVector>
struct SecondaryVariable final
{
std::string const name; //!< Name of the variable; used, e.g., for output.
const unsigned n_components; //!< Number of components of the variable.
//! Functions used for computing the secondary variable.
- SecondaryVariableFunctions<GlobalVector> fcts;
+ SecondaryVariableFunctions fcts;
};
//! Handles configuration of several secondary variables from the project file.
-template<typename GlobalVector>
class SecondaryVariableCollection final
{
public:
@@ -109,43 +106,13 @@ public:
*/
SecondaryVariableCollection(
boost::optional<BaseLib::ConfigTree> const& config,
- std::initializer_list<std::string> tag_names)
- {
- if (!config) return;
-
- // read which variables are defined in the config
- for (auto const& tag_name : tag_names) {
- if (!_all_secondary_variables.insert(tag_name).second) {
- OGS_FATAL("Tag name <%s> has been specified twice as a secondary variable.");
- }
-
- //! \ogs_file_special
- if (auto var_name = config->getConfigParameterOptional<std::string>(tag_name))
- {
- // TODO check primary vars, too
- BaseLib::insertIfKeyValueUniqueElseError(
- _map_tagname_to_varname, tag_name, *var_name,
- "Secondary variable names must be unique.");
- }
- }
- }
+ std::initializer_list<std::string> tag_names);
/*! Tells if a secondary variable with the specified name has been set up.
*
* \note \c variable_name is not the tag name in the project file!
*/
- bool variableExists(std::string const& variable_name) const
- {
- auto pred = [&variable_name](std::pair<std::string, std::string> const& p) {
- return p.second == variable_name;
- };
-
- // check if out_var is a secondary variable
- auto const& var = std::find_if(
- _map_tagname_to_varname.cbegin(), _map_tagname_to_varname.cend(), pred);
-
- return var != _map_tagname_to_varname.cend();
- }
+ bool variableExists(std::string const& variable_name) const;
/*! Set up a secondary variable.
*
@@ -161,47 +128,17 @@ public:
*/
void addSecondaryVariable(std::string const& tag_name,
const unsigned num_components,
- SecondaryVariableFunctions<GlobalVector>&& fcts)
- {
- auto it = _map_tagname_to_varname.find(tag_name);
-
- // get user-supplied var_name for the given tag_name
- if (it != _map_tagname_to_varname.end())
- {
- auto const& var_name = it->first;
-
- if (!_configured_secondary_variables
- .emplace(std::make_pair(
- var_name,
- SecondaryVariable<GlobalVector>{
- var_name, num_components, std::move(fcts)}))
- .second)
- {
- OGS_FATAL("The secondary variable with name `%s' has already been "
- "set up.",
- var_name.c_str());
- }
- }
- else if (_all_secondary_variables.find(tag_name) ==
- _all_secondary_variables.end())
- {
- OGS_FATAL("The tag <%s> has not been registered to mark a secondary "
- "variable.",
- tag_name.c_str());
- }
- }
+ SecondaryVariableFunctions&& fcts);
//! Returns an iterator to the first secondary variable.
- typename std::map<std::string,
- SecondaryVariable<GlobalVector>>::const_iterator
+ std::map<std::string, SecondaryVariable>::const_iterator
begin() const
{
return _configured_secondary_variables.begin();
}
//! Returns an iterator past the last secondary variable.
- typename std::map<std::string,
- SecondaryVariable<GlobalVector>>::const_iterator
+ std::map<std::string, SecondaryVariable>::const_iterator
end() const
{
return _configured_secondary_variables.end();
@@ -213,7 +150,7 @@ private:
//! Collection of all configured secondary variables.
//! Maps the variable name to the corresponding SecondaryVariable.
- std::map<std::string, SecondaryVariable<GlobalVector>> _configured_secondary_variables;
+ std::map<std::string, SecondaryVariable> _configured_secondary_variables;
//! Set of all tags available as a secondary variable.
std::set<std::string> _all_secondary_variables;
@@ -222,16 +159,16 @@ private:
//! Creates an object that computes a secondary variable via extrapolation
//! of integration point values.
-template<typename GlobalVector, typename PropertyEnum, typename LocalAssembler>
-SecondaryVariableFunctions<GlobalVector>
+template<typename PropertyEnum, typename LocalAssembler>
+SecondaryVariableFunctions
makeExtrapolator(PropertyEnum const property,
- NumLib::Extrapolator<GlobalVector, PropertyEnum, LocalAssembler>&
+ NumLib::Extrapolator<PropertyEnum, LocalAssembler>&
extrapolator,
- typename NumLib::Extrapolator<GlobalVector, PropertyEnum,
+ typename NumLib::Extrapolator<PropertyEnum,
LocalAssembler>::LocalAssemblers const& local_assemblers)
{
static_assert(std::is_base_of<
- NumLib::Extrapolatable<GlobalVector, PropertyEnum>, LocalAssembler>::value,
+ NumLib::Extrapolatable<PropertyEnum>, LocalAssembler>::value,
"The passed local assembler type (i.e. the local assembler interface) must"
" derive from NumLib::Extrapolatable<>.");
diff --git a/ProcessLib/TES/TESLocalAssembler-impl.h b/ProcessLib/TES/TESLocalAssembler-impl.h
index 4d4aa2ed9712869369ea2fdd6b1d656e6f2a24b4..524089b91208ad2f263fe0cf220e1bd4d5c7eac8 100644
--- a/ProcessLib/TES/TESLocalAssembler-impl.h
+++ b/ProcessLib/TES/TESLocalAssembler-impl.h
@@ -105,9 +105,9 @@ namespace ProcessLib
namespace TES
{
template <typename ShapeFunction_, typename IntegrationMethod_,
- typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
+ unsigned GlobalDim>
TESLocalAssembler<
- ShapeFunction_, IntegrationMethod_, GlobalMatrix, GlobalVector,
+ ShapeFunction_, IntegrationMethod_,
GlobalDim>::TESLocalAssembler(MeshLib::Element const& e,
std::size_t const /*local_matrix_size*/,
unsigned const integration_order,
@@ -131,9 +131,8 @@ TESLocalAssembler<
}
template <typename ShapeFunction_, typename IntegrationMethod_,
- typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
-void TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalMatrix,
- GlobalVector,
+ unsigned GlobalDim>
+void TESLocalAssembler<ShapeFunction_, IntegrationMethod_,
GlobalDim>::assemble(const double /*t*/,
std::vector<double> const& local_x)
{
@@ -186,9 +185,8 @@ void TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalMatrix,
}
template <typename ShapeFunction_, typename IntegrationMethod_,
- typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
-void TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalMatrix,
- GlobalVector, GlobalDim>::
+ unsigned GlobalDim>
+void TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::
addToGlobal(
NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) const
@@ -199,9 +197,9 @@ void TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalMatrix,
}
template <typename ShapeFunction_, typename IntegrationMethod_,
- typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
+ unsigned GlobalDim>
std::vector<double> const& TESLocalAssembler<
- ShapeFunction_, IntegrationMethod_, GlobalMatrix, GlobalVector,
+ ShapeFunction_, IntegrationMethod_,
GlobalDim>::getIntegrationPointValues(TESIntPtVariables const var,
std::vector<double>& cache) const
{
@@ -209,9 +207,9 @@ std::vector<double> const& TESLocalAssembler<
}
template <typename ShapeFunction_, typename IntegrationMethod_,
- typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
+ unsigned GlobalDim>
bool TESLocalAssembler<
- ShapeFunction_, IntegrationMethod_, GlobalMatrix, GlobalVector,
+ ShapeFunction_, IntegrationMethod_,
GlobalDim>::checkBounds(std::vector<double> const& local_x,
std::vector<double> const& local_x_prev_ts)
{
diff --git a/ProcessLib/TES/TESLocalAssembler.h b/ProcessLib/TES/TESLocalAssembler.h
index cc80fc33e7fce097d3e79eef1e7db13b97b2e54c..d3e12d37e7addf6d3c01486198206959891e1384 100644
--- a/ProcessLib/TES/TESLocalAssembler.h
+++ b/ProcessLib/TES/TESLocalAssembler.h
@@ -22,9 +22,8 @@ namespace ProcessLib
{
namespace TES
{
-template <typename GlobalMatrix, typename GlobalVector>
class TESLocalAssemblerInterface
- : public NumLib::Extrapolatable<GlobalVector, TESIntPtVariables>
+ : public NumLib::Extrapolatable<TESIntPtVariables>
{
public:
virtual ~TESLocalAssemblerInterface() = default;
@@ -41,9 +40,9 @@ public:
};
template <typename ShapeFunction_, typename IntegrationMethod_,
- typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
+ unsigned GlobalDim>
class TESLocalAssembler final
- : public TESLocalAssemblerInterface<GlobalMatrix, GlobalVector>
+ : public TESLocalAssemblerInterface
{
public:
using ShapeFunction = ShapeFunction_;
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
index 0d91e6a7c4db4dc2737080825e250fc59ac8b100..f5240ca0f4d97fe2e387188bc698a97f46907dd4 100644
--- a/ProcessLib/TES/TESProcess.cpp
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -9,6 +9,8 @@
#include "TESProcess.h"
+#include "Utils/CreateLocalAssemblers.h"
+
// TODO Copied from VectorMatrixAssembler. Could be provided by the DOF table.
inline NumLib::LocalToGlobalIndexMap::RowColumnIndices
getRowColumnIndices_(std::size_t const id,
@@ -31,7 +33,6 @@ getRowColumnIndices_(std::size_t const id,
indices);
}
-template <typename GlobalVector>
void getVectorValues(
GlobalVector const& x,
NumLib::LocalToGlobalIndexMap::RowColumnIndices const& r_c_indices,
@@ -48,7 +49,6 @@ void getVectorValues(
}
// TODO that essentially duplicates code which is also present in ProcessOutput.
-template <typename GlobalVector>
double getNodalValue(GlobalVector const& x, MeshLib::Mesh const& mesh,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::size_t const node_id,
@@ -73,8 +73,8 @@ TESProcess::TESProcess(
std::unique_ptr<Process::TimeDiscretization>&&
time_discretization,
std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
- SecondaryVariableCollection<GlobalVector>&& secondary_variables,
- ProcessOutput<GlobalVector>&& process_output,
+ SecondaryVariableCollection&& secondary_variables,
+ ProcessOutput&& process_output,
const BaseLib::ConfigTree& config)
: Process(
mesh, nonlinear_solver, std::move(time_discretization),
@@ -202,12 +202,12 @@ void TESProcess::initializeConcreteProcess(
// secondary variables
auto add2nd = [&](std::string const& var_name, unsigned const n_components,
- SecondaryVariableFunctions<GlobalVector>&& fcts) {
+ SecondaryVariableFunctions&& fcts) {
this->_secondary_variables.addSecondaryVariable(var_name, n_components,
std::move(fcts));
};
auto makeEx =
- [&](TESIntPtVariables var) -> SecondaryVariableFunctions<GlobalVector> {
+ [&](TESIntPtVariables var) -> SecondaryVariableFunctions {
return ProcessLib::makeExtrapolator(var, *_extrapolator,
_local_assemblers);
};
@@ -304,7 +304,7 @@ NumLib::IterationResult TESProcess::postIteration(
std::vector<double> local_x_prev_ts_cache;
auto check_variable_bounds = [&](std::size_t id,
- LocalAssembler& loc_asm) {
+ TESLocalAssemblerInterface& loc_asm) {
auto const r_c_indices = getRowColumnIndices_(
id, *this->_local_to_global_index_map, indices_cache);
getVectorValues(x, r_c_indices, local_x_cache);
diff --git a/ProcessLib/TES/TESProcess.h b/ProcessLib/TES/TESProcess.h
index 0570def18bf786b6f773f019fca58419284d845c..58a0e617c5d3cbdd1d3926627d12c2cb5b7db742 100644
--- a/ProcessLib/TES/TESProcess.h
+++ b/ProcessLib/TES/TESProcess.h
@@ -41,8 +41,8 @@ public:
time_discretization,
std::vector<std::reference_wrapper<ProcessVariable>>&&
process_variables,
- SecondaryVariableCollection<GlobalVector>&& secondary_variables,
- ProcessOutput<GlobalVector>&& process_output,
+ SecondaryVariableCollection&& secondary_variables,
+ ProcessOutput&& process_output,
BaseLib::ConfigTree const& config);
void preTimestep(GlobalVector const& x, const double t,
@@ -52,18 +52,15 @@ public:
bool isLinear() const override { return false; }
private:
- using LocalAssembler =
- TESLocalAssemblerInterface<GlobalMatrix, GlobalVector>;
-
using GlobalAssembler = NumLib::VectorMatrixAssembler<
- GlobalMatrix, GlobalVector, LocalAssembler,
+ TESLocalAssemblerInterface,
NumLib::ODESystemTag::FirstOrderImplicitQuasilinear>;
using ExtrapolatorInterface =
- NumLib::Extrapolator<GlobalVector, TESIntPtVariables, LocalAssembler>;
+ NumLib::Extrapolator<TESIntPtVariables, TESLocalAssemblerInterface>;
using ExtrapolatorImplementation =
NumLib::LocalLinearLeastSquaresExtrapolator<
- GlobalVector, TESIntPtVariables, LocalAssembler>;
+ TESIntPtVariables, TESLocalAssemblerInterface>;
void initializeConcreteProcess(
NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -89,7 +86,7 @@ private:
std::unique_ptr<GlobalVector>& result_cache);
std::unique_ptr<GlobalAssembler> _global_assembler;
- std::vector<std::unique_ptr<LocalAssembler>> _local_assemblers;
+ std::vector<std::unique_ptr<TESLocalAssemblerInterface>> _local_assemblers;
AssemblyParams _assembly_params;
@@ -119,7 +116,7 @@ inline std::unique_ptr<TESProcess> createTESProcess(
variables, config,
{"fluid_pressure", "temperature", "vapour_mass_fraction"});
- SecondaryVariableCollection<GlobalVector>
+ SecondaryVariableCollection
secondary_variables{
config.getConfigSubtreeOptional("secondary_variables"),
{"solid_density", "reaction_rate", "velocity_x", "velocity_y",
@@ -127,7 +124,7 @@ inline std::unique_ptr<TESProcess> createTESProcess(
"vapour_partial_pressure", "relative_humidity",
"equilibrium_loading"}};
- ProcessOutput<GlobalVector> process_output{
+ ProcessOutput process_output{
config.getConfigSubtree("output"), process_variables,
secondary_variables};
diff --git a/ProcessLib/Utils/CreateLocalAssemblers.h b/ProcessLib/Utils/CreateLocalAssemblers.h
index f4fae98b4582ab010a9b8c62674b6bf626c96fa3..e14aae36025f633d0d71cc71af825c0a23d45408 100644
--- a/ProcessLib/Utils/CreateLocalAssemblers.h
+++ b/ProcessLib/Utils/CreateLocalAssemblers.h
@@ -24,7 +24,7 @@ namespace detail
{
template<unsigned GlobalDim,
- template <typename, typename, typename, typename, unsigned> class
+ template <typename, typename, unsigned> class
LocalAssemblerImplementation,
typename LocalAssemblerInterface,
typename... ExtraCtorArgs>
@@ -40,8 +40,6 @@ void createLocalAssemblers(
using LocalDataInitializer = LocalDataInitializer<
LocalAssemblerInterface,
LocalAssemblerImplementation,
- GlobalMatrix,
- GlobalVector,
GlobalDim,
ExtraCtorArgs...>;
@@ -74,7 +72,7 @@ void createLocalAssemblers(
* The first two template parameters cannot be deduced from the arguments.
* Therefore they always have to be provided manually.
*/
-template<template <typename, typename, typename, typename, unsigned> class
+template<template <typename, typename, unsigned> class
LocalAssemblerImplementation,
typename LocalAssemblerInterface,
typename... ExtraCtorArgs>
diff --git a/ProcessLib/Utils/LocalDataInitializer.h b/ProcessLib/Utils/LocalDataInitializer.h
index 3f644c31556085363cf65783cded9febd83377f9..5a2ecba8b801110c3728c579553e8bcd19d31bc1 100644
--- a/ProcessLib/Utils/LocalDataInitializer.h
+++ b/ProcessLib/Utils/LocalDataInitializer.h
@@ -122,9 +122,7 @@ namespace ProcessLib
/// NumLib::ShapeLine2 is created.
template <
typename LocalAssemblerInterface,
- template <typename, typename, typename, typename, unsigned> class LocalAssemblerData,
- typename GlobalMatrix,
- typename GlobalVector,
+ template <typename, typename, unsigned> class LocalAssemblerData,
unsigned GlobalDim,
typename... ConstructorArgs>
class LocalDataInitializer final
@@ -285,7 +283,7 @@ private:
using LAData = LocalAssemblerData<
ShapeFunction,
IntegrationMethod<ShapeFunction>,
- GlobalMatrix, GlobalVector, GlobalDim>;
+ GlobalDim>;
/// Generates a function that creates a new LocalAssembler of type LAData<SHAPE_FCT>
template<typename ShapeFct>
diff --git a/Tests/NumLib/SteadyDiffusion2DExample1.h b/Tests/NumLib/SteadyDiffusion2DExample1.h
index 5e72dce9b35b50b7e7ca6ed638feb91127c9ba19..213aed5fac92b35d56953fc3305f164de1195281 100644
--- a/Tests/NumLib/SteadyDiffusion2DExample1.h
+++ b/Tests/NumLib/SteadyDiffusion2DExample1.h
@@ -29,7 +29,6 @@ template<typename IndexType>struct SteadyDiffusion2DExample1
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
using LocalVectorType = Eigen::VectorXd;
- template <typename GlobalMatrix, typename GlobalVector>
class LocalAssemblerData
{
public:
@@ -72,14 +71,13 @@ template<typename IndexType>struct SteadyDiffusion2DExample1
LocalVectorType const* _localRhs = nullptr;
};
- template <typename GlobalMatrix, typename GlobalVector>
static
void initializeLocalData(const MeshLib::Element& e,
- LocalAssemblerData<GlobalMatrix, GlobalVector>*& data_ptr,
+ LocalAssemblerData*& data_ptr,
std::size_t const local_matrix_size,
SteadyDiffusion2DExample1 const& example)
{
- data_ptr = new LocalAssemblerData<GlobalMatrix, GlobalVector>;
+ data_ptr = new LocalAssemblerData;
data_ptr->init(e, local_matrix_size, example._localA, example._localRhs);
}
diff --git a/Tests/NumLib/TestExtrapolation.cpp b/Tests/NumLib/TestExtrapolation.cpp
index 9e509a690683388cc601bf3ca1e8a274897bdee0..afe1fb2667dea4e4f8507ab02a2315fcaa068e20 100644
--- a/Tests/NumLib/TestExtrapolation.cpp
+++ b/Tests/NumLib/TestExtrapolation.cpp
@@ -69,9 +69,8 @@ enum class IntegrationPointValue
DerivedQuantity // a quantity computed for each integration point on-the-fly
};
-template<typename GlobalMatrix, typename GlobalVector>
class LocalAssemblerDataInterface
- : public NumLib::Extrapolatable<GlobalVector, IntegrationPointValue>
+ : public NumLib::Extrapolatable<IntegrationPointValue>
{
public:
virtual void interpolateNodalValuesToIntegrationPoints(
@@ -81,11 +80,9 @@ public:
template<typename ShapeFunction,
typename IntegrationMethod,
- typename GlobalMatrix,
- typename GlobalVector,
unsigned GlobalDim>
class LocalAssemblerData
- : public LocalAssemblerDataInterface<GlobalMatrix, GlobalVector>
+ : public LocalAssemblerDataInterface
{
using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
@@ -150,17 +147,17 @@ private:
class TestProcess
{
public:
- using LocalAssembler = LocalAssemblerDataInterface<GlobalMatrix, GlobalVector>;
+ using LocalAssembler = LocalAssemblerDataInterface;
using GlobalAssembler = NumLib::VectorMatrixAssembler<
- GlobalMatrix, GlobalVector, LocalAssembler,
+ LocalAssembler,
// The exact tag does not matter here.
NumLib::ODESystemTag::FirstOrderImplicitQuasilinear>;
using ExtrapolatorInterface =
- NumLib::Extrapolator<GlobalVector, IntegrationPointValue, LocalAssembler>;
+ NumLib::Extrapolator<IntegrationPointValue, LocalAssembler>;
using ExtrapolatorImplementation =
NumLib::LocalLinearLeastSquaresExtrapolator<
- GlobalVector, IntegrationPointValue, LocalAssembler>;
+ IntegrationPointValue, LocalAssembler>;
TestProcess(MeshLib::Mesh const& mesh, unsigned const integration_order)
: _integration_order(integration_order)
@@ -235,7 +232,7 @@ private:
{
using LocalDataInitializer = ProcessLib::LocalDataInitializer<
LocalAssembler, LocalAssemblerData,
- GlobalMatrix, GlobalVector, GlobalDim>;
+ GlobalDim>;
_local_assemblers.resize(mesh.getNumberOfElements());
diff --git a/Tests/NumLib/TestSerialLinearSolver.cpp b/Tests/NumLib/TestSerialLinearSolver.cpp
index cc495b12b777ed9d9bf9140e6ae805529e3db696..14fac011b9f7f0800a952c516d8609fa443e9551 100644
--- a/Tests/NumLib/TestSerialLinearSolver.cpp
+++ b/Tests/NumLib/TestSerialLinearSolver.cpp
@@ -73,7 +73,7 @@ TEST(NumLibSerialLinearSolver, Steady2DdiffusionQuadElem)
auto x = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(ms);
// TODO no setZero() for rhs, x?
- using LocalAssembler = Example::LocalAssemblerData<GlobalMatrix, GlobalVector>;
+ using LocalAssembler = Example::LocalAssemblerData;
// Initializer of the local assembler data.
std::vector<LocalAssembler*> local_assembler_data;
local_assembler_data.resize(ex1.msh->getNumberOfElements());
@@ -87,7 +87,7 @@ TEST(NumLibSerialLinearSolver, Steady2DdiffusionQuadElem)
auto const num_local_dof = local_to_global_index_map.getNumElementDOF(id);
- Example::initializeLocalData<GlobalMatrix, GlobalVector>(
+ Example::initializeLocalData(
item, item_data, num_local_dof, ex1);
};
@@ -100,7 +100,7 @@ TEST(NumLibSerialLinearSolver, Steady2DdiffusionQuadElem)
// TODO in the future use simpler NumLib::ODESystemTag
// Local and global assemblers.
typedef NumLib::VectorMatrixAssembler<
- GlobalMatrix, GlobalVector, LocalAssembler,
+ LocalAssembler,
NumLib::ODESystemTag::FirstOrderImplicitQuasilinear> GlobalAssembler;
GlobalAssembler assembler(local_to_global_index_map);