Skip to content
Snippets Groups Projects
Commit c8376b44 authored by Dmitri Naumov's avatar Dmitri Naumov
Browse files

Merge pull request #982 from endJunction/MoveUpProcessBcs 

Move up process bcs
parents e4ba29d8 7ed9f2e2
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
ProcessLib/GroundwaterFlowProcess.h
+ 4
− 88
View file @ c8376b44
......@@ -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
......
ProcessLib/Process.h
+ 91
− 12
View file @ c8376b44
......@@ -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
......
ProcessLib/ProcessVariable.h
+ 3
− 2
View file @ c8376b44
......@@ -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)...)};
}
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment