From fb74a3a39b15b9ab8c72cae985ab79dde6b5120a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Felix=20M=C3=BCller?= <felix.mueller2@mailbox.tu-dresden.de>
Date: Fri, 23 Mar 2018 15:11:38 +0100
Subject: [PATCH] Added class Marker based on AMDiS Marker class
---
src/amdis/CMakeLists.txt | 2 +
src/amdis/Marker.cpp | 278 ++++++++++++++++++++++++++++++++
src/amdis/Marker.hpp | 333 +++++++++++++++++++++++++++++++++++++++
3 files changed, 613 insertions(+)
create mode 100644 src/amdis/Marker.cpp
create mode 100644 src/amdis/Marker.hpp
diff --git a/src/amdis/CMakeLists.txt b/src/amdis/CMakeLists.txt
index 45e875a5..d2067008 100644
--- a/src/amdis/CMakeLists.txt
+++ b/src/amdis/CMakeLists.txt
@@ -8,6 +8,7 @@ dune_library_add_sources(amdis SOURCES
AMDiS.cpp
Initfile.cpp
InitfileParser.cpp
+ Marker.cpp
ProblemInstatBase.cpp
# ProblemInstat.cpp
ProblemStat.cpp
@@ -36,6 +37,7 @@ install(FILES
LinearAlgebra.hpp
LocalAssembler.hpp
LocalAssemblerBase.hpp
+ Marker.hpp
Mesh.hpp
Operations.hpp
Operators.hpp
diff --git a/src/amdis/Marker.cpp b/src/amdis/Marker.cpp
new file mode 100644
index 00000000..877c30b5
--- /dev/null
+++ b/src/amdis/Marker.cpp
@@ -0,0 +1,278 @@
+// TODO: Cleanup of copied comments
+
+#include "Marker.hpp"
+
+#include <amdis/common/Math.hpp>
+
+namespace AMDiS {
+
+ using std::pow;
+
+ template <class Grid>
+ Marker<Grid>* Marker<Grid>::createMarker(std::string name, int row, const EstType& est, Grid* grid)
+ {
+ int strategy = 0;
+ Parameters::get(name + "->strategy", strategy);
+
+ Marker *marker = NULL;
+
+ switch (strategy) {
+ case 0: // no refinement/coarsening
+ break;
+ case 1:
+ msg("Creating global refinement (GR) marker\n");
+ marker = new GRMarker<Grid>(name, row, est, grid);
+ break;
+ case 2:
+ msg("Creating maximum strategy (MS) marker\n");
+ marker = new MSMarker<Grid>(name, row, est, grid);
+ break;
+ case 3:
+ msg("Creating equidistribution strategy (ES) marker\n");
+ marker = new ESMarker<Grid>(name, row, est, grid);
+ break;
+ case 4:
+ msg("Creating guaranteed error reduction strategy (GERS) marker\n");
+ marker = new GERSMarker<Grid>(name, row, est, grid);
+ break;
+ default:
+ error_exit("invalid strategy\n");
+ }
+
+ return marker;
+ }
+
+
+
+ template <class Grid>
+ void Marker<Grid>::initMarking(AdaptInfo& adaptInfo)
+ {
+ int row_ = row == -1 ? 0 : row;
+
+ elMarkRefine = 0;
+ elMarkCoarsen = 0;
+ estSum = pow(adaptInfo.getEstSum(row_), p);
+ estMax = adaptInfo.getEstMax(row_);
+ }
+
+
+
+ template <class Grid>
+ void Marker<Grid>::finishMarking(AdaptInfo& adaptInfo)
+ {
+ msg(elMarkRefine, " elements marked for refinement\n");
+ msg(elMarkCoarsen, " elements marked for coarsening\n");
+ }
+
+
+
+ template <class Grid>
+ void Marker<Grid>::markElement(AdaptInfo& adaptInfo, const Element& elem)
+ {
+ msg("DEBUG: Marker::markElement\n");
+ int row_ = row == -1 ? 0 : row;
+ const auto& index = grid->leafIndexSet().index(elem);
+ double lError = est[index][row_];
+
+ if (adaptInfo.isRefinementAllowed(row_) && lError > markRLimit) {
+ if (maxRefineLevel == -1 || elem.level() < maxRefineLevel) {
+ this->mark(elem, 1);
+ }
+ } else {
+ if (adaptInfo.isCoarseningAllowed(row_) && lError <= markCLimit) {
+ if (minRefineLevel == -1 || elem.level() > minRefineLevel) {
+ if (/*!elem->getElementData()->getElementData(COARSENABLE) ||*/
+ lError /*+ elem->getCoarseningEstimation(row)*/ <= markCLimit) {
+ this->mark(elem, -1);
+ }
+ }
+ }
+ }
+ }
+
+
+
+ template <class Grid>
+ Flag Marker<Grid>::markGrid(AdaptInfo& adaptInfo)
+ {
+ msg("DEBUG: Marker::markGrid\n");
+
+ test_exit(grid, "No grid!\n");
+
+ initMarking(adaptInfo);
+
+ int row_ = row == -1 ? 0 : row;
+
+ if (!adaptInfo.isCoarseningAllowed(row_) &&
+ !adaptInfo.isRefinementAllowed(row_)) {
+ 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 Grid>
+ void MSMarker<Grid>::initMarking(AdaptInfo& adaptInfo)
+ {
+ Marker<Grid>::initMarking(adaptInfo);
+
+ int row_ = this->row == -1 ? 0 : this->row;
+
+ double MSGammaP = pow(MSGamma, this->p);
+ double MSGammaCP = pow(MSGammaC, this->p);
+
+ this->markRLimit = MSGammaP * adaptInfo.getEstMax(row_);
+ this->markCLimit = MSGammaCP * adaptInfo.getEstMax(row_);
+
+ msg("start max_est: ", adaptInfo.getEstMax(row_), ", mark_limits: ", this->markRLimit, ", " , this->markCLimit, "\n");
+ }
+
+
+ template <class Grid>
+ void ESMarker<Grid>::initMarking(AdaptInfo& adaptInfo)
+ {
+ Marker<Grid>::initMarking(adaptInfo);
+
+ int row_ = this->row == -1 ? 0 : this->row;
+
+ double ESThetaP = pow(ESTheta, this->p);
+ double ESThetaCP = pow(ESThetaC, this->p);
+ double epsP = pow(adaptInfo.getSpaceTolerance(row_), 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 , "\n");
+ }
+
+
+ template <class Grid>
+ Flag GERSMarker<Grid>::markGrid(AdaptInfo& adaptInfo)
+ {
+ msg("DEBUG: GERSMarker::markGrid\n");
+
+ Marker<Grid>::initMarking(adaptInfo);
+
+ int row_ = this->row == -1 ? 0 : this->row;
+
+ if (!adaptInfo.isCoarseningAllowed(row_) &&
+ !adaptInfo.isRefinementAllowed(row_))
+ return 0;
+
+ GERSSum = 0.0;
+
+ double LTheta = pow(1.0 - GERSThetaStar, this->p);
+ double epsP = pow(adaptInfo.getSpaceTolerance(row_), 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 = ", pow(LTheta, 1.0 / this->p), "\n");
+ }
+ }
+
+ oldErrSum = this->estSum;
+ double GERSGamma = 1.0;
+
+ if (adaptInfo.isRefinementAllowed(row_)) {
+ 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, "\n");
+ }
+
+ if (adaptInfo.isCoarseningAllowed(row_)) {
+ 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, "\n");
+ } while(GERSSum > LTheta);
+
+ msg("GERS coarsening with gamma = ", GERSGamma, "\n");
+ }
+
+ Marker<Grid>::finishMarking(adaptInfo);
+
+ Flag markFlag;
+ if (this->elMarkRefine)
+ markFlag = 1;
+ if (this->elMarkCoarsen)
+ markFlag |= 2;
+
+ return(markFlag);
+ }
+
+
+ template <class Grid>
+ void GERSMarker<Grid>::markElementForRefinement(AdaptInfo& adaptInfo, const Element& elem)
+ {
+ int row_ = this->row == -1 ? 0 : this->row;
+ double lError = (this->est)[(this->grid->leafIndexSet()).index(elem)][row_];
+
+ if (lError > this->markRLimit) {
+ GERSSum += lError;
+ this->mark(elem, 1);
+ }
+ }
+
+
+ template <class Grid>
+ void GERSMarker<Grid>::markElementForCoarsening(AdaptInfo& adaptInfo, const Element& elem)
+ {
+ int row_ = this->row == -1 ? 0 : this->row;
+ double lError = (this->est)[(this->grid->leafIndexSet()).index(elem)][row_];
+
+ if ((this->grid)->getMark(elem) <= 0) {
+/* if (elem->getElementData()->getElementData(COARSENABLE))*/
+/* lError += elem->getCoarseningEstimation(row);*/
+
+ if (lError <= this->markCLimit) {
+ GERSSum += lError;
+ this->mark(elem, -1);
+ } else {
+ this->mark(elem, 0);
+ }
+ }
+ }
+
+
+}
diff --git a/src/amdis/Marker.hpp b/src/amdis/Marker.hpp
new file mode 100644
index 00000000..9817e216
--- /dev/null
+++ b/src/amdis/Marker.hpp
@@ -0,0 +1,333 @@
+#pragma once
+// TODO: Cleanup of copied comments
+#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 Grid>
+ class Marker
+ {
+
+ using GridView = typename Grid::LeafGridView;
+ using IndexSet = typename GridView::IndexSet;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using EstType = std::vector<std::vector<double> >;
+
+ public:
+ /// Constructor.
+ Marker() {}
+
+ /// Constructor.
+ Marker(std::string name_, int row_, const EstType& est_, Grid* grid_)
+ : name(name_),
+ row(row_),
+ grid(grid_),
+ est(est_),
+ maximumMarking(false),
+ p(2),
+ info(10),
+ maxRefineLevel(-1),
+ minRefineLevel(-1)
+ {
+ 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.
+ inline void mark(const Element& elem, int newMark)
+ {
+ int oldMark = grid->getMark(elem);
+
+ if (!maximumMarking || (newMark > oldMark)) {
+ bool (marked) = grid->mark(newMark, elem);
+ msg("DEBUG: Was marked? ", marked, "; oldMark = ", oldMark, ", newMark = ", newMark, "\n");
+ 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, const Element& elem);
+
+ /// Marking of the mesh.
+ virtual Flag markGrid(AdaptInfo& adaptInfo);
+
+ /// Sets \ref maximumMarking.
+ inline void setMaximumMarking(bool maxMark)
+ {
+ maximumMarking = maxMark;
+ }
+
+ inline int getElMarkRefine()
+ {
+ return elMarkRefine;
+ }
+
+ inline int getElMarkCoarsen()
+ {
+ return elMarkCoarsen;
+ }
+
+ /// Returns \ref name of the Marker
+ inline std::string getName() const
+ {
+ return name;
+ }
+
+ /// Creates a scalar marker depending on the strategy set in parameters.
+ static Marker<Grid>* createMarker(std::string name, int row_, const EstType& est_, Grid* grid_);
+
+ protected:
+ /// Name of the scalar marker.
+ std::string name;
+
+ /// Equal to -1 for scalar problems. Component number if marker is
+ /// part of a vector valued marker.
+ int row;
+
+ /// Pointer to the grid
+ Grid* grid;
+
+ /// Pointer to the local error estimates
+ EstType est;
+
+ /// estimation sum
+ double estSum;
+
+ /// estmation maximum
+ double estMax;
+
+ /// 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
+ double p;
+
+ /// Info level.
+ int info;
+
+ /// Counter for elements marked for refinement
+ int elMarkRefine;
+
+ /// Counter for elements marked for coarsening
+ int elMarkCoarsen;
+
+ /// Maximal level of all elements.
+ int maxRefineLevel;
+
+ /// Minimal level of all elements.
+ int minRefineLevel;
+ };
+
+
+ /**
+ * \ingroup Adaption
+ *
+ * \brief
+ * Global refinement.
+ */
+ template <class Grid>
+ class GRMarker : public Marker<Grid>
+ {
+ using GridView = typename Grid::LeafGridView;
+ using IndexSet = typename GridView::IndexSet;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using EstType = std::vector<std::vector<double> >;
+
+ public:
+ /// Constructor.
+ GRMarker(std::string name_, int row_, const EstType& est_, Grid* grid_)
+ : Marker<Grid>(name_, row_, est_, grid_)
+ {}
+
+ /// Implementation of Marker::markElement().
+ virtual void markElement(AdaptInfo& adaptInfo, const Element& elem)
+ {
+ msg("DEBUG: GRMarker::markElement\n");
+ if (adaptInfo.isRefinementAllowed(this->row == -1 ? 0 : this->row))
+ this->mark(elem, 1);
+ }
+ };
+
+
+ /**
+ * \ingroup Adaption
+ *
+ * \brief
+ * Maximum strategy.
+ */
+
+ template <class Grid>
+ class MSMarker : public Marker<Grid>
+ {
+ using GridView = typename Grid::LeafGridView;
+ using IndexSet = typename GridView::IndexSet;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using EstType = std::vector<std::vector<double> >;
+
+ public:
+ /// Constructor.
+ MSMarker(std::string name_, int row_, const EstType& est_, Grid* grid_)
+ : Marker<Grid>(name_, row_, est_, grid_),
+ MSGamma(0.5),
+ MSGammaC(0.1)
+ {
+ 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;
+
+ /// Marking parameter.
+ double MSGammaC;
+ };
+
+
+ /**
+ * \ingroup Adaption
+ *
+ * \brief
+ * Equidistribution strategy
+ */
+
+ template <class Grid>
+ class ESMarker : public Marker<Grid>
+ {
+ using GridView = typename Grid::LeafGridView;
+ using IndexSet = typename GridView::IndexSet;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using EstType = std::vector<std::vector<double> >;
+
+ public:
+ /// Constructor.
+ ESMarker(std::string name_, int row_, const EstType& est_, Grid* grid_)
+ : Marker<Grid>(name_, row_, est_, grid_),
+ ESTheta(0.9),
+ ESThetaC(0.2)
+ {
+ 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;
+
+ /// Marking parameter.
+ double ESThetaC;
+ };
+
+
+ /**
+ * \ingroup Adaption
+ *
+ * \brief
+ * Guaranteed error reduction strategy
+ */
+
+ template <class Grid>
+ class GERSMarker : public Marker<Grid>
+ {
+ using GridView = typename Grid::LeafGridView;
+ using IndexSet = typename GridView::IndexSet;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using EstType = std::vector<std::vector<double> >;
+
+ public:
+ /// Constructor.
+ GERSMarker(std::string name_, int row_, const EstType& est_, Grid* grid_)
+ : Marker<Grid>(name_, row_, est_, grid_),
+ oldErrSum(0.0),
+ GERSThetaStar(0.6),
+ GERSNu(0.1),
+ GERSThetaC(0.1)
+ {
+ 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, const Element& elem);
+
+ /// Coarsening marking function.
+ void markElementForCoarsening(AdaptInfo& adaptInfo, const Element& elem);
+
+ protected:
+ /// Marking parameter.
+ double GERSSum;
+
+ /// Marking parameter.
+ double oldErrSum;
+
+ /// Marking parameter.
+ double GERSThetaStar;
+
+ /// Marking parameter.
+ double GERSNu;
+
+ /// Marking parameter.
+ double GERSThetaC;
+ };
+
+}
--
GitLab