diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 25bd20b0ffffdec4a2479a728a629d470f6cbd31..2a072bbb88a316320fe67b71f608960547bec21c 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -84,9 +84,6 @@ ProjectData::~ProjectData()
{
delete _geoObjects;
- for(ProcessLib::Process* p : _processes)
- delete p;
-
for (MeshLib::Mesh* m : _mesh_vec)
delete m;
}
diff --git a/Applications/ApplicationsLib/ProjectData.h b/Applications/ApplicationsLib/ProjectData.h
index 2fbf73f076049e162d7e6a6cb92a828652a9d57c..444169fc360beef9747f9418e5d841852e475284 100644
--- a/Applications/ApplicationsLib/ProjectData.h
+++ b/Applications/ApplicationsLib/ProjectData.h
@@ -97,7 +97,7 @@ public:
// TODO at the moment we have only one mesh, later there can be
// several meshes. Then we have to assign the referenced mesh
// here.
- _processes.push_back(
+ _processes.emplace_back(
new ProcessLib::GroundwaterFlowProcess<GlobalSetupType>(
*_mesh_vec[0], _process_variables, _parameters, pc));
}
@@ -110,24 +110,26 @@ public:
/// Iterator access for processes.
/// Provides read access to the process container.
- std::vector<ProcessLib::Process*>::const_iterator
+ std::vector<
+ std::unique_ptr<ProcessLib::Process<GlobalSetupType>>>::const_iterator
processesBegin() const
{
return _processes.begin();
}
- std::vector<ProcessLib::Process*>::iterator
+ std::vector<std::unique_ptr<ProcessLib::Process<GlobalSetupType>>>::iterator
processesBegin()
{
return _processes.begin();
}
/// Iterator access for processes as in processesBegin().
- std::vector<ProcessLib::Process*>::const_iterator
+ std::vector<
+ std::unique_ptr<ProcessLib::Process<GlobalSetupType>>>::const_iterator
processesEnd() const
{
return _processes.end();
}
- std::vector<ProcessLib::Process*>::iterator
+ std::vector<std::unique_ptr<ProcessLib::Process<GlobalSetupType>>>::iterator
processesEnd()
{
return _processes.end();
@@ -187,7 +189,8 @@ private:
private:
GeoLib::GEOObjects *_geoObjects = new GeoLib::GEOObjects();
std::vector<MeshLib::Mesh*> _mesh_vec;
- std::vector<ProcessLib::Process*> _processes;
+ std::vector<std::unique_ptr<ProcessLib::Process<GlobalSetupType>>>
+ _processes;
std::vector<ProcessLib::ProcessVariable> _process_variables;
/// Buffer for each process' config used in the process building function.
diff --git a/ProcessLib/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlowProcess.h
index 097b85c1d1727df527ce4ddda3a18fc964cb9b06..8f91add480230a7f7936cfa07f72067e2f41712b 100644
--- a/ProcessLib/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlowProcess.h
@@ -13,30 +13,20 @@
#include <cassert>
#include <memory>
-#include <boost/optional.hpp>
#include <boost/algorithm/string/erase.hpp>
+#include <boost/optional.hpp>
-#include "logog/include/logog.hpp"
-#include "BaseLib/ConfigTree.h"
-#ifdef USE_PETSC
-#include "MeshLib/NodePartitionedMesh.h"
-#include "MathLib/LinAlg/PETSc/PETScMatrixOption.h"
-#endif
+#include "logog/include/logog.hpp"
#include "AssemblerLib/LocalAssemblerBuilder.h"
-#include "AssemblerLib/VectorMatrixAssembler.h"
#include "AssemblerLib/LocalDataInitializer.h"
-#include "AssemblerLib/LocalToGlobalIndexMap.h"
-#include "AssemblerLib/ComputeSparsityPattern.h"
#include "FileIO/VtkIO/VtuInterface.h"
#include "MathLib/LinAlg/ApplyKnownSolution.h"
-#include "MathLib/LinAlg/SetMatrixSparsity.h"
#include "MeshLib/MeshSubset.h"
-#include "MeshLib/MeshSubsets.h"
#include "MeshGeoToolsLib/MeshNodeSearcher.h"
#include "UniformDirichletBoundaryCondition.h"
@@ -46,7 +36,6 @@
#include "NeumannBc.h"
#include "Parameter.h"
#include "Process.h"
-#include "ProcessVariable.h"
namespace MeshLib
{
@@ -59,7 +48,7 @@ namespace ProcessLib
{
template<typename GlobalSetup>
-class GroundwaterFlowProcess : public Process
+class GroundwaterFlowProcess : public Process<GlobalSetup>
{
unsigned const _integration_order = 2;
@@ -68,7 +57,7 @@ public:
std::vector<ProcessVariable> const& variables,
std::vector<std::unique_ptr<ParameterBase>> const& parameters,
BaseLib::ConfigTree const& config)
- : Process(mesh)
+ : Process<GlobalSetup>(mesh)
{
DBUG("Create GroundwaterFlowProcess.");
@@ -131,14 +120,10 @@ public:
}
// Linear solver options
- {
- auto const par = config.get_child_optional("linear_solver");
-
- if (par)
- {
- _linear_solver_options.reset(new BaseLib::ConfigTree(*par));
- }
- }
+ if (auto const& linear_solver_options =
+ config.get_child_optional("linear_solver"))
+ Process<GlobalSetup>::setLinearSolverOptions(
+ *linear_solver_options);
}
template <unsigned GlobalDim>
@@ -146,7 +131,7 @@ public:
{
DBUG("Create local assemblers.");
// Populate the vector of local assemblers.
- _local_assemblers.resize(_mesh.getNElements());
+ _local_assemblers.resize(this->_mesh.getNElements());
// Shape matrices initializer
using LocalDataInitializer = AssemblerLib::LocalDataInitializer<
GroundwaterFlow::LocalAssemblerDataInterface,
@@ -163,20 +148,16 @@ public:
LocalDataInitializer>;
LocalAssemblerBuilder local_asm_builder(
- initializer, *_local_to_global_index_map);
+ initializer, *this->_local_to_global_index_map);
DBUG("Calling local assembler builder for all mesh elements.");
- _global_setup.execute(
+ this->_global_setup.execute(
local_asm_builder,
- _mesh.getElements(),
+ this->_mesh.getElements(),
_local_assemblers,
*_hydraulic_conductivity,
_integration_order);
- DBUG("Create global assembler.");
- _global_assembler.reset(
- new GlobalAssembler(*_A, *_rhs, *_local_to_global_index_map));
-
DBUG("Initialize boundary conditions.");
MeshGeoToolsLib::MeshNodeSearcher& hydraulic_head_mesh_node_searcher =
MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher(
@@ -184,7 +165,7 @@ public:
_hydraulic_head->initializeDirichletBCs(
hydraulic_head_mesh_node_searcher,
- *_local_to_global_index_map, 0,
+ *this->_local_to_global_index_map, 0,
_dirichlet_bc.global_ids, _dirichlet_bc.values);
//
@@ -200,102 +181,65 @@ public:
_hydraulic_head->createNeumannBcs(
std::back_inserter(_neumann_bcs),
hydraulic_head_mesh_element_searcher,
- _global_setup,
+ this->_global_setup,
_integration_order,
- *_local_to_global_index_map,
+ *this->_local_to_global_index_map,
0,
*_mesh_subset_all_nodes);
}
- for (auto bc : _neumann_bcs)
- bc->initialize(_global_setup, *_A, *_rhs, _mesh.getDimension());
-
+ Process<GlobalSetup>::initializeNeumannBcs(_neumann_bcs);
}
- void initialize() override
+ void initializeMeshSubsets(MeshLib::Mesh const& mesh) override
{
- DBUG("Initialize GroundwaterFlowProcess.");
-
- DBUG("Construct dof mappings.");
// Create single component dof in every of the mesh's nodes.
- _mesh_subset_all_nodes = new MeshLib::MeshSubset(_mesh, &_mesh.getNodes());
+ _mesh_subset_all_nodes =
+ new MeshLib::MeshSubset(mesh, &mesh.getNodes());
// Collect the mesh subsets in a vector.
- _all_mesh_subsets.push_back(new MeshLib::MeshSubsets(_mesh_subset_all_nodes));
-
- _local_to_global_index_map.reset(
- new AssemblerLib::LocalToGlobalIndexMap(_all_mesh_subsets, AssemblerLib::ComponentOrder::BY_COMPONENT));
+ this->_all_mesh_subsets.push_back(
+ new MeshLib::MeshSubsets(_mesh_subset_all_nodes));
+ }
-#ifdef USE_PETSC
- DBUG("Allocate global matrix, vectors, and linear solver.");
- MathLib::PETScMatrixOption mat_opt;
- const MeshLib::NodePartitionedMesh& pmesh =
- static_cast<const MeshLib::NodePartitionedMesh&>(_mesh);
- mat_opt.d_nz = pmesh.getMaximumNConnectedNodesToNode();
- mat_opt.o_nz = mat_opt.d_nz;
- const std::size_t num_unknowns =
- _local_to_global_index_map->dofSizeGlobal();
- _A.reset(_global_setup.createMatrix(num_unknowns, mat_opt));
-#else
- DBUG("Compute sparsity pattern");
- _sparsity_pattern = std::move(
- AssemblerLib::computeSparsityPattern(
- *_local_to_global_index_map, _mesh));
-
- DBUG("Allocate global matrix, vectors, and linear solver.");
- const std::size_t num_unknowns = _local_to_global_index_map->dofSize();
- _A.reset(_global_setup.createMatrix(num_unknowns));
-#endif
+ std::string getLinearSolverName() const override
+ {
+ return "gw_";
+ }
- _x.reset(_global_setup.createVector(num_unknowns));
- _rhs.reset(_global_setup.createVector(num_unknowns));
- _linear_solver.reset(new typename GlobalSetup::LinearSolver(
- *_A, "gw_", _linear_solver_options.get()));
+ void init() override
+ {
+ DBUG("Initialize GroundwaterFlowProcess.");
- setInitialConditions(*_hydraulic_head);
+ Process<GlobalSetup>::setInitialConditions(*_hydraulic_head, 0);
- if (_mesh.getDimension()==1)
+ if (this->_mesh.getDimension()==1)
createLocalAssemblers<1>();
- else if (_mesh.getDimension()==2)
+ else if (this->_mesh.getDimension()==2)
createLocalAssemblers<2>();
- else if (_mesh.getDimension()==3)
+ else if (this->_mesh.getDimension()==3)
createLocalAssemblers<3>();
else
assert(false);
}
- void setInitialConditions(ProcessVariable const& variable)
- {
- std::size_t const n = _mesh.getNNodes();
- for (std::size_t i = 0; i < n; ++i)
- {
- MeshLib::Location const l(_mesh.getID(),
- MeshLib::MeshItemType::Node, i);
- auto const global_index = // 0 is the component id.
- std::abs( _local_to_global_index_map->getGlobalIndex(l, 0) );
- _x->set(global_index,
- variable.getInitialConditionValue(*_mesh.getNode(i)));
- }
- }
-
- bool solve(const double /*delta_t*/) override
+ bool assemble(const double /*delta_t*/) override
{
- DBUG("Solve GroundwaterFlowProcess.");
+ DBUG("Assemble GroundwaterFlowProcess.");
- _A->setZero();
- MathLib::setMatrixSparsity(*_A, _sparsity_pattern);
- *_rhs = 0; // This resets the whole vector.
+ *this->_rhs = 0; // This resets the whole vector.
// Call global assembler for each local assembly item.
- _global_setup.execute(*_global_assembler, _local_assemblers);
+ this->_global_setup.execute(*this->_global_assembler,
+ _local_assemblers);
// Call global assembler for each Neumann boundary local assembler.
for (auto bc : _neumann_bcs)
- bc->integrate(_global_setup);
+ bc->integrate(this->_global_setup);
- MathLib::applyKnownSolution(*_A, *_rhs, *_x, _dirichlet_bc.global_ids,
- _dirichlet_bc.values);
- _linear_solver->solve(*_rhs, *_x);
+ MathLib::applyKnownSolution(*this->_A, *this->_rhs, *this->_x,
+ _dirichlet_bc.global_ids,
+ _dirichlet_bc.values);
return true;
}
@@ -308,41 +252,41 @@ public:
// Get or create a property vector for results.
boost::optional<MeshLib::PropertyVector<double>&> result;
- if (_mesh.getProperties().hasPropertyVector(property_name))
+ if (this->_mesh.getProperties().hasPropertyVector(property_name))
{
- result = _mesh.getProperties().template
+ result = this->_mesh.getProperties().template
getPropertyVector<double>(property_name);
}
else
{
- result = _mesh.getProperties().template
+ result = this->_mesh.getProperties().template
createNewPropertyVector<double>(property_name,
MeshLib::MeshItemType::Node);
- result->resize(_x->size());
+ result->resize(this->_x->size());
}
- assert(result && result->size() == _x->size());
+ assert(result && result->size() == this->_x->size());
#ifdef USE_PETSC
- std::unique_ptr<double[]> u(new double[_x->size()]);
- _x->getGlobalVector(u.get()); // get the global solution
+ std::unique_ptr<double[]> u(new double[this->_x->size()]);
+ this->_x->getGlobalVector(u.get()); // get the global solution
- std::size_t const n = _mesh.getNNodes();
+ std::size_t const n = this->_mesh.getNNodes();
for (std::size_t i = 0; i < n; ++i)
{
- MeshLib::Location const l(_mesh.getID(),
+ MeshLib::Location const l(this->_mesh.getID(),
MeshLib::MeshItemType::Node, i);
auto const global_index = std::abs( // 0 is the component id.
- _local_to_global_index_map->getGlobalIndex(l, 0));
+ this->_local_to_global_index_map->getGlobalIndex(l, 0));
(*result)[i] = u[global_index];
}
#else
// Copy result
- for (std::size_t i = 0; i < _x->size(); ++i)
- (*result)[i] = (*_x)[i];
+ for (std::size_t i = 0; i < this->_x->size(); ++i)
+ (*result)[i] = (*this->_x)[i];
#endif
// Write output file
- FileIO::VtuInterface vtu_interface(&_mesh, vtkXMLWriter::Binary, true);
+ FileIO::VtuInterface vtu_interface(&this->_mesh, vtkXMLWriter::Binary, true);
vtu_interface.writeToFile(file_name);
}
@@ -359,9 +303,6 @@ public:
for (auto p : _local_assemblers)
delete p;
- for (auto p : _all_mesh_subsets)
- delete p;
-
delete _mesh_subset_all_nodes;
}
@@ -371,30 +312,12 @@ private:
Parameter<double, MeshLib::Element const&> const* _hydraulic_conductivity = nullptr;
MeshLib::MeshSubset const* _mesh_subset_all_nodes = nullptr;
- std::vector<MeshLib::MeshSubsets*> _all_mesh_subsets;
-
- GlobalSetup _global_setup;
- std::unique_ptr<BaseLib::ConfigTree> _linear_solver_options;
- std::unique_ptr<typename GlobalSetup::LinearSolver> _linear_solver;
- std::unique_ptr<typename GlobalSetup::MatrixType> _A;
- std::unique_ptr<typename GlobalSetup::VectorType> _rhs;
- std::unique_ptr<typename GlobalSetup::VectorType> _x;
using LocalAssembler = GroundwaterFlow::LocalAssemblerDataInterface<
typename GlobalSetup::MatrixType, typename GlobalSetup::VectorType>;
std::vector<LocalAssembler*> _local_assemblers;
- using GlobalAssembler =
- AssemblerLib::VectorMatrixAssembler<
- typename GlobalSetup::MatrixType,
- typename GlobalSetup::VectorType>;
-
-
- std::unique_ptr<AssemblerLib::LocalToGlobalIndexMap> _local_to_global_index_map;
-
- std::unique_ptr<GlobalAssembler> _global_assembler;
-
/// Global ids in the global matrix/vector where the dirichlet bc is
/// imposed and their corresponding values.
struct DirichletBC {
@@ -403,8 +326,6 @@ private:
} _dirichlet_bc;
std::vector<NeumannBc<GlobalSetup>*> _neumann_bcs;
-
- AssemblerLib::SparsityPattern _sparsity_pattern;
};
} // namespace ProcessLib
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 47a09b12737b3a57abe7e80919956382a8d6c971..2a7f7fc75464aa9db4548462ebdd2e4e67a09841 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -12,35 +12,178 @@
#include <string>
+#include "AssemblerLib/ComputeSparsityPattern.h"
+#include "AssemblerLib/VectorMatrixAssembler.h"
+#include "AssemblerLib/LocalToGlobalIndexMap.h"
+#include "BaseLib/ConfigTree.h"
+#include "MathLib/LinAlg/SetMatrixSparsity.h"
+#include "MeshLib/MeshSubsets.h"
+
+#ifdef USE_PETSC
+#include "MeshLib/NodePartitionedMesh.h"
+#include "MathLib/LinAlg/PETSc/PETScMatrixOption.h"
+#endif
+
+#include "ProcessVariable.h"
+
namespace MeshLib
{
- class Mesh;
+class Mesh;
}
namespace ProcessLib
{
-
+template <typename GlobalSetup>
class Process
{
public:
- Process(MeshLib::Mesh& mesh)
- : _mesh(mesh)
- { }
+ Process(MeshLib::Mesh& mesh) : _mesh(mesh) {}
+ virtual ~Process()
+ {
+ for (auto p : _all_mesh_subsets)
+ delete p;
+ }
+
+ /// Process specific initialization called by initialize().
+ virtual void init() = 0;
+ virtual bool assemble(const double delta_t) = 0;
+
+ virtual std::string getLinearSolverName() const = 0;
+
+ /// Postprocessing after solve().
+ /// The file_name is indicating the name of possible output file.
+ virtual void post(std::string const& file_name) = 0;
+ virtual void postTimestep(std::string const& file_name,
+ const unsigned timestep) = 0;
+
+ /// Creates mesh subsets, i.e. components, for given mesh.
+ virtual void initializeMeshSubsets(MeshLib::Mesh const& mesh) = 0;
+
+ void initialize()
+ {
+ DBUG("Initialize process.");
+
+ DBUG("Construct dof mappings.");
+ initializeMeshSubsets(_mesh);
+
+ _local_to_global_index_map.reset(
+ new AssemblerLib::LocalToGlobalIndexMap(
+ _all_mesh_subsets, AssemblerLib::ComponentOrder::BY_COMPONENT));
+
+#ifndef USE_PETSC
+ DBUG("Compute sparsity pattern");
+ computeSparsityPattern();
+#endif
+
+ // create global vectors and linear solver
+ createLinearSolver(getLinearSolverName());
+
+ DBUG("Create global assembler.");
+ _global_assembler.reset(
+ new GlobalAssembler(*_A, *_rhs, *_local_to_global_index_map));
- virtual ~Process() = default;
+ init(); // Execute proces specific initialization.
+ }
- virtual void initialize() = 0;
- virtual bool solve(const double delta_t) = 0;
+ void initializeNeumannBcs(std::vector<NeumannBc<GlobalSetup>*> const& bcs)
+ {
+ for (auto bc : bcs)
+ bc->initialize(_global_setup, *_A, *_rhs, _mesh.getDimension());
+ }
- /// Postprocessing after solve().
- /// The file_name is indicating the name of possible output file.
- virtual void post(std::string const& file_name) = 0;
- virtual void postTimestep(std::string const& file_name, const unsigned timestep) = 0;
+ bool solve(const double delta_t)
+ {
+ _A->setZero();
+ MathLib::setMatrixSparsity(*_A, _sparsity_pattern);
+
+ bool const result = assemble(delta_t);
+
+ _linear_solver->solve(*_rhs, *_x);
+ return result;
+ }
+
+protected:
+ /// Set linear solver options; called by the derived process which is
+ /// parsing the configuration.
+ void setLinearSolverOptions(const BaseLib::ConfigTree& config)
+ {
+ _linear_solver_options.reset(new BaseLib::ConfigTree(config));
+ }
+
+ /// Sets the initial condition values in the solution vector x for a given
+ /// process variable and component.
+ void setInitialConditions(ProcessVariable const& variable,
+ int const component_id)
+ {
+ std::size_t const n = _mesh.getNNodes();
+ for (std::size_t i = 0; i < n; ++i)
+ {
+ MeshLib::Location const l(_mesh.getID(),
+ MeshLib::MeshItemType::Node, i);
+ auto const global_index = std::abs(
+ _local_to_global_index_map->getGlobalIndex(l, component_id));
+ _x->set(global_index,
+ variable.getInitialConditionValue(*_mesh.getNode(i)));
+ }
+ }
+
+private:
+ /// Creates global matrix, rhs and solution vectors, and the linear solver.
+ void createLinearSolver(std::string const& solver_name)
+ {
+ DBUG("Allocate global matrix, vectors, and linear solver.");
+#ifdef USE_PETSC
+ MathLib::PETScMatrixOption mat_opt;
+ const MeshLib::NodePartitionedMesh& pmesh =
+ static_cast<const MeshLib::NodePartitionedMesh&>(_mesh);
+ mat_opt.d_nz = pmesh.getMaximumNConnectedNodesToNode();
+ mat_opt.o_nz = mat_opt.d_nz;
+ const std::size_t num_unknowns =
+ _local_to_global_index_map->dofSizeGlobal();
+ _A.reset(_global_setup.createMatrix(num_unknowns, mat_opt));
+#else
+ const std::size_t num_unknowns = _local_to_global_index_map->dofSize();
+ _A.reset(_global_setup.createMatrix(num_unknowns));
+#endif
+ _x.reset(_global_setup.createVector(num_unknowns));
+ _rhs.reset(_global_setup.createVector(num_unknowns));
+ _linear_solver.reset(new typename GlobalSetup::LinearSolver(
+ *_A, solver_name, _linear_solver_options.get()));
+ }
+
+ /// Computes and stores global matrix' sparsity pattern from given
+ /// DOF-table.
+ void computeSparsityPattern()
+ {
+ _sparsity_pattern = std::move(AssemblerLib::computeSparsityPattern(
+ *_local_to_global_index_map, _mesh));
+ }
protected:
- MeshLib::Mesh& _mesh;
+ MeshLib::Mesh& _mesh;
+ std::vector<MeshLib::MeshSubsets*> _all_mesh_subsets;
+
+ GlobalSetup _global_setup;
+
+ using GlobalAssembler =
+ AssemblerLib::VectorMatrixAssembler<typename GlobalSetup::MatrixType,
+ typename GlobalSetup::VectorType>;
+
+ std::unique_ptr<GlobalAssembler> _global_assembler;
+
+ std::unique_ptr<AssemblerLib::LocalToGlobalIndexMap>
+ _local_to_global_index_map;
+
+ std::unique_ptr<BaseLib::ConfigTree> _linear_solver_options;
+ std::unique_ptr<typename GlobalSetup::LinearSolver> _linear_solver;
+
+ std::unique_ptr<typename GlobalSetup::MatrixType> _A;
+ std::unique_ptr<typename GlobalSetup::VectorType> _rhs;
+ std::unique_ptr<typename GlobalSetup::VectorType> _x;
+
+ AssemblerLib::SparsityPattern _sparsity_pattern;
};
-} // namespace ProcessLib
+} // namespace ProcessLib
#endif // PROCESS_LIB_PROCESS_H_