diff --git a/src/amdis/Marker.hpp b/src/amdis/Marker.hpp index 628c6b4278d050be6b949ed46435de916e974d23..75bcad474b7e6e016476b911744a9f6b47e9dfb8 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 346d51e5335c2cb10a5a69cb3cdb1c4d01ff9e34..7f06e94d711ef5514d9ba1b33242da5394ca03c7 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 12833cebc30d674634e01eef9c3b2380ccc34080..9352089e83f137088cde0b6136cc5d3922baf007 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 c40f6d83b0fbd32bb953bd6d9185fe10c93b9a83..263b5db3d35308019eb2b2eb6146853381c79815 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 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)... }; } }