From 9379904a1f61e2cb2b5262b8d0ab5ccaf7ebf956 Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Thu, 24 May 2018 12:42:19 +0200
Subject: [PATCH] rebased on develop and cleaned up
---
src/amdis/Marker.hpp | 234 ++++++++++-----------
src/amdis/Marker.inc.hpp | 370 ++++++++++++++++-----------------
src/amdis/ProblemStat.hpp | 70 +------
src/amdis/ProblemStat.inc.hpp | 109 ++--------
src/amdis/utility/TreePath.hpp | 39 ++--
5 files changed, 334 insertions(+), 488 deletions(-)
diff --git a/src/amdis/Marker.hpp b/src/amdis/Marker.hpp
index 628c6b42..75bcad47 100644
--- a/src/amdis/Marker.hpp
+++ b/src/amdis/Marker.hpp
@@ -1,5 +1,5 @@
#pragma once
-// TODO: Cleanup of copied comments
+
#include <string>
#include <dune/grid/common/grid.hh>
@@ -10,7 +10,7 @@
namespace AMDiS {
-
+
/**
* \ingroup Adaption
*
@@ -20,60 +20,56 @@ namespace AMDiS {
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 EstType = std::vector<double>;
+ using Estimates = std::vector<double>;
public:
/// Constructor.
Marker() {}
/// Constructor.
- Marker(std::string name_, std::string component_, const EstType& est_, std::shared_ptr<Grid> grid_)
- : name(name_),
- component(component_),
- grid(grid_),
- est(est_),
- maximumMarking(false),
- p(2),
- info(10),
- maxRefineLevel(-1),
- minRefineLevel(-1),
- refineAllowed(true),
- coarsenAllowed(false)
+ 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);
+ 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
+ /// 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.
- inline void mark(const Element& elem, int newMark)
+ void mark(Element const& elem, int newMark)
{
- int oldMark = grid->getMark(elem);
+ int oldMark = grid_->getMark(elem);
- if (!maximumMarking || (newMark > oldMark)) {
- bool (marked) = grid->mark(newMark, elem);
+ if (!maximumMarking_ || (newMark > oldMark)) {
+ bool marked = grid_->mark(newMark, elem);
if (marked) {
if (oldMark > 0) {
- elMarkRefine--;
+ elMarkRefine_--;
} else if (oldMark < 0) {
- elMarkCoarsen--;
+ elMarkCoarsen_--;
}
if (newMark > 0) {
- elMarkRefine++;
+ elMarkRefine_++;
} else if (newMark < 0) {
- elMarkCoarsen++;
+ elMarkCoarsen_++;
}
} else {
msg("Marking failed");
@@ -88,88 +84,88 @@ namespace AMDiS {
virtual void finishMarking(AdaptInfo& adaptInfo);
/// Marks one element.
- virtual void markElement(AdaptInfo& adaptInfo, const Element& elem);
+ virtual void markElement(AdaptInfo& adaptInfo, Element const& elem);
/// Marking of the mesh.
virtual Flag markGrid(AdaptInfo& adaptInfo);
/// Sets \ref maximumMarking.
- inline void setMaximumMarking(bool maxMark)
+ void setMaximumMarking(bool maxMark)
{
- maximumMarking = maxMark;
+ maximumMarking_ = maxMark;
}
- inline int getElMarkRefine()
+ int getElMarkRefine()
{
- return elMarkRefine;
+ return elMarkRefine_;
}
- inline int getElMarkCoarsen()
+ int getElMarkCoarsen()
{
- return elMarkCoarsen;
+ return elMarkCoarsen_;
}
/// Returns \ref name of the Marker
- inline std::string getName() const
+ std::string getName() const
{
- return name;
+ 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_, const EstType& est_, std::shared_ptr<Grid> grid_);
+ 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;
+ std::string name_;
/// Problem component to work on
- std::string component;
+ std::string component_;
/// Pointer to the grid
- std::shared_ptr<Grid> grid;
+ std::shared_ptr<Grid> grid_;
/// Pointer to the local error estimates
- EstType est;
+ Estimates est_;
/// estimation sum
- double estSum;
+ double estSum_ = 0.0;
/// estmation maximum
- double estMax;
+ double estMax_ = 0.0;
/// If true, elements are marked only if the new value is bigger than
/// the current marking.
- bool maximumMarking;
+ bool maximumMarking_;
/// Lower limit for error estimation, from which an element is marked for
/// refinement
- double markRLimit;
+ double markRLimit_;
/// Upper limit for error estimation, from which an element is marked for
/// coarsening
- double markCLimit;
+ double markCLimit_;
- /// power in estimator norm
- double p;
+ /// power in estimator norm TODO: is this info necessary in marker?
+ double p_ = 2.0;
/// Info level.
- int info;
+ int info_ = 10;
/// Counter for elements marked for refinement
- int elMarkRefine;
+ int elMarkRefine_ = 0;
/// Counter for elements marked for coarsening
- int elMarkCoarsen;
+ int elMarkCoarsen_ = 0;
/// Maximal level of all elements.
- int maxRefineLevel;
+ int maxRefineLevel_ = -1;
/// Minimal level of all elements.
- int minRefineLevel;
+ int minRefineLevel_ = -1;
- bool refineAllowed;
+ bool refineAllowed_ = true;
- bool coarsenAllowed;
+ bool coarsenAllowed_ = false;
};
@@ -180,25 +176,25 @@ namespace AMDiS {
* Global refinement.
*/
template <class Traits>
- class GRMarker : public Marker<Traits>
+ class GRMarker
+ : public Marker<Traits>
{
- 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 EstType = std::vector<double>;
+ 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_, const EstType& est_, std::shared_ptr<Grid> grid_)
- : Marker<Traits>(name_, component_, est_, grid_)
+ 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, const Element& elem)
+ virtual void markElement(AdaptInfo& adaptInfo, Element const& elem)
{
- if (this->refineAllowed)
+ if (this->refineAllowed_)
this->mark(elem, 1);
}
};
@@ -210,26 +206,23 @@ namespace AMDiS {
* \brief
* Maximum strategy.
*/
-
+
template <class Traits>
- class MSMarker : public Marker<Traits>
+ class MSMarker
+ : public Marker<Traits>
{
- 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 EstType = std::vector<double>;
+ using Super = Marker<Traits>;
+ using Grid = typename Super::Grid;
+ using Estimates = typename Super::Estimates;
public:
/// Constructor.
- MSMarker(std::string name_, std::string component_, const EstType& est_, std::shared_ptr<Grid> grid_)
- : Marker<Traits>(name_, component_, est_, grid_),
- MSGamma(0.5),
- MSGammaC(0.1)
+ 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);
+ Parameters::get(this->name_ + "->MSGamma", msGamma_);
+ Parameters::get(this->name_ + "->MSGammaC", msGammaC_);
}
/// Implementation of Marker::initMarking().
@@ -237,10 +230,10 @@ namespace AMDiS {
protected:
/// Marking parameter.
- double MSGamma;
+ double msGamma_ = 0.5;
/// Marking parameter.
- double MSGammaC;
+ double msGammaC_ = 0.1;
};
@@ -250,26 +243,23 @@ namespace AMDiS {
* \brief
* Equidistribution strategy
*/
-
+
template <class Traits>
- class ESMarker : public Marker<Traits>
+ class ESMarker
+ : public Marker<Traits>
{
- 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 EstType = std::vector<double>;
+ using Super = Marker<Traits>;
+ using Grid = typename Super::Grid;
+ using Estimates = typename Super::Estimates;
public:
/// Constructor.
- ESMarker(std::string name_, std::string component_, const EstType& est_, std::shared_ptr<Grid> grid_)
- : Marker<Traits>(name_, component_, est_, grid_),
- ESTheta(0.9),
- ESThetaC(0.2)
+ 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);
+ Parameters::get(this->name_ + "->ESTheta", esTheta_);
+ Parameters::get(this->name_ + "->ESThetaC", esThetaC_);
}
/// Implementation of Marker::initMarking().
@@ -277,10 +267,10 @@ namespace AMDiS {
protected:
/// Marking parameter.
- double ESTheta;
+ double esTheta_ = 0.9;
/// Marking parameter.
- double ESThetaC;
+ double esThetaC_ = 0.2;
};
@@ -290,29 +280,25 @@ namespace AMDiS {
* \brief
* Guaranteed error reduction strategy
*/
-
+
template <class Traits>
- class GERSMarker : public Marker<Traits>
+ class GERSMarker
+ : public Marker<Traits>
{
- 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 EstType = std::vector<double>;
+ 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_, const EstType& est_, std::shared_ptr<Grid> grid_)
- : Marker<Traits>(name_, component_, est_, grid_),
- oldErrSum(0.0),
- GERSThetaStar(0.6),
- GERSNu(0.1),
- GERSThetaC(0.1)
+ 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);
+ Parameters::get(this->name_ + "->GERSThetaStar", gersThetaStar_);
+ Parameters::get(this->name_ + "->GERSNu", gersNu_);
+ Parameters::get(this->name_ + "->GERSThetaC", gersThetaC_);
}
/// Implementation of Marker::markGrid().
@@ -320,28 +306,28 @@ namespace AMDiS {
protected:
/// Refinement marking function.
- void markElementForRefinement(AdaptInfo& adaptInfo, const Element& elem);
+ void markElementForRefinement(AdaptInfo& adaptInfo, Element const& elem);
/// Coarsening marking function.
- void markElementForCoarsening(AdaptInfo& adaptInfo, const Element& elem);
+ void markElementForCoarsening(AdaptInfo& adaptInfo, Element const& elem);
protected:
/// Marking parameter.
- double GERSSum;
+ double gersSum_ = 0.0;
/// Marking parameter.
- double oldErrSum;
+ double oldErrSum_ = 0.0;
/// Marking parameter.
- double GERSThetaStar;
+ double gersThetaStar_ = 0.6;
/// Marking parameter.
- double GERSNu;
+ double gersNu_ = 0.1;
/// Marking parameter.
- double GERSThetaC;
+ double gersThetaC_ = 0.1;
};
}
-#include "Marker.inc.hpp"
\ No newline at end of file
+#include "Marker.inc.hpp"
diff --git a/src/amdis/Marker.inc.hpp b/src/amdis/Marker.inc.hpp
index 346d51e5..7f06e94d 100644
--- a/src/amdis/Marker.inc.hpp
+++ b/src/amdis/Marker.inc.hpp
@@ -1,258 +1,240 @@
-// TODO: Cleanup of copied comments
-
-#include <amdis/common/Math.hpp>
+#include <cmath>
namespace AMDiS {
- using std::pow;
-
- template <class Traits>
- std::shared_ptr<Marker<Traits> > Marker<Traits>::createMarker(std::string name, std::string component, const EstType& est, std::shared_ptr<Grid> 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>
+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 = pow(adaptInfo.getEstSum(component), p);
- estMax = adaptInfo.getEstMax(component);
- refineAllowed = adaptInfo.isRefinementAllowed(component);
- coarsenAllowed = adaptInfo.isCoarseningAllowed(component);
- }
+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>::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))
- if (/*!elem->getElementData()->getElementData(COARSENABLE) ||*/
- lError /*+ elem->getCoarseningEstimation(index)*/ <= markCLimit)
- this->mark(elem, -1);
- }
+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!");
+template <class Traits>
+Flag Marker<Traits>::markGrid(AdaptInfo& adaptInfo)
+{
+ if (!grid_)
+ error_exit("No grid!");
- initMarking(adaptInfo);
+ initMarking(adaptInfo);
- if (!coarsenAllowed &&
- !refineAllowed) {
- return 0;
- }
- for (const auto& elem : Dune::elements(grid->leafGridView())) {
- markElement(adaptInfo, elem);
- }
+ if (!coarsenAllowed_ && !refineAllowed_)
+ return 0;
- finishMarking(adaptInfo);
+ for (const auto& elem : Dune::elements(grid_->leafGridView()))
+ markElement(adaptInfo, elem);
- Flag markFlag;
- if (elMarkRefine)
- markFlag = 1;
- if (elMarkCoarsen)
- markFlag |= 2;
+ finishMarking(adaptInfo);
- return markFlag;
- }
+ Flag markFlag;
+ if (elMarkRefine_)
+ markFlag = 1;
+ if (elMarkCoarsen_)
+ markFlag |= 2;
+ return markFlag;
+}
- template <class Traits>
- void MSMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
- {
- Marker<Traits>::initMarking(adaptInfo);
- double MSGammaP = pow(MSGamma, this->p);
- double MSGammaCP = pow(MSGammaC, this->p);
+template <class Traits>
+void MSMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
+{
+ Super::initMarking(adaptInfo);
- this->markRLimit = MSGammaP * adaptInfo.getEstMax(this->component);
- this->markCLimit = MSGammaCP * adaptInfo.getEstMax(this->component);
+ double msGammaP = std::pow(msGamma_, this->p_);
+ double msGammaCP = std::pow(msGammaC_, this->p_);
- msg("start max_est: ", adaptInfo.getEstMax(this->component), ", mark_limits: ", this->markRLimit, ", " , this->markCLimit);
- }
+ 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)
- {
- Marker<Traits>::initMarking(adaptInfo);
+template <class Traits>
+void ESMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
+{
+ Super::initMarking(adaptInfo);
- double ESThetaP = pow(ESTheta, this->p);
- double ESThetaCP = pow(ESThetaC, this->p);
- double epsP = pow(adaptInfo.getSpaceTolerance(this->component), this->p);
+ 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);
+ int nLeaves = (this->grid_->leafGridView()).size(0);
/*#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
- Parallel::mpi::globalAdd(nLeaves);
+ Parallel::mpi::globalAdd(nLeaves);
#endif*/
- this->markRLimit = ESThetaP * epsP / nLeaves;
- this->markCLimit = ESThetaCP * epsP / nLeaves;
+ this->markRLimit_ = esThetaP * epsP / nLeaves;
+ this->markCLimit_ = esThetaCP * epsP / nLeaves;
- msg("start mark_limits: ", this->markRLimit, ", " , this->markCLimit, "; nt = ", nLeaves);
- }
+ msg("start mark_limits: ", this->markRLimit_, ", " , this->markCLimit_, "; nt = ", nLeaves);
+}
- template <class Traits>
- Flag GERSMarker<Traits>::markGrid(AdaptInfo& adaptInfo)
- {
- Marker<Traits>::initMarking(adaptInfo);
+template <class Traits>
+Flag GERSMarker<Traits>::markGrid(AdaptInfo& adaptInfo)
+{
+ Super::initMarking(adaptInfo);
- if (!this->coarsenAllowed &&
- !this->refineAllowed)
- return 0;
+ if (!this->coarsenAllowed_ && !this->refineAllowed_)
+ return 0;
- GERSSum = 0.0;
+ gersSum_ = 0.0;
- double LTheta = pow(1.0 - GERSThetaStar, this->p);
- double epsP = pow(adaptInfo.getSpaceTolerance(this->component), this->p);
+ 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 (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 = ", pow(LTheta, 1.0 / this->p));
- }
+ 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;
+ oldErrSum_ = this->estSum_;
+ double GERSGamma = 1.0;
- for (const auto& elem : Dune::elements(this->grid->leafGridView()))
- markElementForRefinement(adaptInfo, elem);
+ if (this->refineAllowed_) {
+ if (LTheta > 0) {
+ do {
+ gersSum_ = 0.0;
+ GERSGamma -= gersNu_;
+ this->markRLimit_ = GERSGamma * this->estMax_;
- } while((GERSGamma > 0) && (GERSSum < LTheta * this->estSum));
- }
+ for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
+ markElementForRefinement(adaptInfo, elem);
- msg("GERS refinement with gamma = ", GERSGamma);
+ } while((GERSGamma > 0) && (gersSum_ < LTheta * this->estSum_));
}
- if (this->coarsenAllowed) {
- GERSGamma = 0.3;
- LTheta = GERSThetaC * epsP;
+ msg("GERS refinement with gamma = ", GERSGamma);
+ }
- do {
- GERSSum = 0.0;
- GERSGamma -= GERSNu;
- this->markCLimit = GERSGamma * this->estMax;
+ if (this->coarsenAllowed_) {
+ GERSGamma = 0.3;
+ LTheta = gersThetaC_ * epsP;
- for (const auto& elem : Dune::elements(this->grid->leafGridView()))
- markElementForCoarsening(adaptInfo, elem);
+ do {
+ gersSum_ = 0.0;
+ GERSGamma -= gersNu_;
+ this->markCLimit_ = GERSGamma * this->estMax_;
- msg("coarse loop: gamma = ", GERSGamma, ", sum = ", GERSSum, ", limit = ", LTheta);
- } while(GERSSum > LTheta);
+ for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
+ markElementForCoarsening(adaptInfo, elem);
- msg("GERS coarsening with gamma = ", GERSGamma);
- }
+ msg("coarse loop: gamma = ", GERSGamma, ", sum = ", gersSum_, ", limit = ", LTheta);
+ } while(gersSum_ > LTheta);
+
+ msg("GERS coarsening with gamma = ", GERSGamma);
+ }
- Marker<Traits>::finishMarking(adaptInfo);
+ Super::finishMarking(adaptInfo);
- Flag markFlag;
- if (this->elMarkRefine)
- markFlag = 1;
- if (this->elMarkCoarsen)
- markFlag |= 2;
+ Flag markFlag;
+ if (this->elMarkRefine_)
+ markFlag = 1;
+ if (this->elMarkCoarsen_)
+ markFlag |= 2;
- return(markFlag);
- }
+ return markFlag;
+}
- template <class Traits>
- void GERSMarker<Traits>::markElementForRefinement(AdaptInfo& adaptInfo, const Element& elem)
- {
- double lError = (this->est)[(this->grid->leafIndexSet()).index(elem)];
+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);
- }
+ 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 (elem->getElementData()->getElementData(COARSENABLE))*/
-/* lError += elem->getCoarseningEstimation(index);*/
+template <class Traits>
+void GERSMarker<Traits>::markElementForCoarsening(AdaptInfo& adaptInfo, const Element& elem)
+{
+ double lError = this->est_[this->grid_->leafIndexSet().index(elem)];
- if (lError <= this->markCLimit) {
- GERSSum += lError;
- this->mark(elem, -1);
- } else {
- this->mark(elem, 0);
- }
+ 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 12833ceb..9352089e 100644
--- a/src/amdis/ProblemStat.hpp
+++ b/src/amdis/ProblemStat.hpp
@@ -16,7 +16,6 @@
#include <amdis/Assembler.hpp>
#include <amdis/CreatorInterface.hpp>
#include <amdis/DirichletBC.hpp>
-//#include <amdis/Estimator.hpp>
#include <amdis/Flag.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/LinearAlgebra.hpp>
@@ -145,26 +144,12 @@ namespace AMDiS
bool createMatrixData = true,
bool storeMatrixData = false) override;
-
- /// Implementation of \ref ProblemStatBase::estimate.
- virtual void estimate(AdaptInfo& adaptInfo) override;
-
-
- /// Implementation of \ref ProblemStatBase::markElements.
- virtual Flag markElements(AdaptInfo& adaptInfo) override;
-
-
- /// Implementation of \ref ProblemStatBase::refineMesh.
- virtual Flag refineMesh(AdaptInfo& adaptInfo) 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);
@@ -225,13 +210,7 @@ namespace AMDiS
createFileWriter();
}
-/*
- /// 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_; }
@@ -279,16 +258,12 @@ namespace AMDiS
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");
-
- auto localView = globalBasis->localView();
- forEachNode(localView.tree(), [&,this](auto const& node, auto treePath)
+ 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
+ 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");
@@ -305,30 +280,6 @@ namespace AMDiS
linearSolver_ = solverCreator->create(name_ + "->solver");
}
- void createEstimator()
- {/* COPIED FROM AMDiS
- for (std::size_t i = 0, i < nComponents, i++) {
- std::string estName = "";
- Parameters::get(name + "->estimator->name[" + std::to_string(i) + "]", estName);
- auto estimatorCreator
- = named(CreatorMap<Estimator>::getCreator(estName, name + "->estimator->name[" + std::to_string(i) + "]"));
- if (estimatorCreator) {
- estimatorCreator->setName(estName);
- estimatorCreator->setRow(i);
- estimatorCreator->setSolution(solution->getDOFVector(i));
- estimator[i] = estimatorCreator->create();
- }
-
- if (estimator[i]) {
- for (std::size_t j = 0; j < nComponents; j++)
- estimator[i]->addSystem((*systemMatrix)[i][j],
- solution->getDOFVector(j),
- rhs->getDOFVector(i));
- }
- }*/
- }
-
-
void createMarker();
void createFileWriter();
@@ -366,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:
@@ -390,10 +341,7 @@ namespace AMDiS
std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
/// Pointer to the adaptation markers
- std::list<std::shared_ptr<Marker<Traits> > > marker;
-
- /// Pointer to the estimators for this problem
-// std::vector<Estimator*> estimator;
+ std::list<std::shared_ptr<Marker<Traits> > > marker_;
/// An object of the linearSolver Interface
std::shared_ptr<LinearSolverType> linearSolver_;
@@ -405,8 +353,8 @@ namespace AMDiS
std::shared_ptr<SystemVector> solution_;
/// A vector with the local element error estimates
- /// reverse indexed by [component index][element index]
- std::map<std::string, std::vector<double> > 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
@@ -418,8 +366,6 @@ namespace AMDiS
MatrixOperators<GlobalBasis> matrixOperators_;
VectorOperators<GlobalBasis> rhsOperators_;
Constraints<GlobalBasis> constraints_;
-
-// VectorData<GlobalBasis, TODO: NodeData> estimates_;
};
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp
index c40f6d83..263b5db3 100644
--- a/src/amdis/ProblemStat.inc.hpp
+++ b/src/amdis/ProblemStat.inc.hpp
@@ -79,7 +79,7 @@ void ProblemStat<Traits>::initialize(
if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
solution_ = adoptProblem->solution_;
- estimates = adoptProblem->estimates;
+ estimates_ = adoptProblem->estimates_;
rhs_ = adoptProblem->rhs_;
systemMatrix_ = adoptProblem->systemMatrix_;
}
@@ -103,23 +103,12 @@ void ProblemStat<Traits>::initialize(
warning("no solver created\n");
}
-/*
- // create estimator
- estimator.resize(nComponents, NULL);
-
- if (initFlag.isSet(INIT_ESTIMATOR))
- createEstimator();
-
- if (adoptProblem && adoptFlag.isSet(INIT_ESTIMATOR))
- estimator = adoptProblem->estimator;
-*/
-
// create marker
if (initFlag.isSet(INIT_MARKER))
createMarker();
if (adoptProblem && adoptFlag.isSet(INIT_MARKER))
- marker = adoptProblem->marker;
+ marker_ = adoptProblem->marker_;
// create file writer
@@ -133,27 +122,23 @@ void ProblemStat<Traits>::initialize(
template <class Traits>
void ProblemStat<Traits>::createMarker()
{
- auto localView = globalBasis->localView();
- forEachNode(localView.tree(), [&,this](auto const& node, auto treePath)
+ AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
{
- std::string componentName = name + "->marker[" + to_string(treePath) + "]";
+ std::string componentName = name_ + "->marker[" + to_string(treePath) + "]";
if (!Parameters::get<std::string>(componentName + "->strategy"))
return;
- int i = std::stoi(to_string(treePath)); // TODO: To be removed
- // replace i with treePath once supported
- auto newMarker = Marker<Traits>::createMarker(
- componentName, std::to_string(i), estimates[std::to_string(i)], grid);
+ std::string tp = to_string(treePath);
+ auto newMarker = Marker<Traits>::createMarker(componentName, tp, estimates_[tp], grid_);
if (newMarker)
- {
- marker.push_back(std::move(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: What about two markers each for two grids?
- if (marker.size() > 1 && nGrids == 1)
- marker.back()->setMaximumMarking(true);
+ // TODO: needs to be checked!
+ if (marker_.size() > 1)
+ marker_.back()->setMaximumMarking(true);
});
}
@@ -328,43 +313,13 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
}
-template <class Traits>
-void ProblemStat<Traits>::
-estimate(AdaptInfo& adaptInfo)
-{
- Dune::Timer t;
-/*
- for (std::size_t i = 0; i < nComponents; i++) {
- Estimator *scalEstimator = estimator[i];
-
- if (scalEstimator) {
- scalEstimator->estimate(adaptInfo->getTimestep());
-
- adaptInfo->setEstSum(scalEstimator->getErrorSum(), i);
- adaptInfo->setEstMax(scalEstimator->getErrorMax(), i);
-
- if (adaptInfo->getRosenbrockMode() == false) {
- adaptInfo->setTimeEstSum(scalEstimator->getTimeEst(), i);
- adaptInfo->setTimeEstMax(scalEstimator->getTimeEstMax(), i);
- }
- }
- }
-
-//#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-// MPI::COMM_WORLD.Barrier();
-//#endif
-*/
- msg("estimation of the error needed ", t.elapsed(), "seconds\n");
-}
-
-
template <class Traits>
Flag ProblemStat<Traits>::markElements(AdaptInfo& adaptInfo)
{
Dune::Timer t;
Flag markFlag = 0;
- for (auto currentMarker : marker)
+ for (auto& currentMarker : marker_)
markFlag |= currentMarker->markGrid(adaptInfo);
msg("markElements needed ", t.elapsed(), " seconds");
@@ -373,44 +328,6 @@ Flag ProblemStat<Traits>::markElements(AdaptInfo& adaptInfo)
}
-template <class Traits>
-Flag ProblemStat<Traits>::refineMesh(AdaptInfo& adaptInfo)
-{
- // TODO: actual element data, other components
-
- grid->preAdapt();
-
- std::map<typename Grid::LocalIdSet::IdType, double> vertexData;
-// std::map<typename Grid::LocalIdSet::IdType, double> elementData;
- const auto& gridView = globalBasis->gridView();
- const auto& indexSet = gridView.indexSet();
- const auto& idSet = grid->localIdSet();
-
- // Save data in container during adaptation
- for (const auto& v : vertices(gridView)) {
- vertexData[idSet.id(v)] = (*solution)[indexSet.index(v)];
- }/*
- for (const auto& e : elements(gridView)) {
- elementData[idSet.id(e)] = something[indexSet.index(e)];
- }*/
-
- Flag flag = grid->adapt();
-/*
- // Unpack data
- solution->compress();
- for (const auto& v : vertices(gridView)) {
- (*solution)[indexSet.index(v)] = vertexData[idSet.id(v)];
- }/*
- for (const auto& e : elements(gridView)) {
- something[indexSet.index(e)] = elementData[idSet.id(e)];
- }*/
-
- grid->postAdapt();
-
- return flag;
-}
-
-
template <class Traits>
void ProblemStat<Traits>::
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
@@ -418,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);
diff --git a/src/amdis/utility/TreePath.hpp b/src/amdis/utility/TreePath.hpp
index 9d29ead4..4949a928 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)...
};
}
}
--
GitLab