Added class Marker based on AMDiS Marker class

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

src/amdis/Marker.cpp

+// 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);
+      }
+    }
+  }
+
+
+}

src/amdis/Marker.hpp

+#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;
+  };
+
+}