diff --git a/src/amdis/CMakeLists.txt b/src/amdis/CMakeLists.txt
index 45e875a5359a865729daeeb9077b3abe1f06eed2..6ac8d9c5e3c396ced25e78011c0270ea5399619c 100644
--- a/src/amdis/CMakeLists.txt
+++ b/src/amdis/CMakeLists.txt
@@ -36,6 +36,8 @@ install(FILES
LinearAlgebra.hpp
LocalAssembler.hpp
LocalAssemblerBase.hpp
+ Marker.hpp
+ Marker.inc.hpp
Mesh.hpp
Operations.hpp
Operators.hpp
diff --git a/src/amdis/Marker.hpp b/src/amdis/Marker.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..75bcad474b7e6e016476b911744a9f6b47e9dfb8
--- /dev/null
+++ b/src/amdis/Marker.hpp
@@ -0,0 +1,333 @@
+#pragma once
+
+#include <string>
+
+#include <dune/grid/common/grid.hh>
+
+#include <amdis/AdaptInfo.hpp>
+#include <amdis/Flag.hpp>
+#include <amdis/Initfile.hpp>
+
+namespace AMDiS {
+
+
+ /**
+ * \ingroup Adaption
+ *
+ * \brief
+ * Base class for all scalar markers.
+ */
+ template <class Traits>
+ class Marker
+ {
+ protected:
+ using GlobalBasis = typename Traits::GlobalBasis;
+ using GridView = typename GlobalBasis::GridView;
+ using Grid = typename GridView::Grid;
+ using IndexSet = typename GridView::IndexSet;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Estimates = std::vector<double>;
+
+ public:
+ /// Constructor.
+ Marker() {}
+
+ /// Constructor.
+ Marker(std::string name, std::string component, Estimates const& est,
+ std::shared_ptr<Grid> const& grid)
+ : name_(name)
+ , component_(component)
+ , grid_(grid)
+ , est_(est)
+ , maximumMarking_(false)
+ {
+ Parameters::get(name_ + "->p", p_);
+ Parameters::get(name_ + "->info", info_);
+ Parameters::get(name_ + "->max refinement level", maxRefineLevel_);
+ Parameters::get(name_ + "->min refinement level", minRefineLevel_);
+ }
+
+ /// destructor
+ virtual ~Marker() {}
+
+ /// Marks element with newMark. If \ref maximumMarking_ is set, the element
+ /// is marked only if the new mark is bigger than the old one. The return
+ /// value specifies whether the element has been marked, or not.
+ void mark(Element const& elem, int newMark)
+ {
+ int oldMark = grid_->getMark(elem);
+
+ if (!maximumMarking_ || (newMark > oldMark)) {
+ bool marked = grid_->mark(newMark, elem);
+ if (marked) {
+ if (oldMark > 0) {
+ elMarkRefine_--;
+ } else if (oldMark < 0) {
+ elMarkCoarsen_--;
+ }
+
+ if (newMark > 0) {
+ elMarkRefine_++;
+ } else if (newMark < 0) {
+ elMarkCoarsen_++;
+ }
+ } else {
+ msg("Marking failed");
+ }
+ }
+ }
+
+ /// Can be used by sub classes. Called before traversal.
+ virtual void initMarking(AdaptInfo& adaptInfo);
+
+ /// Can be used by sub classes. Called after traversal.
+ virtual void finishMarking(AdaptInfo& adaptInfo);
+
+ /// Marks one element.
+ virtual void markElement(AdaptInfo& adaptInfo, Element const& elem);
+
+ /// Marking of the mesh.
+ virtual Flag markGrid(AdaptInfo& adaptInfo);
+
+ /// Sets \ref maximumMarking.
+ void setMaximumMarking(bool maxMark)
+ {
+ maximumMarking_ = maxMark;
+ }
+
+ int getElMarkRefine()
+ {
+ return elMarkRefine_;
+ }
+
+ int getElMarkCoarsen()
+ {
+ return elMarkCoarsen_;
+ }
+
+ /// Returns \ref name of the Marker
+ std::string getName() const
+ {
+ return name_;
+ }
+
+ /// Creates a scalar marker depending on the strategy set in parameters.
+ static std::shared_ptr<Marker<Traits> > createMarker(std::string name, std::string component, Estimates const& est, std::shared_ptr<Grid> const& grid);
+
+ protected:
+ /// Name of the scalar marker.
+ std::string name_;
+
+ /// Problem component to work on
+ std::string component_;
+
+ /// Pointer to the grid
+ std::shared_ptr<Grid> grid_;
+
+ /// Pointer to the local error estimates
+ Estimates est_;
+
+ /// estimation sum
+ double estSum_ = 0.0;
+
+ /// estmation maximum
+ double estMax_ = 0.0;
+
+ /// If true, elements are marked only if the new value is bigger than
+ /// the current marking.
+ bool maximumMarking_;
+
+ /// Lower limit for error estimation, from which an element is marked for
+ /// refinement
+ double markRLimit_;
+
+ /// Upper limit for error estimation, from which an element is marked for
+ /// coarsening
+ double markCLimit_;
+
+ /// power in estimator norm TODO: is this info necessary in marker?
+ double p_ = 2.0;
+
+ /// Info level.
+ int info_ = 10;
+
+ /// Counter for elements marked for refinement
+ int elMarkRefine_ = 0;
+
+ /// Counter for elements marked for coarsening
+ int elMarkCoarsen_ = 0;
+
+ /// Maximal level of all elements.
+ int maxRefineLevel_ = -1;
+
+ /// Minimal level of all elements.
+ int minRefineLevel_ = -1;
+
+ bool refineAllowed_ = true;
+
+ bool coarsenAllowed_ = false;
+ };
+
+
+ /**
+ * \ingroup Adaption
+ *
+ * \brief
+ * Global refinement.
+ */
+ template <class Traits>
+ class GRMarker
+ : public Marker<Traits>
+ {
+ using Super = Marker<Traits>;
+ using Grid = typename Super::Grid;
+ using Element = typename Super::Element;
+ using Estimates = typename Super::Estimates;
+
+ public:
+ /// Constructor.
+ GRMarker(std::string name, std::string component, Estimates const& est,
+ std::shared_ptr<Grid> const& grid)
+ : Marker<Traits>(std::move(name), std::move(component), est, grid)
+ {}
+
+ /// Implementation of Marker::markElement().
+ virtual void markElement(AdaptInfo& adaptInfo, Element const& elem)
+ {
+ if (this->refineAllowed_)
+ this->mark(elem, 1);
+ }
+ };
+
+
+ /**
+ * \ingroup Adaption
+ *
+ * \brief
+ * Maximum strategy.
+ */
+
+ template <class Traits>
+ class MSMarker
+ : public Marker<Traits>
+ {
+ using Super = Marker<Traits>;
+ using Grid = typename Super::Grid;
+ using Estimates = typename Super::Estimates;
+
+ public:
+ /// Constructor.
+ MSMarker(std::string name, std::string component, Estimates const& est,
+ std::shared_ptr<Grid> const& grid)
+ : Super{std::move(name), std::move(component), est, grid}
+ {
+ Parameters::get(this->name_ + "->MSGamma", msGamma_);
+ Parameters::get(this->name_ + "->MSGammaC", msGammaC_);
+ }
+
+ /// Implementation of Marker::initMarking().
+ void initMarking(AdaptInfo& adaptInfo);
+
+ protected:
+ /// Marking parameter.
+ double msGamma_ = 0.5;
+
+ /// Marking parameter.
+ double msGammaC_ = 0.1;
+ };
+
+
+ /**
+ * \ingroup Adaption
+ *
+ * \brief
+ * Equidistribution strategy
+ */
+
+ template <class Traits>
+ class ESMarker
+ : public Marker<Traits>
+ {
+ using Super = Marker<Traits>;
+ using Grid = typename Super::Grid;
+ using Estimates = typename Super::Estimates;
+
+ public:
+ /// Constructor.
+ ESMarker(std::string name, std::string component, Estimates const& est,
+ std::shared_ptr<Grid> const& grid)
+ : Super{std::move(name), std::move(component), est, grid}
+ {
+ Parameters::get(this->name_ + "->ESTheta", esTheta_);
+ Parameters::get(this->name_ + "->ESThetaC", esThetaC_);
+ }
+
+ /// Implementation of Marker::initMarking().
+ virtual void initMarking(AdaptInfo& adaptInfo);
+
+ protected:
+ /// Marking parameter.
+ double esTheta_ = 0.9;
+
+ /// Marking parameter.
+ double esThetaC_ = 0.2;
+ };
+
+
+ /**
+ * \ingroup Adaption
+ *
+ * \brief
+ * Guaranteed error reduction strategy
+ */
+
+ template <class Traits>
+ class GERSMarker
+ : public Marker<Traits>
+ {
+ using Super = Marker<Traits>;
+ using Grid = typename Super::Grid;
+ using Element = typename Super::Element;
+ using Estimates = typename Super::Estimates;
+
+ public:
+ /// Constructor.
+ GERSMarker(std::string name, std::string component, Estimates const& est,
+ std::shared_ptr<Grid> const& grid)
+ : Super{std::move(name), std::move(component), est, grid}
+ {
+ Parameters::get(this->name_ + "->GERSThetaStar", gersThetaStar_);
+ Parameters::get(this->name_ + "->GERSNu", gersNu_);
+ Parameters::get(this->name_ + "->GERSThetaC", gersThetaC_);
+ }
+
+ /// Implementation of Marker::markGrid().
+ virtual Flag markGrid(AdaptInfo& adaptInfo);
+
+ protected:
+ /// Refinement marking function.
+ void markElementForRefinement(AdaptInfo& adaptInfo, Element const& elem);
+
+ /// Coarsening marking function.
+ void markElementForCoarsening(AdaptInfo& adaptInfo, Element const& elem);
+
+ protected:
+ /// Marking parameter.
+ double gersSum_ = 0.0;
+
+ /// Marking parameter.
+ double oldErrSum_ = 0.0;
+
+ /// Marking parameter.
+ double gersThetaStar_ = 0.6;
+
+ /// Marking parameter.
+ double gersNu_ = 0.1;
+
+ /// Marking parameter.
+ double gersThetaC_ = 0.1;
+ };
+
+}
+
+#include "Marker.inc.hpp"
diff --git a/src/amdis/Marker.inc.hpp b/src/amdis/Marker.inc.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7f06e94d711ef5514d9ba1b33242da5394ca03c7
--- /dev/null
+++ b/src/amdis/Marker.inc.hpp
@@ -0,0 +1,240 @@
+#include <cmath>
+
+namespace AMDiS {
+
+template <class Traits>
+std::shared_ptr<Marker<Traits> > Marker<Traits>::createMarker(std::string name, std::string component, Estimates const& est, std::shared_ptr<Grid> const& grid)
+{
+ int strategy = 0;
+ Parameters::get(name + "->strategy", strategy);
+
+ std::shared_ptr<Marker<Traits> > marker;
+
+ switch (strategy) {
+ case 0: // no refinement/coarsening
+ break;
+ case 1:
+ msg("Creating global refinement (GR) marker");
+ marker = std::make_shared<GRMarker<Traits> >(name, component, est, grid);
+ break;
+ case 2:
+ msg("Creating maximum strategy (MS) marker");
+ marker = std::make_shared<MSMarker<Traits> >(name, component, est, grid);
+ break;
+ case 3:
+ msg("Creating equidistribution strategy (ES) marker");
+ marker = std::make_shared<ESMarker<Traits> >(name, component, est, grid);
+ break;
+ case 4:
+ msg("Creating guaranteed error reduction strategy (GERS) marker");
+ marker = std::make_shared<GERSMarker<Traits> >(name, component, est, grid);
+ break;
+ default:
+ error_exit("invalid strategy");
+ }
+
+ return marker;
+}
+
+
+template <class Traits>
+void Marker<Traits>::initMarking(AdaptInfo& adaptInfo)
+{
+ elMarkRefine_ = 0;
+ elMarkCoarsen_ = 0;
+ estSum_ = std::pow(adaptInfo.getEstSum(component_), p_);
+ estMax_ = adaptInfo.getEstMax(component_);
+ refineAllowed_ = adaptInfo.isRefinementAllowed(component_);
+ coarsenAllowed_ = adaptInfo.isCoarseningAllowed(component_);
+}
+
+
+template <class Traits>
+void Marker<Traits>::finishMarking(AdaptInfo& adaptInfo)
+{
+ msg(elMarkRefine_, " elements marked for refinement");
+ msg(elMarkCoarsen_, " elements marked for coarsening");
+}
+
+
+template <class Traits>
+void Marker<Traits>::markElement(AdaptInfo& adaptInfo, const Element& elem)
+{
+ const auto& index = grid_->leafIndexSet().index(elem);
+ double lError = est_[index];
+
+ if (lError > markRLimit_ && refineAllowed_
+ && (maxRefineLevel_ == -1 || elem.level() < maxRefineLevel_)) {
+ this->mark(elem, 1);
+ } else if (lError <= markCLimit_ && coarsenAllowed_
+ && (minRefineLevel_ == -1 || elem.level() > minRefineLevel_)) {
+ this->mark(elem, -1);
+ }
+}
+
+
+template <class Traits>
+Flag Marker<Traits>::markGrid(AdaptInfo& adaptInfo)
+{
+ if (!grid_)
+ error_exit("No grid!");
+
+ initMarking(adaptInfo);
+
+ if (!coarsenAllowed_ && !refineAllowed_)
+ return 0;
+
+ for (const auto& elem : Dune::elements(grid_->leafGridView()))
+ markElement(adaptInfo, elem);
+
+ finishMarking(adaptInfo);
+
+ Flag markFlag;
+ if (elMarkRefine_)
+ markFlag = 1;
+ if (elMarkCoarsen_)
+ markFlag |= 2;
+
+ return markFlag;
+}
+
+
+template <class Traits>
+void MSMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
+{
+ Super::initMarking(adaptInfo);
+
+ double msGammaP = std::pow(msGamma_, this->p_);
+ double msGammaCP = std::pow(msGammaC_, this->p_);
+
+ this->markRLimit_ = msGammaP * adaptInfo.getEstMax(this->component_);
+ this->markCLimit_ = msGammaCP * adaptInfo.getEstMax(this->component_);
+
+ msg("start max_est: ", adaptInfo.getEstMax(this->component_), ", mark_limits: ", this->markRLimit_, ", " , this->markCLimit_);
+}
+
+
+template <class Traits>
+void ESMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
+{
+ Super::initMarking(adaptInfo);
+
+ double esThetaP = std::pow(esTheta_, this->p_);
+ double esThetaCP = std::pow(esThetaC_, this->p_);
+ double epsP = std::pow(adaptInfo.getSpaceTolerance(this->component_), this->p_);
+
+ int nLeaves = (this->grid_->leafGridView()).size(0);
+/*#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ Parallel::mpi::globalAdd(nLeaves);
+#endif*/
+
+ this->markRLimit_ = esThetaP * epsP / nLeaves;
+ this->markCLimit_ = esThetaCP * epsP / nLeaves;
+
+ msg("start mark_limits: ", this->markRLimit_, ", " , this->markCLimit_, "; nt = ", nLeaves);
+}
+
+
+template <class Traits>
+Flag GERSMarker<Traits>::markGrid(AdaptInfo& adaptInfo)
+{
+ Super::initMarking(adaptInfo);
+
+ if (!this->coarsenAllowed_ && !this->refineAllowed_)
+ return 0;
+
+ gersSum_ = 0.0;
+
+ double LTheta = std::pow(1.0 - gersThetaStar_, this->p_);
+ double epsP = std::pow(adaptInfo.getSpaceTolerance(this->component_), this->p_);
+
+ if (this->estSum_ < oldErrSum_) {
+ double improv = this->estSum_ / oldErrSum_;
+ double wanted = 0.8 * epsP / this->estSum_;
+ double redfac = std::min((1.0 - wanted) / (1.0 - improv), 1.0);
+ redfac = std::max(redfac, 0.0);
+
+ if (redfac < 1.0) {
+ LTheta *= redfac;
+ msg("GERS: use extrapolated theta_star = ", std::pow(LTheta, 1.0 / this->p_));
+ }
+ }
+
+ oldErrSum_ = this->estSum_;
+ double GERSGamma = 1.0;
+
+ if (this->refineAllowed_) {
+ if (LTheta > 0) {
+ do {
+ gersSum_ = 0.0;
+ GERSGamma -= gersNu_;
+ this->markRLimit_ = GERSGamma * this->estMax_;
+
+ for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
+ markElementForRefinement(adaptInfo, elem);
+
+ } while((GERSGamma > 0) && (gersSum_ < LTheta * this->estSum_));
+ }
+
+ msg("GERS refinement with gamma = ", GERSGamma);
+ }
+
+ if (this->coarsenAllowed_) {
+ GERSGamma = 0.3;
+ LTheta = gersThetaC_ * epsP;
+
+ do {
+ gersSum_ = 0.0;
+ GERSGamma -= gersNu_;
+ this->markCLimit_ = GERSGamma * this->estMax_;
+
+ for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
+ markElementForCoarsening(adaptInfo, elem);
+
+ msg("coarse loop: gamma = ", GERSGamma, ", sum = ", gersSum_, ", limit = ", LTheta);
+ } while(gersSum_ > LTheta);
+
+ msg("GERS coarsening with gamma = ", GERSGamma);
+ }
+
+ Super::finishMarking(adaptInfo);
+
+ Flag markFlag;
+ if (this->elMarkRefine_)
+ markFlag = 1;
+ if (this->elMarkCoarsen_)
+ markFlag |= 2;
+
+ return markFlag;
+}
+
+
+template <class Traits>
+void GERSMarker<Traits>::markElementForRefinement(AdaptInfo& adaptInfo, const Element& elem)
+{
+ double lError = this->est_[this->grid_->leafIndexSet().index(elem)];
+
+ if (lError > this->markRLimit_) {
+ gersSum_ += lError;
+ this->mark(elem, 1);
+ }
+}
+
+
+template <class Traits>
+void GERSMarker<Traits>::markElementForCoarsening(AdaptInfo& adaptInfo, const Element& elem)
+{
+ double lError = this->est_[this->grid_->leafIndexSet().index(elem)];
+
+ if (this->grid_->getMark(elem) <= 0) {
+ if (lError <= this->markCLimit_) {
+ gersSum_ += lError;
+ this->mark(elem, -1);
+ } else {
+ this->mark(elem, 0);
+ }
+ }
+}
+
+
+} // end namespace AMDiS
diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp
index 28ebe21718a25a447d4c647f5642be3c8c65534f..9352089e83f137088cde0b6136cc5d3922baf007 100644
--- a/src/amdis/ProblemStat.hpp
+++ b/src/amdis/ProblemStat.hpp
@@ -19,6 +19,7 @@
#include <amdis/Flag.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/LinearAlgebra.hpp>
+#include <amdis/Marker.hpp>
#include <amdis/Mesh.hpp>
#include <amdis/ProblemStatBase.hpp>
#include <amdis/ProblemStatTraits.hpp>
@@ -33,6 +34,7 @@
#include <amdis/io/FileWriterInterface.hpp>
+#include <amdis/utility/TreeData.hpp>
#include <amdis/utility/TreePath.hpp>
namespace AMDiS
@@ -52,7 +54,7 @@ namespace AMDiS
using Element = typename GridView::template Codim<0>::Entity;
using WorldVector = typename Element::Geometry::GlobalCoordinate;
- /// Dimension of the mesh
+ /// Dimension of the grid
static constexpr int dim = Grid::dimension;
/// Dimension of the world
@@ -66,15 +68,15 @@ namespace AMDiS
public:
/**
* \brief Constructor. Takes the name of the problem that is used to
- * access values correpsonding to this püroblem in the parameter file.
+ * access values corresponding to this problem in the parameter file.
**/
explicit ProblemStat(std::string name)
: StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
, name_(std::move(name))
{}
- /// Constructor taking additionally a reference to a mesh that is used
- /// instead of the default created mesh, \ref ProblemStat
+ /// Constructor taking additionally a reference to a grid that is used
+ /// instead of the default created grid, \ref ProblemStat
ProblemStat(std::string name, Grid& grid)
: ProblemStat(std::move(name))
{
@@ -142,14 +144,12 @@ namespace AMDiS
bool createMatrixData = true,
bool storeMatrixData = false) override;
-
/// Implementation of \ref ProblemStatBase::buildAfterCoarse
virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Flag flag,
bool asmMatrix = true,
bool asmVector = true) override;
-
/// Writes output files.
void writeFiles(AdaptInfo& adaptInfo, bool force = false);
@@ -199,7 +199,7 @@ namespace AMDiS
/// Return a pointer to the grid, \ref grid
std::shared_ptr<Grid> getGrid() { return grid_; }
- /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
+ /// Set the grid. Stores pointer to passed reference and initializes feSpaces
/// matrices and vectors, as well as the file-writer.
void setGrid(Grid& grid)
{
@@ -211,12 +211,6 @@ namespace AMDiS
}
- /// Return the gridView of the leaf-level
- 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_; }
@@ -233,7 +227,7 @@ namespace AMDiS
{
gridName_ = "";
Parameters::get(name_ + "->mesh", gridName_);
- test_exit(!gridName_.empty(), "No mesh name specified for '", name_, "->mesh'!");
+ test_exit(!gridName_.empty(), "No grid name specified for '", name_, "->mesh'!");
grid_ = MeshCreator<Grid>::create(gridName_);
@@ -263,6 +257,16 @@ namespace AMDiS
systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat");
solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution");
rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
+
+ AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
+ {
+ std::string i = to_string(treePath);
+ estimates_[i].resize(globalBasis_->gridView().indexSet().size(0));
+ for (std::size_t j = 0; j < estimates_[i].size(); j++)
+ estimates_[i][j] = 0.0; // TODO: Remove when estimate() is implemented
+ });
+
+ rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
}
void createSolver()
@@ -276,6 +280,8 @@ namespace AMDiS
linearSolver_ = solverCreator->create(name_ + "->solver");
}
+ void createMarker();
+
void createFileWriter();
@@ -311,7 +317,7 @@ namespace AMDiS
virtual Flag coarsenMesh(AdaptInfo& adaptInfo) override { return 0; }
/// Implementation of \ref ProblemStatBase::markElements.
- virtual Flag markElements(AdaptInfo& adaptInfo) override { return 0; }
+ virtual Flag markElements(AdaptInfo& adaptInfo) override;
private:
@@ -321,7 +327,10 @@ namespace AMDiS
/// Grid of this problem.
std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems
- /// Name of the mesh
+ /// Number of grids
+ int nGrids = 1;
+
+ /// Name of the grid
std::string gridName_ = "none";
/// FE spaces of this problem.
@@ -331,6 +340,9 @@ namespace AMDiS
/// A FileWriter object
std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
+ /// Pointer to the adaptation markers
+ std::list<std::shared_ptr<Marker<Traits> > > marker_;
+
/// An object of the linearSolver Interface
std::shared_ptr<LinearSolverType> linearSolver_;
@@ -340,6 +352,10 @@ namespace AMDiS
/// A block-vector with the solution components
std::shared_ptr<SystemVector> solution_;
+ /// A vector with the local element error estimates
+ /// for each node in the basis tree, indexed by [to_string(treePath)][element index]
+ std::map<std::string, std::vector<double> > estimates_;
+
/// A block-vector (load-vector) corresponding to the right.hand side
/// of the equation, filled during assembling
std::shared_ptr<SystemVector> rhs_;
diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp
index 4499009bdc5a041383647220988da8a43dd196e8..263b5db3d35308019eb2b2eb6146853381c79815 100644
--- a/src/amdis/ProblemStat.inc.hpp
+++ b/src/amdis/ProblemStat.inc.hpp
@@ -13,6 +13,7 @@
#include <amdis/LocalAssembler.hpp>
#include <amdis/GridFunctionOperator.hpp>
#include <amdis/common/Loops.hpp>
+// #include <amdis/Estimator.hpp>
namespace AMDiS {
@@ -22,7 +23,7 @@ void ProblemStat<Traits>::initialize(
Self* adoptProblem,
Flag adoptFlag)
{
- // create grides
+ // create grids
if (grid_) {
warning("grid already created");
}
@@ -78,6 +79,7 @@ void ProblemStat<Traits>::initialize(
if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
solution_ = adoptProblem->solution_;
+ estimates_ = adoptProblem->estimates_;
rhs_ = adoptProblem->rhs_;
systemMatrix_ = adoptProblem->systemMatrix_;
}
@@ -101,6 +103,13 @@ void ProblemStat<Traits>::initialize(
warning("no solver created\n");
}
+ // create marker
+ if (initFlag.isSet(INIT_MARKER))
+ createMarker();
+
+ if (adoptProblem && adoptFlag.isSet(INIT_MARKER))
+ marker_ = adoptProblem->marker_;
+
// create file writer
if (initFlag.isSet(INIT_FILEWRITER))
@@ -110,6 +119,30 @@ void ProblemStat<Traits>::initialize(
}
+template <class Traits>
+void ProblemStat<Traits>::createMarker()
+{
+ AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
+ {
+ std::string componentName = name_ + "->marker[" + to_string(treePath) + "]";
+
+ if (!Parameters::get<std::string>(componentName + "->strategy"))
+ return;
+
+ std::string tp = to_string(treePath);
+ auto newMarker = Marker<Traits>::createMarker(componentName, tp, estimates_[tp], grid_);
+ if (newMarker)
+ marker_.push_back(std::move(newMarker));
+
+ // If there is more than one marker, and all components are defined
+ // on the same grid, the maximum marking has to be enabled.
+ // TODO: needs to be checked!
+ if (marker_.size() > 1)
+ marker_.back()->setMaximumMarking(true);
+ });
+}
+
+
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
{
@@ -280,6 +313,21 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
}
+template <class Traits>
+Flag ProblemStat<Traits>::markElements(AdaptInfo& adaptInfo)
+{
+ Dune::Timer t;
+
+ Flag markFlag = 0;
+ for (auto& currentMarker : marker_)
+ markFlag |= currentMarker->markGrid(adaptInfo);
+
+ msg("markElements needed ", t.elapsed(), " seconds");
+
+ return markFlag;
+}
+
+
template <class Traits>
void ProblemStat<Traits>::
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
@@ -287,7 +335,7 @@ buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool
Dune::Timer t;
Assembler<Traits> assembler(*globalBasis_, matrixOperators_, rhsOperators_, constraints_);
- auto gv = leafGridView();
+ auto gv = grid_->leafGridView();
globalBasis_->update(gv);
assembler.assemble(*systemMatrix_, *solution_, *rhs_, asmMatrix, asmVector);
@@ -305,4 +353,5 @@ writeFiles(AdaptInfo& adaptInfo, bool force)
msg("writeFiles needed ", t.elapsed(), " seconds");
}
+
} // end namespace AMDiS
diff --git a/src/amdis/utility/TreePath.hpp b/src/amdis/utility/TreePath.hpp
index 9d29ead4fa48046b9b6652b8a6ec233c24c39df7..4949a9281ee4025edd45171809451ffddd1a2603 100644
--- a/src/amdis/utility/TreePath.hpp
+++ b/src/amdis/utility/TreePath.hpp
@@ -8,8 +8,6 @@
#include <dune/typetree/treepath.hh>
#include <dune/typetree/typetraits.hh>
-#include <amdis/common/Mpl.hpp>
-
namespace AMDiS
{
struct RootTreePath {};
@@ -24,21 +22,26 @@ namespace AMDiS
{
template <class TP>
struct IsPreTreePath
- : any_of_t<std::is_same<TP, int>::value, std::is_same<TP,std::size_t>::value>
+ : std::is_integral<TP>
{};
template <int I>
- struct IsPreTreePath<int_t<I>>
+ struct IsPreTreePath<std::integral_constant<int, I>>
: std::true_type
{};
template <std::size_t I>
- struct IsPreTreePath<index_t<I>>
+ struct IsPreTreePath<std::integral_constant<std::size_t, I>>
+ : std::true_type
+ {};
+
+ template <int... I>
+ struct IsPreTreePath<std::integer_sequence<int, I...>>
: std::true_type
{};
template <std::size_t... I>
- struct IsPreTreePath<Indices<I...>>
+ struct IsPreTreePath<std::index_sequence<I...>>
: std::true_type
{};
@@ -87,16 +90,28 @@ namespace AMDiS
template <int I>
- auto makeTreePath(int_t<I>) { return Dune::TypeTree::hybridTreePath(index_<std::size_t(I)>); }
+ auto makeTreePath(std::integral_constant<int,I>)
+ {
+ return Dune::TypeTree::hybridTreePath(std::integral_constant<std::size_t, std::size_t(I)>{});
+ }
template <int... I>
- auto makeTreePath(Ints<I...>) { return Dune::TypeTree::hybridTreePath(index_<std::size_t(I)>...); }
+ auto makeTreePath(std::integer_sequence<int, I...>)
+ {
+ return Dune::TypeTree::hybridTreePath(std::integral_constant<std::size_t, std::size_t(I)>{}...);
+ }
template <std::size_t I>
- auto makeTreePath(index_t<I> _i) { return Dune::TypeTree::hybridTreePath(_i); }
+ auto makeTreePath(std::integral_constant<std::size_t,I> _i)
+ {
+ return Dune::TypeTree::hybridTreePath(_i);
+ }
template <std::size_t... I>
- auto makeTreePath(Indices<I...>) { return Dune::TypeTree::hybridTreePath(index_<I>...); }
+ auto makeTreePath(std::index_sequence<I...>)
+ {
+ return Dune::TypeTree::hybridTreePath(std::integral_constant<std::size_t, I>{}...);
+ }
template <class... T>
@@ -108,7 +123,7 @@ namespace AMDiS
template <std::size_t... I>
auto makeTreePath(Dune::TypeTree::TreePath<I...>)
{
- return Dune::TypeTree::hybridTreePath(index_<I>...);
+ return Dune::TypeTree::hybridTreePath(std::integral_constant<std::size_t, I>{}...);
}
template <class TP>
@@ -128,7 +143,7 @@ namespace AMDiS
void printTreePath(std::ostream& os, TreePath const& tp, std::index_sequence<I...>)
{
(void)std::initializer_list<int>{
- ((void)(os << treePathIndex(tp, index_<I>) << ","), 0)...
+ ((void)(os << treePathIndex(tp, std::integral_constant<std::size_t, I>{}) << ","), 0)...
};
}
}