Skip to content
GitLab
Commit 481178b5 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

some changes from feature/block_matrix adopted

parent 5a687fad
Hide whitespace changes
Inline Side-by-side
src/amdis/ProblemStat.hpp
View file @ 481178b5
......@@ -49,8 +49,8 @@ namespace AMDiS
using GlobalBasis = typename Traits::GlobalBasis;
using GridView = typename GlobalBasis::GridView;
using Grid = typename GridView::Grid;
using Element = typename GridView::template Codim<0>::Entity;
using WorldVector = typename Element::Geometry::GlobalCoordinate;
/// Dimension of the mesh
static constexpr int dim = Grid::dimension;
......@@ -58,8 +58,6 @@ namespace AMDiS
/// Dimension of the world
static constexpr int dow = Grid::dimensionworld;
using WorldVector = typename Element::Geometry::GlobalCoordinate;
using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
using SystemVector = DOFVector<Traits, double>;
......@@ -72,16 +70,25 @@ namespace AMDiS
**/
explicit ProblemStat(std::string name)
: StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
, name(name)
, name_(std::move(name))
{}
/// Constructor taking additionally a reference to a mesh that is used
/// instead of the default created mesh, \ref ProblemStat
ProblemStat(std::string name, Grid& grid)
: ProblemStat(name)
: ProblemStat(std::move(name))
{
this->grid = Dune::stackobject_to_shared_ptr(grid);
Parameters::get(name + "->mesh", gridName);
grid_ = Dune::stackobject_to_shared_ptr(grid);
Parameters::get(name_ + "->mesh", gridName_);
}
/// \brief Constructor taking a grid reference and a basis reference.
/// Stores pointers to both.
ProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
: ProblemStat(std::move(name), grid)
{
globalBasis_ = Dune::stackobject_to_shared_ptr(globalBasis);
initGlobalBasis(*globalBasis_);
}
......@@ -147,57 +154,16 @@ namespace AMDiS
void writeFiles(AdaptInfo& adaptInfo, bool force = false);
public: // implementation of iteration interface methods
/**
* \brief Determines the execution order of the single adaption steps.
*
* If adapt is true, mesh adaption will be performed. This allows to avoid
* mesh adaption, e.g. in timestep adaption loops of timestep adaptive
* strategies.
*
* Implementation of \ref StandardProblemIteration::oneIteration.
**/
virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
{
// for (std::size_t i = 0; i < getNumComponents(); ++i)
// if (adaptInfo.spaceToleranceReached(i))
// adaptInfo.allowRefinement(false, i);
// else
// adaptInfo.allowRefinement(true, i);
return StandardProblemIteration::oneIteration(adaptInfo, toDo);
}
/// Implementation of \ref ProblemStatBase::buildBeforeRefine.
virtual void buildBeforeRefine(AdaptInfo&, Flag) override { /* Does nothing. */ }
/// Implementation of \ref ProblemStatBase::buildBeforeCoarsen.
virtual void buildBeforeCoarsen(AdaptInfo&, Flag) override { /* Does nothing. */ }
/// Implementation of \ref ProblemStatBase::estimate.
virtual void estimate(AdaptInfo& adaptInfo) override { /* Does nothing. */ }
/// Implementation of \ref ProblemStatBase::refineMesh.
virtual Flag refineMesh(AdaptInfo& adaptInfo) override { return 0; }
/// Implementation of \ref ProblemStatBase::coarsenMesh.
virtual Flag coarsenMesh(AdaptInfo& adaptInfo) override { return 0; }
/// Implementation of \ref ProblemStatBase::markElements.
virtual Flag markElements(AdaptInfo& adaptInfo) override { return 0; }
public: // get-methods
/// Returns a pointer to system-matrix, \ref systemMatrix
auto getSystemMatrix() { return systemMatrix; }
auto getSystemMatrix() const { return systemMatrix; }
/// Returns a pointer to system-matrix, \ref systemMatrix_
std::shared_ptr<SystemMatrix> getSystemMatrix() { return systemMatrix_; }
std::shared_ptr<SystemMatrix> getSystemMatrix() const { return systemMatrix_; }
/// Returns a pointer to the solution vector, \ref solution_
std::shared_ptr<SystemVector> getSolutionVector() const
{
return solution;
return solution_;
}
/// Return a mutable view to a solution component
......@@ -205,7 +171,7 @@ namespace AMDiS
auto getSolution(TreePath const& path = {})
{
auto&& tp = makeTreePath(path);
return makeDOFVectorView(*solution, tp);
return makeDOFVectorView(*solution_, tp);
}
/// Return a const view to a solution component
......@@ -213,31 +179,32 @@ namespace AMDiS
auto getSolution(TreePath const& path = {}) const
{
auto&& tp = makeTreePath(path);
return makeDOFVectorView(*solution, tp);
return makeDOFVectorView(*solution_, tp);
}
/// Return a point to the rhs system-vector, \ref rhs
auto getRhsVector() { return rhs; }
auto getRhsVector() const { return rhs; }
std::shared_ptr<SystemVector> getRhsVector() { return rhs_; }
std::shared_ptr<SystemVector> getRhsVector() const { return rhs_; }
/// Return a pointer to the linear solver, \ref linearSolver
std::shared_ptr<LinearSolverType> getSolver() { return linearSolver; }
std::shared_ptr<LinearSolverType> getSolver() { return linearSolver_; }
void setSolver(std::shared_ptr<LinearSolverType> const& solver)
{
linearSolver = solver;
linearSolver_ = solver;
}
/// Return a pointer to the grid, \ref grid
std::shared_ptr<Grid> getGrid() { return grid; }
std::shared_ptr<Grid> getGrid() { return grid_; }
/// Set the mesh. Stores pointer to passed reference and initializes feSpaces
/// matrices and vectors, as well as the file-writer.
void setGrid(Grid& grid_)
void setGrid(Grid& grid)
{
grid = Dune::stackobject_to_shared_ptr(grid_);
grid_ = Dune::stackobject_to_shared_ptr(grid);
createGlobalBasis();
createMatricesAndVectors();
createFileWriter();
......@@ -245,109 +212,161 @@ namespace AMDiS
/// Return the gridView of the leaf-level
auto const& gridView() { return globalBasis->gridView(); }
auto leafGridView() { return grid_->leafGridView(); }
/// Return the gridView of levle `level`
auto levelGridView(int level) { return grid_->levelGridView(level); }
/// Return the \ref feSpaces
std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis; }
std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis_; }
/// Implementation of \ref ProblemStatBase::getName
virtual std::string getName() const override
{
return name;
return name_;
}
protected: // initialization methods
void createGrid()
{
gridName = "";
Parameters::get(name + "->mesh", gridName);
test_exit(!gridName.empty(), "No mesh name specified for '", name, "->mesh'!");
gridName_ = "";
Parameters::get(name_ + "->mesh", gridName_);
test_exit(!gridName_.empty(), "No mesh name specified for '", name_, "->mesh'!");
grid = MeshCreator<Grid>::create(gridName);
grid_ = MeshCreator<Grid>::create(gridName_);
msg("Create grid:");
msg("#elements = " , grid->size(0));
msg("#faces/edges = ", grid->size(1));
msg("#vertices = " , grid->size(dim));
msg("#elements = " , grid_->size(0));
msg("#faces/edges = ", grid_->size(1));
msg("#vertices = " , grid_->size(dim));
msg("");
}
void createGlobalBasis()
{
globalBasis = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(grid->leafGridView()));
initGlobalBasis(*globalBasis);
globalBasis_ = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(grid_->leafGridView()));
initGlobalBasis(*globalBasis_);
}
void initGlobalBasis(GlobalBasis const& globalBasis)
{
globalTree = std::make_shared<typename GlobalBasis::LocalView::Tree>(globalBasis.localView().tree());
matrixOperators.init(*globalTree);
rhsOperators.init(*globalTree);
constraints.init(*globalTree);
localView_ = std::make_shared<typename GlobalBasis::LocalView>(globalBasis.localView());
matrixOperators_.init(localView_->tree(), tag::store{});
rhsOperators_.init(localView_->tree(), tag::store{});
constraints_.init(localView_->tree(), tag::store{});
}
void createMatricesAndVectors()
{
systemMatrix = std::make_shared<SystemMatrix>(*globalBasis, *globalBasis, "mat");
solution = std::make_shared<SystemVector>(*globalBasis, "solution");
rhs = std::make_shared<SystemVector>(*globalBasis, "rhs");
systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat");
solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution");
rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
}
void createSolver()
{
std::string solverName = "cg";
Parameters::get(name + "->solver->name", solverName);
Parameters::get(name_ + "->solver->name", solverName);
auto solverCreator
= named(CreatorMap<LinearSolverType>::getCreator(solverName, name + "->solver->name"));
= named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));
linearSolver = solverCreator->create(name + "->solver");
linearSolver_ = solverCreator->create(name_ + "->solver");
}
void createFileWriter();
public: // implementation of iteration interface methods
/**
* \brief Determines the execution order of the single adaption steps.
*
* If adapt is true, mesh adaption will be performed. This allows to avoid
* mesh adaption, e.g. in timestep adaption loops of timestep adaptive
* strategies.
*
* Implementation of \ref StandardProblemIteration::oneIteration.
**/
virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
{
return StandardProblemIteration::oneIteration(adaptInfo, toDo);
}
/// Implementation of \ref ProblemStatBase::buildBeforeRefine.
virtual void buildBeforeRefine(AdaptInfo&, Flag) override { /* Does nothing. */ }
/// Implementation of \ref ProblemStatBase::buildBeforeCoarsen.
virtual void buildBeforeCoarsen(AdaptInfo&, Flag) override { /* Does nothing. */ }
/// Implementation of \ref ProblemStatBase::estimate.
virtual void estimate(AdaptInfo& adaptInfo) override { /* Does nothing. */ }
/// Implementation of \ref ProblemStatBase::refineMesh.
virtual Flag refineMesh(AdaptInfo& adaptInfo) override { return 0; }
/// Implementation of \ref ProblemStatBase::coarsenMesh.
virtual Flag coarsenMesh(AdaptInfo& adaptInfo) override { return 0; }
/// Implementation of \ref ProblemStatBase::markElements.
virtual Flag markElements(AdaptInfo& adaptInfo) override { return 0; }
private:
/// Name of this problem.
std::string name;
std::string name_;
/// Grid of this problem.
std::shared_ptr<Grid> grid; // TODO: generalize to multi-mesh problems
std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems
/// Name of the mesh
std::string gridName = "none";
std::string gridName_ = "none";
/// FE spaces of this problem.
std::shared_ptr<GlobalBasis> globalBasis;
std::shared_ptr<typename GlobalBasis::LocalView::Tree> globalTree;
std::shared_ptr<GlobalBasis> globalBasis_;
std::shared_ptr<typename GlobalBasis::LocalView> localView_;
/// A FileWriter object
std::list<std::shared_ptr<FileWriterInterface>> filewriter;
std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
/// An object of the linearSolver Interface
std::shared_ptr<LinearSolverType> linearSolver;
std::shared_ptr<LinearSolverType> linearSolver_;
/// A block-matrix that is filled during assembling
std::shared_ptr<SystemMatrix> systemMatrix;
std::shared_ptr<SystemMatrix> systemMatrix_;
/// A block-vector with the solution components
std::shared_ptr<SystemVector> solution;
std::shared_ptr<SystemVector> solution_;
/// A block-vector (load-vector) corresponding to the right.hand side
/// of the equation, filled during assembling
std::shared_ptr<SystemVector> rhs;
std::shared_ptr<SystemVector> rhs_;
private: // some internal data-structures
MatrixOperators<GlobalBasis> matrixOperators;
VectorOperators<GlobalBasis> rhsOperators;
Constraints<GlobalBasis> constraints;
MatrixOperators<GlobalBasis> matrixOperators_;
VectorOperators<GlobalBasis> rhsOperators_;
Constraints<GlobalBasis> constraints_;
};
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
// Deduction rule
template <class Grid, class GlobalBasis>
ProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
-> ProblemStat<DefaultProblemTraits<GlobalBasis>>;
#endif
// Generator for ProblemStat with given Grid and GlobalBasis
template <class Grid, class GlobalBasis>
ProblemStat<DefaultProblemTraits<GlobalBasis>>
makeProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
{
return {std::move(name), grid, globalBasis};
}
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
......
src/amdis/ProblemStat.inc.hpp
View file @ 481178b5
......@@ -23,7 +23,7 @@ void ProblemStat<Traits>::initialize(
Flag adoptFlag)
{
// create grides
if (grid) {
if (grid_) {
warning("grid already created");
}
else {
......@@ -37,22 +37,22 @@ void ProblemStat<Traits>::initialize(
(adoptFlag.isSet(INIT_MESH) ||
adoptFlag.isSet(INIT_SYSTEM) ||
adoptFlag.isSet(INIT_FE_SPACE))) {
grid = adoptProblem->getGrid();
grid_ = adoptProblem->grid_;
}
}
if (!grid)
if (!grid_)
warning("no grid created");
int globalRefinements = 0;
Parameters::get(gridName + "->global refinements", globalRefinements);
Parameters::get(gridName_ + "->global refinements", globalRefinements);
if (globalRefinements > 0) {
grid->globalRefine(globalRefinements);
grid_->globalRefine(globalRefinements);
}
// create fespace
if (globalBasis) {
if (globalBasis_) {
warning("globalBasis already created");
}
else {
......@@ -63,11 +63,12 @@ void ProblemStat<Traits>::initialize(
if (adoptProblem &&
(adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
globalBasis = adoptProblem->globalBasis;
globalBasis_ = adoptProblem->globalBasis_;
initGlobalBasis(*globalBasis_);
}
}
if (!globalBasis)
if (!globalBasis_)
warning("no globalBasis created\n");
......@@ -76,14 +77,14 @@ void ProblemStat<Traits>::initialize(
createMatricesAndVectors();
if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
solution = adoptProblem->solution;
rhs = adoptProblem->rhs;
systemMatrix = adoptProblem->systemMatrix;
solution_ = adoptProblem->solution_;
rhs_ = adoptProblem->rhs_;
systemMatrix_ = adoptProblem->systemMatrix_;
}
// create solver
if (linearSolver) {
if (linearSolver_) {
warning("solver already created\n");
}
else {
......@@ -91,12 +92,12 @@ void ProblemStat<Traits>::initialize(
createSolver();
if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
test_exit(!linearSolver, "solver already created\n");
linearSolver = adoptProblem->linearSolver;
test_exit(!linearSolver_, "solver already created\n");
linearSolver_ = adoptProblem->linearSolver_;
}
}
if (!linearSolver) {
if (!linearSolver_) {
warning("no solver created\n");
}
......@@ -105,23 +106,22 @@ void ProblemStat<Traits>::initialize(
if (initFlag.isSet(INIT_FILEWRITER))
createFileWriter();
solution->compress();
solution_->compress();
}
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
{
auto localView = globalBasis->localView();
AMDiS::forEachNode_(localView.tree(), [&,this](auto const& node, auto treePath)
AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
{
std::string componentName = name + "->output[" + to_string(treePath) + "]";
std::string componentName = name_ + "->output[" + to_string(treePath) + "]";
if (!Parameters::get<std::string>(componentName + "->filename"))
return;
auto writer = makeFileWriterPtr(componentName, this->getSolution(treePath));
filewriter.push_back(std::move(writer));
filewriter_.push_back(std::move(writer));
});
}
......@@ -140,14 +140,14 @@ void ProblemStat<Traits>::addMatrixOperator(
static_assert( Concepts::PreTreePath<ColTreePath>,
"col must be a valid treepath, or an integer/index-constant");
auto i = child(globalBasis->localView().tree(), makeTreePath(row));
auto j = child(globalBasis->localView().tree(), makeTreePath(col));
auto i = child(localView_->tree(), makeTreePath(row));
auto j = child(localView_->tree(), makeTreePath(col));
auto op = makeGridOperator(preOp, globalBasis->gridView());
auto op = makeGridOperator(preOp, globalBasis_->gridView());
auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i, j);
matrixOperators[i][j].element.push_back({localAssembler, factor, estFactor});
matrixOperators[i][j].changing = true;
matrixOperators_[i][j].element.push_back({localAssembler, factor, estFactor});
matrixOperators_[i][j].changing = true;
}
......@@ -164,15 +164,15 @@ void ProblemStat<Traits>::addMatrixOperator(
static_assert( Concepts::PreTreePath<ColTreePath>,
"col must be a valid treepath, or an integer/index-constant");
auto i = child(globalBasis->localView().tree(), makeTreePath(row));
auto j = child(globalBasis->localView().tree(), makeTreePath(col));
auto i = child(localView_->tree(), makeTreePath(row));
auto j = child(localView_->tree(), makeTreePath(col));
using Intersection = typename GridView::Intersection;
auto op = makeGridOperator(preOp, globalBasis->gridView());
auto op = makeGridOperator(preOp, globalBasis_->gridView());
auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i, j);
matrixOperators[i][j].boundary.push_back({localAssembler, factor, estFactor, b});
matrixOperators[i][j].changing = true;
matrixOperators_[i][j].boundary.push_back({localAssembler, factor, estFactor, b});
matrixOperators_[i][j].changing = true;
}
......@@ -186,13 +186,13 @@ void ProblemStat<Traits>::addVectorOperator(
static_assert( Concepts::PreTreePath<TreePath>,
"path must be a valid treepath, or an integer/index-constant");
auto i = child(globalBasis->localView().tree(), makeTreePath(path));
auto i = child(localView_->tree(), makeTreePath(path));
auto op = makeGridOperator(preOp, globalBasis->gridView());
auto op = makeGridOperator(preOp, globalBasis_->gridView());
auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i);
rhsOperators[i].element.push_back({localAssembler, factor, estFactor});
rhsOperators[i].changing = true;
rhsOperators_[i].element.push_back({localAssembler, factor, estFactor});
rhsOperators_[i].changing = true;
}
......@@ -207,14 +207,14 @@ void ProblemStat<Traits>::addVectorOperator(
static_assert( Concepts::PreTreePath<TreePath>,
"path must be a valid treepath, or an integer/index-constant");
auto i = child(globalBasis->localView().tree(), makeTreePath(path));
auto i = child(localView_->tree(), makeTreePath(path));
using Intersection = typename GridView::Intersection;
auto op = makeGridOperator(preOp, globalBasis->gridView());
auto op = makeGridOperator(preOp, globalBasis_->gridView());
auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i);
rhsOperators[i].boundary.push_back({localAssembler, factor, estFactor, b});
rhsOperators[i].changing = true;
rhsOperators_[i].boundary.push_back({localAssembler, factor, estFactor, b});
rhsOperators_[i].changing = true;
}
......@@ -227,15 +227,15 @@ addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Val
static_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
"Function passed to addDirichletBC for `predicate` does not model the Functor<bool(WorldVector)> concept");
auto i = child(globalBasis->localView().tree(), makeTreePath(row));
auto j = child(globalBasis->localView().tree(), makeTreePath(col));
auto i = child(localView_->tree(), makeTreePath(row));
auto j = child(localView_->tree(), makeTreePath(col));
auto valueGridFct = makeGridFunction(values, globalBasis->gridView());
auto valueGridFct = makeGridFunction(values, globalBasis_->gridView());
using Range = RangeType_t<decltype(i)>;
using BcType = DirichletBC<WorldVector,Range>;
auto bc = std::make_shared<BcType>(predicate, valueGridFct);
constraints[i][j].push_back(bc);
constraints_[i][j].push_back(bc);
// TODO: make DirichletBC an abstract class and add specialization with gridfunction type
}
......@@ -246,14 +246,13 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
Dune::Timer t;
SolverInfo solverInfo(name + "->solver");
SolverInfo solverInfo(name_ + "->solver");
solverInfo.setCreateMatrixData(createMatrixData);
solverInfo.setStoreMatrixData(storeMatrixData);
solution->compress();
solution_->compress();
linearSolver->solve(systemMatrix->getMatrix(),
solution->getVector(), rhs->getVector(),
linearSolver_->solve(systemMatrix_->getMatrix(), solution_->getVector(), rhs_->getVector(),
solverInfo);
if (solverInfo.getInfo() > 0) {
......@@ -287,8 +286,10 @@ buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool
{
Dune::Timer t;
Assembler<Traits> assembler(*globalBasis, matrixOperators, rhsOperators, constraints);
assembler.assemble(*systemMatrix, *solution, *rhs, asmMatrix, asmVector);
Assembler<Traits> assembler(*globalBasis_, matrixOperators_, rhsOperators_, constraints_);
auto gv = leafGridView();
globalBasis_->update(gv);
assembler.assemble(*systemMatrix_, *solution_, *rhs_, asmMatrix, asmVector);
msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
}
......@@ -299,7 +300,7 @@ void ProblemStat<Traits>::
writeFiles(AdaptInfo& adaptInfo, bool force)
{
Dune::Timer t;
for (auto writer : filewriter)
for (auto writer : filewriter_)
writer->writeFiles(adaptInfo, force);
msg("writeFiles needed ", t.elapsed(), " seconds");
}
......
src/amdis/gridfunctions/DOFVectorView.inc.hpp
View file @ 481178b5
......@@ -13,7 +13,7 @@ LocalFunction::operator()(LocalDomain const& x) const
auto&& coefficients = *globalFunction_->dofVector_;
auto&& nodeToRangeEntry = globalFunction_->nodeToRangeEntry_;
forEachLeafNode(*subTree_, [&,this](auto const& node, auto)
AMDiS::forEachLeafNode_(*subTree_, [&,this](auto const& node, auto)