diff --git a/ProcessLib/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlowProcess.h
index b9bbdf2df5a755663683a24f43e788d79e26efe4..c4db0b63676e4e9bc16c6f0f74797c24e45d938f 100644
--- a/ProcessLib/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlowProcess.h
@@ -11,30 +11,15 @@
#define PROCESS_LIB_GROUNDWATERFLOWPROCESS_H_
#include <cassert>
-#include <memory>
-#include <boost/algorithm/string/erase.hpp>
#include <boost/optional.hpp>
-
-#include "logog/include/logog.hpp"
-
#include "AssemblerLib/LocalAssemblerBuilder.h"
#include "AssemblerLib/LocalDataInitializer.h"
#include "FileIO/VtkIO/VtuInterface.h"
-#include "MathLib/LinAlg/ApplyKnownSolution.h"
-
-#include "MeshLib/MeshSubset.h"
-#include "MeshGeoToolsLib/MeshNodeSearcher.h"
-
-#include "UniformDirichletBoundaryCondition.h"
-
#include "GroundwaterFlowFEM.h"
-#include "NeumannBcAssembler.h"
-#include "NeumannBc.h"
-#include "DirichletBc.h"
#include "Parameter.h"
#include "Process.h"
@@ -51,8 +36,6 @@ namespace ProcessLib
template<typename GlobalSetup>
class GroundwaterFlowProcess : public Process<GlobalSetup>
{
- unsigned const _integration_order = 2;
-
public:
GroundwaterFlowProcess(MeshLib::Mesh& mesh,
std::vector<ProcessVariable> const& variables,
@@ -80,7 +63,8 @@ public:
DBUG("Associate hydraulic_head with process variable \'%s\'.",
name.c_str());
- _hydraulic_head = const_cast<ProcessVariable*>(&*variable);
+ this->_process_variables.emplace_back(
+ const_cast<ProcessVariable*>(&*variable));
}
// Hydraulic conductivity parameter.
@@ -152,51 +136,7 @@ public:
this->_mesh.getElements(),
_local_assemblers,
*_hydraulic_conductivity,
- _integration_order);
-
- DBUG("Initialize boundary conditions.");
- MeshGeoToolsLib::MeshNodeSearcher& hydraulic_head_mesh_node_searcher =
- MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher(
- _hydraulic_head->getMesh());
-
- _hydraulic_head->initializeDirichletBCs(
- std::back_inserter(_dirichlet_bcs),
- hydraulic_head_mesh_node_searcher,
- *this->_local_to_global_index_map,
- 0);
-
- //
- // Neumann boundary conditions.
- //
- {
- // Find mesh nodes.
- MeshGeoToolsLib::BoundaryElementsSearcher hydraulic_head_mesh_element_searcher(
- _hydraulic_head->getMesh(), hydraulic_head_mesh_node_searcher);
-
- // Create a neumann BC for the hydraulic head storing them in the
- // _neumann_bcs vector.
- _hydraulic_head->createNeumannBcs(
- std::back_inserter(_neumann_bcs),
- hydraulic_head_mesh_element_searcher,
- this->_global_setup,
- _integration_order,
- *this->_local_to_global_index_map,
- 0,
- *_mesh_subset_all_nodes);
- }
-
- Process<GlobalSetup>::initializeNeumannBcs(_neumann_bcs);
- }
-
- void initializeMeshSubsets(MeshLib::Mesh const& mesh) override
- {
- // Create single component dof in every of the mesh's nodes.
- _mesh_subset_all_nodes =
- new MeshLib::MeshSubset(mesh, &mesh.getNodes());
-
- // Collect the mesh subsets in a vector.
- this->_all_mesh_subsets.push_back(
- new MeshLib::MeshSubsets(_mesh_subset_all_nodes));
+ this->_integration_order);
}
std::string getLinearSolverName() const override
@@ -204,12 +144,8 @@ public:
return "gw_";
}
- void init() override
+ void createLocalAssemblers() override
{
- DBUG("Initialize GroundwaterFlowProcess.");
-
- Process<GlobalSetup>::setInitialConditions(*_hydraulic_head, 0);
-
if (this->_mesh.getDimension()==1)
createLocalAssemblers<1>();
else if (this->_mesh.getDimension()==2)
@@ -230,14 +166,6 @@ public:
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(this->_global_setup);
-
- for (auto const& bc : _dirichlet_bcs)
- MathLib::applyKnownSolution(*this->_A, *this->_rhs, *this->_x,
- bc.global_ids, bc.values);
-
return true;
}
@@ -294,29 +222,17 @@ public:
~GroundwaterFlowProcess()
{
- for (auto p : _neumann_bcs)
- delete p;
-
for (auto p : _local_assemblers)
delete p;
-
- delete _mesh_subset_all_nodes;
}
private:
- ProcessVariable* _hydraulic_head = nullptr;
-
Parameter<double, MeshLib::Element const&> const* _hydraulic_conductivity = nullptr;
- MeshLib::MeshSubset const* _mesh_subset_all_nodes = nullptr;
-
using LocalAssembler = GroundwaterFlow::LocalAssemblerDataInterface<
typename GlobalSetup::MatrixType, typename GlobalSetup::VectorType>;
std::vector<LocalAssembler*> _local_assemblers;
-
- std::vector<DirichletBc<GlobalIndexType>> _dirichlet_bcs;
- std::vector<NeumannBc<GlobalSetup>*> _neumann_bcs;
};
} // namespace ProcessLib
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 9de6e4ac5a50206e2377a2bb7861c542afc18681..fc86eb09239410556da60805634f067c2e64dc10 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -10,13 +10,20 @@
#ifndef PROCESS_LIB_PROCESS_H_
#define PROCESS_LIB_PROCESS_H_
+#include <memory>
#include <string>
+#include <logog/include/logog.hpp>
+
+
#include "AssemblerLib/ComputeSparsityPattern.h"
-#include "AssemblerLib/VectorMatrixAssembler.h"
#include "AssemblerLib/LocalToGlobalIndexMap.h"
+#include "AssemblerLib/VectorMatrixAssembler.h"
#include "BaseLib/ConfigTreeNew.h"
+#include "MathLib/LinAlg/ApplyKnownSolution.h"
#include "MathLib/LinAlg/SetMatrixSparsity.h"
+#include "MeshGeoToolsLib/MeshNodeSearcher.h"
+#include "MeshLib/MeshSubset.h"
#include "MeshLib/MeshSubsets.h"
#ifdef USE_PETSC
@@ -24,7 +31,11 @@
#include "MathLib/LinAlg/PETSc/PETScMatrixOption.h"
#endif
+#include "DirichletBc.h"
+#include "NeumannBc.h"
+#include "NeumannBcAssembler.h"
#include "ProcessVariable.h"
+#include "UniformDirichletBoundaryCondition.h"
namespace MeshLib
{
@@ -42,10 +53,11 @@ public:
{
for (auto p : _all_mesh_subsets)
delete p;
+ delete _mesh_subset_all_nodes;
}
/// Process specific initialization called by initialize().
- virtual void init() = 0;
+ virtual void createLocalAssemblers() = 0;
virtual bool assemble(const double delta_t) = 0;
virtual std::string getLinearSolverName() const = 0;
@@ -56,15 +68,12 @@ public:
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);
+ initializeMeshSubsets();
_local_to_global_index_map.reset(
new AssemblerLib::LocalToGlobalIndexMap(
@@ -82,12 +91,20 @@ public:
_global_assembler.reset(
new GlobalAssembler(*_A, *_rhs, *_local_to_global_index_map));
- init(); // Execute proces specific initialization.
- }
+ createLocalAssemblers();
- void initializeNeumannBcs(std::vector<NeumannBc<GlobalSetup>*> const& bcs)
- {
- for (auto bc : bcs)
+ for (auto const& pv : _process_variables)
+ {
+ DBUG("Set initial conditions.");
+ setInitialConditions(*pv, 0); // 0 is the component id
+
+ DBUG("Initialize boundary conditions.");
+ createDirichletBcs(*pv, 0); // 0 is the component id
+
+ createNeumannBcs(*pv, 0); // 0 is the component id
+ }
+
+ for (auto& bc : _neumann_bcs)
bc->initialize(_global_setup, *_A, *_rhs, _mesh.getDimension());
}
@@ -98,6 +115,14 @@ public:
bool const result = assemble(delta_t);
+ // Call global assembler for each Neumann boundary local assembler.
+ for (auto const& bc : _neumann_bcs)
+ bc->integrate(_global_setup);
+
+ for (auto const& bc : _dirichlet_bcs)
+ MathLib::applyKnownSolution(*_A, *_rhs, *_x,
+ bc.global_ids, bc.values);
+
_linear_solver->solve(*_rhs, *_x);
return result;
}
@@ -111,6 +136,19 @@ protected:
std::move(config)));
}
+private:
+ /// Creates mesh subsets, i.e. components, for given mesh.
+ void initializeMeshSubsets()
+ {
+ // Create single component dof in every of the mesh's nodes.
+ _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));
+ }
+
/// Sets the initial condition values in the solution vector x for a given
/// process variable and component.
void setInitialConditions(ProcessVariable const& variable,
@@ -128,7 +166,39 @@ protected:
}
}
-private:
+ void createDirichletBcs(ProcessVariable& variable,
+ 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,
+ component_id);
+ }
+
+ void createNeumannBcs(ProcessVariable& variable, 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,
+ _global_setup,
+ _integration_order,
+ *_local_to_global_index_map,
+ component_id,
+ *_mesh_subset_all_nodes);
+ }
+
/// Creates global matrix, rhs and solution vectors, and the linear solver.
void createLinearSolver(std::string const& solver_name)
{
@@ -162,7 +232,10 @@ private:
}
protected:
+ unsigned const _integration_order = 2;
+
MeshLib::Mesh& _mesh;
+ MeshLib::MeshSubset const* _mesh_subset_all_nodes = nullptr;
std::vector<MeshLib::MeshSubsets*> _all_mesh_subsets;
GlobalSetup _global_setup;
@@ -184,6 +257,12 @@ protected:
std::unique_ptr<typename GlobalSetup::VectorType> _x;
AssemblerLib::SparsityPattern _sparsity_pattern;
+
+ std::vector<DirichletBc<GlobalIndexType>> _dirichlet_bcs;
+ std::vector<std::unique_ptr<NeumannBc<GlobalSetup>>> _neumann_bcs;
+
+ /// Variables used by this process.
+ std::vector<ProcessVariable*> _process_variables;
};
} // namespace ProcessLib
diff --git a/ProcessLib/ProcessVariable.h b/ProcessLib/ProcessVariable.h
index f67506765f04301d32c642058ceb3feacfff25ec..565768a807075e5a47bef920bd8d4e1b696f4d2a 100644
--- a/ProcessLib/ProcessVariable.h
+++ b/ProcessLib/ProcessVariable.h
@@ -81,8 +81,9 @@ public:
for (auto& config : _neumann_bc_configs)
{
config->initialize(searcher);
- bcs++ = new NeumannBc<GlobalSetup>(*config,
- std::forward<Args>(args)...);
+ bcs++ = std::unique_ptr<NeumannBc<GlobalSetup>>{
+ new NeumannBc<GlobalSetup>(*config,
+ std::forward<Args>(args)...)};
}
}