Commit dd335839 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Merge branch 'issue/update_changes_from_feature_branch' into 'develop'

some changes from feature/block_matrix adopted

See merge request !19
parents 5a687fad 481178b5
Pipeline #1083 failed with stage
in 2 minutes and 32 seconds
...@@ -49,8 +49,8 @@ namespace AMDiS ...@@ -49,8 +49,8 @@ namespace AMDiS
using GlobalBasis = typename Traits::GlobalBasis; using GlobalBasis = typename Traits::GlobalBasis;
using GridView = typename GlobalBasis::GridView; using GridView = typename GlobalBasis::GridView;
using Grid = typename GridView::Grid; using Grid = typename GridView::Grid;
using Element = typename GridView::template Codim<0>::Entity; using Element = typename GridView::template Codim<0>::Entity;
using WorldVector = typename Element::Geometry::GlobalCoordinate;
/// Dimension of the mesh /// Dimension of the mesh
static constexpr int dim = Grid::dimension; static constexpr int dim = Grid::dimension;
...@@ -58,8 +58,6 @@ namespace AMDiS ...@@ -58,8 +58,6 @@ namespace AMDiS
/// Dimension of the world /// Dimension of the world
static constexpr int dow = Grid::dimensionworld; static constexpr int dow = Grid::dimensionworld;
using WorldVector = typename Element::Geometry::GlobalCoordinate;
using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>; using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
using SystemVector = DOFVector<Traits, double>; using SystemVector = DOFVector<Traits, double>;
...@@ -72,16 +70,25 @@ namespace AMDiS ...@@ -72,16 +70,25 @@ namespace AMDiS
**/ **/
explicit ProblemStat(std::string name) explicit ProblemStat(std::string name)
: StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this)) : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
, name(name) , name_(std::move(name))
{} {}
/// Constructor taking additionally a reference to a mesh that is used /// Constructor taking additionally a reference to a mesh that is used
/// instead of the default created mesh, \ref ProblemStat /// instead of the default created mesh, \ref ProblemStat
ProblemStat(std::string name, Grid& grid) ProblemStat(std::string name, Grid& grid)
: ProblemStat(name) : ProblemStat(std::move(name))
{ {
this->grid = Dune::stackobject_to_shared_ptr(grid); grid_ = Dune::stackobject_to_shared_ptr(grid);
Parameters::get(name + "->mesh", gridName); 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 ...@@ -147,57 +154,16 @@ namespace AMDiS
void writeFiles(AdaptInfo& adaptInfo, bool force = false); 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 public: // get-methods
/// Returns a pointer to system-matrix, \ref systemMatrix /// Returns a pointer to system-matrix, \ref systemMatrix_
auto getSystemMatrix() { return systemMatrix; } std::shared_ptr<SystemMatrix> getSystemMatrix() { return systemMatrix_; }
auto getSystemMatrix() const { 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 std::shared_ptr<SystemVector> getSolutionVector() const
{ {
return solution; return solution_;
} }
/// Return a mutable view to a solution component /// Return a mutable view to a solution component
...@@ -205,7 +171,7 @@ namespace AMDiS ...@@ -205,7 +171,7 @@ namespace AMDiS
auto getSolution(TreePath const& path = {}) auto getSolution(TreePath const& path = {})
{ {
auto&& tp = makeTreePath(path); auto&& tp = makeTreePath(path);
return makeDOFVectorView(*solution, tp); return makeDOFVectorView(*solution_, tp);
} }
/// Return a const view to a solution component /// Return a const view to a solution component
...@@ -213,31 +179,32 @@ namespace AMDiS ...@@ -213,31 +179,32 @@ namespace AMDiS
auto getSolution(TreePath const& path = {}) const auto getSolution(TreePath const& path = {}) const
{ {
auto&& tp = makeTreePath(path); auto&& tp = makeTreePath(path);
return makeDOFVectorView(*solution, tp); return makeDOFVectorView(*solution_, tp);
} }
/// Return a point to the rhs system-vector, \ref rhs /// Return a point to the rhs system-vector, \ref rhs
auto getRhsVector() { return rhs; } std::shared_ptr<SystemVector> getRhsVector() { return rhs_; }
auto getRhsVector() const { return rhs; } std::shared_ptr<SystemVector> getRhsVector() const { return rhs_; }
/// Return a pointer to the linear solver, \ref linearSolver /// 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) void setSolver(std::shared_ptr<LinearSolverType> const& solver)
{ {
linearSolver = solver; linearSolver_ = solver;
} }
/// Return a pointer to the grid, \ref grid /// 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 /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
/// matrices and vectors, as well as the file-writer. /// 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(); createGlobalBasis();
createMatricesAndVectors(); createMatricesAndVectors();
createFileWriter(); createFileWriter();
...@@ -245,109 +212,161 @@ namespace AMDiS ...@@ -245,109 +212,161 @@ namespace AMDiS
/// Return the gridView of the leaf-level /// 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 /// 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 /// Implementation of \ref ProblemStatBase::getName
virtual std::string getName() const override virtual std::string getName() const override
{ {
return name; return name_;
} }
protected: // initialization methods protected: // initialization methods
void createGrid() void createGrid()
{ {
gridName = ""; gridName_ = "";
Parameters::get(name + "->mesh", gridName); Parameters::get(name_ + "->mesh", gridName_);
test_exit(!gridName.empty(), "No mesh name specified for '", name, "->mesh'!"); 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("Create grid:");
msg("#elements = " , grid->size(0)); msg("#elements = " , grid_->size(0));
msg("#faces/edges = ", grid->size(1)); msg("#faces/edges = ", grid_->size(1));
msg("#vertices = " , grid->size(dim)); msg("#vertices = " , grid_->size(dim));
msg(""); msg("");
} }
void createGlobalBasis() void createGlobalBasis()
{ {
globalBasis = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(grid->leafGridView())); globalBasis_ = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(grid_->leafGridView()));
initGlobalBasis(*globalBasis); initGlobalBasis(*globalBasis_);
} }
void initGlobalBasis(GlobalBasis const& globalBasis) void initGlobalBasis(GlobalBasis const& globalBasis)
{ {
globalTree = std::make_shared<typename GlobalBasis::LocalView::Tree>(globalBasis.localView().tree()); localView_ = std::make_shared<typename GlobalBasis::LocalView>(globalBasis.localView());
matrixOperators.init(*globalTree); matrixOperators_.init(localView_->tree(), tag::store{});
rhsOperators.init(*globalTree); rhsOperators_.init(localView_->tree(), tag::store{});
constraints.init(*globalTree); constraints_.init(localView_->tree(), tag::store{});
} }
void createMatricesAndVectors() void createMatricesAndVectors()
{ {
systemMatrix = std::make_shared<SystemMatrix>(*globalBasis, *globalBasis, "mat"); systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat");
solution = std::make_shared<SystemVector>(*globalBasis, "solution"); solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution");
rhs = std::make_shared<SystemVector>(*globalBasis, "rhs"); rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
} }
void createSolver() void createSolver()
{ {
std::string solverName = "cg"; std::string solverName = "cg";
Parameters::get(name + "->solver->name", solverName); Parameters::get(name_ + "->solver->name", solverName);
auto solverCreator 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(); 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: private:
/// Name of this problem. /// Name of this problem.
std::string name; std::string name_;
/// Grid of this problem. /// 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 /// Name of the mesh
std::string gridName = "none"; std::string gridName_ = "none";
/// FE spaces of this problem. /// FE spaces of this problem.
std::shared_ptr<GlobalBasis> globalBasis; std::shared_ptr<GlobalBasis> globalBasis_;
std::shared_ptr<typename GlobalBasis::LocalView::Tree> globalTree; std::shared_ptr<typename GlobalBasis::LocalView> localView_;
/// A FileWriter object /// A FileWriter object
std::list<std::shared_ptr<FileWriterInterface>> filewriter; std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
/// An object of the linearSolver Interface /// An object of the linearSolver Interface
std::shared_ptr<LinearSolverType> linearSolver; std::shared_ptr<LinearSolverType> linearSolver_;
/// A block-matrix that is filled during assembling /// 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 /// 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 /// A block-vector (load-vector) corresponding to the right.hand side
/// of the equation, filled during assembling /// of the equation, filled during assembling
std::shared_ptr<SystemVector> rhs; std::shared_ptr<SystemVector> rhs_;
private: // some internal data-structures private: // some internal data-structures
MatrixOperators<GlobalBasis> matrixOperators; MatrixOperators<GlobalBasis> matrixOperators_;
VectorOperators<GlobalBasis> rhsOperators; VectorOperators<GlobalBasis> rhsOperators_;
Constraints<GlobalBasis> constraints_;
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 #ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
......
...@@ -23,7 +23,7 @@ void ProblemStat<Traits>::initialize( ...@@ -23,7 +23,7 @@ void ProblemStat<Traits>::initialize(
Flag adoptFlag) Flag adoptFlag)
{ {
// create grides // create grides
if (grid) { if (grid_) {
warning("grid already created"); warning("grid already created");
} }
else { else {
...@@ -37,22 +37,22 @@ void ProblemStat<Traits>::initialize( ...@@ -37,22 +37,22 @@ void ProblemStat<Traits>::initialize(
(adoptFlag.isSet(INIT_MESH) || (adoptFlag.isSet(INIT_MESH) ||
adoptFlag.isSet(INIT_SYSTEM) || adoptFlag.isSet(INIT_SYSTEM) ||
adoptFlag.isSet(INIT_FE_SPACE))) { adoptFlag.isSet(INIT_FE_SPACE))) {
grid = adoptProblem->getGrid(); grid_ = adoptProblem->grid_;
} }
} }
if (!grid) if (!grid_)
warning("no grid created"); warning("no grid created");
int globalRefinements = 0; int globalRefinements = 0;
Parameters::get(gridName + "->global refinements", globalRefinements); Parameters::get(gridName_ + "->global refinements", globalRefinements);
if (globalRefinements > 0) { if (globalRefinements > 0) {
grid->globalRefine(globalRefinements); grid_->globalRefine(globalRefinements);
} }
// create fespace // create fespace
if (globalBasis) { if (globalBasis_) {
warning("globalBasis already created"); warning("globalBasis already created");
} }
else { else {
...@@ -63,11 +63,12 @@ void ProblemStat<Traits>::initialize( ...@@ -63,11 +63,12 @@ void ProblemStat<Traits>::initialize(
if (adoptProblem && if (adoptProblem &&
(adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) { (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"); warning("no globalBasis created\n");
...@@ -76,14 +77,14 @@ void ProblemStat<Traits>::initialize( ...@@ -76,14 +77,14 @@ void ProblemStat<Traits>::initialize(
createMatricesAndVectors(); createMatricesAndVectors();
if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) { if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
solution = adoptProblem->solution; solution_ = adoptProblem->solution_;
rhs = adoptProblem->rhs; rhs_ = adoptProblem->rhs_;
systemMatrix = adoptProblem->systemMatrix; systemMatrix_ = adoptProblem->systemMatrix_;
} }
// create solver // create solver
if (linearSolver) { if (linearSolver_) {
warning("solver already created\n"); warning("solver already created\n");
} }
else { else {
...@@ -91,12 +92,12 @@ void ProblemStat<Traits>::initialize( ...@@ -91,12 +92,12 @@ void ProblemStat<Traits>::initialize(
createSolver(); createSolver();
if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) { if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
test_exit(!linearSolver, "solver already created\n"); test_exit(!linearSolver_, "solver already created\n");
linearSolver = adoptProblem->linearSolver; linearSolver_ = adoptProblem->linearSolver_;
} }
} }
if (!linearSolver) { if (!linearSolver_) {
warning("no solver created\n"); warning("no solver created\n");
} }
...@@ -105,23 +106,22 @@ void ProblemStat<Traits>::initialize( ...@@ -105,23 +106,22 @@ void ProblemStat<Traits>::initialize(
if (initFlag.isSet(INIT_FILEWRITER)) if (initFlag.isSet(INIT_FILEWRITER))
createFileWriter(); createFileWriter();
solution->compress(); solution_->compress();
} }
template <class Traits> template <class Traits>
void ProblemStat<Traits>::createFileWriter() 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")) if (!Parameters::get<std::string>(componentName + "->filename"))
return; return;
auto writer = makeFileWriterPtr(componentName, this->getSolution(treePath)); auto writer = makeFileWriterPtr(componentName, this->getSolution(treePath));
filewriter.push_back(std::move(writer));