Merge branch 'feature/marker' into 'develop'

Feature/marker

See merge request spraetor/dune-amdis!20

src/amdis/CMakeLists.txt

@@ -36,6 +36,8 @@ install(FILES
     LinearAlgebra.hpp
     LocalAssembler.hpp
     LocalAssemblerBase.hpp
+    Marker.hpp
+    Marker.inc.hpp
     Mesh.hpp
     Operations.hpp
     Operators.hpp

src/amdis/Marker.hpp

+#pragma once
+
+#include <string>
+
+#include <dune/grid/common/grid.hh>
+
+#include <amdis/AdaptInfo.hpp>
+#include <amdis/Flag.hpp>
+#include <amdis/Initfile.hpp>
+
+namespace AMDiS {
+
+
+  /**
+   * \ingroup Adaption
+   *
+   * \brief
+   * Base class for all scalar markers.
+   */
+  template <class Traits>
+  class Marker
+  {
+  protected:
+    using GlobalBasis = typename Traits::GlobalBasis;
+    using GridView = typename GlobalBasis::GridView;
+    using Grid     = typename GridView::Grid;
+    using IndexSet = typename GridView::IndexSet;
+    using Element  = typename GridView::template Codim<0>::Entity;
+    using Estimates = std::vector<double>;
+
+  public:
+    /// Constructor.
+    Marker() {}
+
+    /// Constructor.
+    Marker(std::string name, std::string component, Estimates const& est,
+           std::shared_ptr<Grid> const& grid)
+      : name_(name)
+      , component_(component)
+      , grid_(grid)
+      , est_(est)
+      , maximumMarking_(false)
+    {
+      Parameters::get(name_ + "->p", p_);
+      Parameters::get(name_ + "->info", info_);
+      Parameters::get(name_ + "->max refinement level", maxRefineLevel_);
+      Parameters::get(name_ + "->min refinement level", minRefineLevel_);
+    }
+
+    /// destructor
+    virtual ~Marker() {}
+
+    /// Marks element with newMark. If \ref maximumMarking_ is set, the element
+    /// is marked only if the new mark is bigger than the old one. The return
+    /// value specifies whether the element has been marked, or not.
+    void mark(Element const& elem, int newMark)
+    {
+      int oldMark = grid_->getMark(elem);
+
+      if (!maximumMarking_ || (newMark > oldMark)) {
+        bool marked = grid_->mark(newMark, elem);
+        if (marked) {
+          if (oldMark > 0) {
+            elMarkRefine_--;
+          } else if (oldMark < 0) {
+            elMarkCoarsen_--;
+          }
+
+          if (newMark > 0) {
+            elMarkRefine_++;
+          } else if (newMark < 0) {
+              elMarkCoarsen_++;
+          }
+        } else {
+          msg("Marking failed");
+        }
+      }
+    }
+
+    /// Can be used by sub classes. Called before traversal.
+    virtual void initMarking(AdaptInfo& adaptInfo);
+
+    /// Can be used by sub classes. Called after traversal.
+    virtual void finishMarking(AdaptInfo& adaptInfo);
+
+    /// Marks one element.
+    virtual void markElement(AdaptInfo& adaptInfo, Element const& elem);
+
+    /// Marking of the mesh.
+    virtual Flag markGrid(AdaptInfo& adaptInfo);
+
+    /// Sets \ref maximumMarking.
+    void setMaximumMarking(bool maxMark)
+    {
+      maximumMarking_ = maxMark;
+    }
+
+    int getElMarkRefine()
+    {
+      return elMarkRefine_;
+    }
+
+    int getElMarkCoarsen()
+    {
+      return elMarkCoarsen_;
+    }
+
+    /// Returns \ref name of the Marker
+    std::string getName() const
+    {
+      return name_;
+    }
+
+    /// Creates a scalar marker depending on the strategy set in parameters.
+    static std::shared_ptr<Marker<Traits> > createMarker(std::string name, std::string component, Estimates const& est, std::shared_ptr<Grid> const& grid);
+
+  protected:
+    /// Name of the scalar marker.
+    std::string name_;
+
+    /// Problem component to work on
+    std::string component_;
+
+    /// Pointer to the grid
+    std::shared_ptr<Grid> grid_;
+
+    /// Pointer to the local error estimates
+    Estimates est_;
+
+    /// estimation sum
+    double estSum_ = 0.0;
+
+    /// estmation maximum
+    double estMax_ = 0.0;
+
+    /// If true, elements are marked only if the new value is bigger than
+    /// the current marking.
+    bool maximumMarking_;
+
+    /// Lower limit for error estimation, from which an element is marked for
+    /// refinement
+    double markRLimit_;
+
+    /// Upper limit for error estimation, from which an element is marked for
+    /// coarsening
+    double markCLimit_;
+
+    /// power in estimator norm TODO: is this info necessary in marker?
+    double p_ = 2.0;
+
+    /// Info level.
+    int info_ = 10;
+
+    /// Counter for elements marked for refinement
+    int elMarkRefine_ = 0;
+
+    /// Counter for elements marked for coarsening
+    int elMarkCoarsen_ = 0;
+
+    /// Maximal level of all elements.
+    int maxRefineLevel_ = -1;
+
+    /// Minimal level of all elements.
+    int minRefineLevel_ = -1;
+
+    bool refineAllowed_ = true;
+
+    bool coarsenAllowed_ = false;
+  };
+
+
+  /**
+   * \ingroup Adaption
+   *
+   * \brief
+   * Global refinement.
+   */
+  template <class Traits>
+  class GRMarker
+      : public Marker<Traits>
+  {
+    using Super = Marker<Traits>;
+    using Grid      = typename Super::Grid;
+    using Element   = typename Super::Element;
+    using Estimates = typename Super::Estimates;
+
+  public:
+    /// Constructor.
+    GRMarker(std::string name, std::string component, Estimates const& est,
+             std::shared_ptr<Grid> const& grid)
+      : Marker<Traits>(std::move(name), std::move(component), est, grid)
+    {}
+
+    /// Implementation of Marker::markElement().
+    virtual void markElement(AdaptInfo& adaptInfo, Element const& elem)
+    {
+      if (this->refineAllowed_)
+        this->mark(elem, 1);
+    }
+  };
+
+
+  /**
+   * \ingroup Adaption
+   *
+   * \brief
+   * Maximum strategy.
+   */
+
+  template <class Traits>
+  class MSMarker
+      : public Marker<Traits>
+  {
+    using Super = Marker<Traits>;
+    using Grid      = typename Super::Grid;
+    using Estimates = typename Super::Estimates;
+
+  public:
+    /// Constructor.
+    MSMarker(std::string name, std::string component, Estimates const& est,
+             std::shared_ptr<Grid> const& grid)
+      : Super{std::move(name), std::move(component), est, grid}
+    {
+      Parameters::get(this->name_ + "->MSGamma", msGamma_);
+      Parameters::get(this->name_ + "->MSGammaC", msGammaC_);
+    }
+
+    /// Implementation of Marker::initMarking().
+    void initMarking(AdaptInfo& adaptInfo);
+
+  protected:
+    /// Marking parameter.
+    double msGamma_ = 0.5;
+
+    /// Marking parameter.
+    double msGammaC_ = 0.1;
+  };
+
+
+  /**
+   * \ingroup Adaption
+   *
+   * \brief
+   * Equidistribution strategy
+   */
+
+  template <class Traits>
+  class ESMarker
+      : public Marker<Traits>
+  {
+    using Super = Marker<Traits>;
+    using Grid      = typename Super::Grid;
+    using Estimates = typename Super::Estimates;
+
+  public:
+    /// Constructor.
+    ESMarker(std::string name, std::string component, Estimates const& est,
+             std::shared_ptr<Grid> const& grid)
+      : Super{std::move(name), std::move(component), est, grid}
+    {
+      Parameters::get(this->name_ + "->ESTheta", esTheta_);
+      Parameters::get(this->name_ + "->ESThetaC", esThetaC_);
+    }
+
+    /// Implementation of Marker::initMarking().
+    virtual void initMarking(AdaptInfo& adaptInfo);
+
+  protected:
+    /// Marking parameter.
+    double esTheta_ = 0.9;
+
+    /// Marking parameter.
+    double esThetaC_ = 0.2;
+  };
+
+
+  /**
+   * \ingroup Adaption
+   *
+   * \brief
+   * Guaranteed error reduction strategy
+   */
+
+  template <class Traits>
+  class GERSMarker
+      : public Marker<Traits>
+  {
+    using Super = Marker<Traits>;
+    using Grid      = typename Super::Grid;
+    using Element   = typename Super::Element;
+    using Estimates = typename Super::Estimates;
+
+  public:
+    /// Constructor.
+    GERSMarker(std::string name, std::string component, Estimates const& est,
+               std::shared_ptr<Grid> const& grid)
+      : Super{std::move(name), std::move(component), est, grid}
+    {
+      Parameters::get(this->name_ + "->GERSThetaStar", gersThetaStar_);
+      Parameters::get(this->name_ + "->GERSNu", gersNu_);
+      Parameters::get(this->name_ + "->GERSThetaC", gersThetaC_);
+    }
+
+    /// Implementation of Marker::markGrid().
+    virtual Flag markGrid(AdaptInfo& adaptInfo);
+
+  protected:
+    /// Refinement marking function.
+    void markElementForRefinement(AdaptInfo& adaptInfo, Element const& elem);
+
+    /// Coarsening marking function.
+    void markElementForCoarsening(AdaptInfo& adaptInfo, Element const& elem);
+
+  protected:
+    /// Marking parameter.
+    double gersSum_ = 0.0;
+
+    /// Marking parameter.
+    double oldErrSum_ = 0.0;
+
+    /// Marking parameter.
+    double gersThetaStar_ = 0.6;
+
+    /// Marking parameter.
+    double gersNu_ = 0.1;
+
+    /// Marking parameter.
+    double gersThetaC_ = 0.1;
+  };
+
+}
+
+#include "Marker.inc.hpp"

src/amdis/Marker.inc.hpp

+#include <cmath>
+
+namespace AMDiS {
+
+template <class Traits>
+std::shared_ptr<Marker<Traits> > Marker<Traits>::createMarker(std::string name, std::string component, Estimates const& est, std::shared_ptr<Grid> const& grid)
+{
+  int strategy = 0;
+  Parameters::get(name + "->strategy", strategy);
+
+  std::shared_ptr<Marker<Traits> > marker;
+
+  switch (strategy) {
+  case 0:  // no refinement/coarsening
+    break;
+  case 1:
+    msg("Creating global refinement (GR) marker");
+    marker = std::make_shared<GRMarker<Traits> >(name, component, est, grid);
+    break;
+  case 2:
+    msg("Creating maximum strategy (MS) marker");
+    marker = std::make_shared<MSMarker<Traits> >(name, component, est, grid);
+    break;
+  case 3:
+    msg("Creating equidistribution strategy (ES) marker");
+    marker = std::make_shared<ESMarker<Traits> >(name, component, est, grid);
+    break;
+  case 4:
+    msg("Creating guaranteed error reduction strategy (GERS) marker");
+    marker = std::make_shared<GERSMarker<Traits> >(name, component, est, grid);
+    break;
+  default:
+    error_exit("invalid strategy");
+  }
+
+  return marker;
+}
+
+
+template <class Traits>
+void Marker<Traits>::initMarking(AdaptInfo& adaptInfo)
+{
+  elMarkRefine_ = 0;
+  elMarkCoarsen_ = 0;
+  estSum_ = std::pow(adaptInfo.getEstSum(component_), p_);
+  estMax_ = adaptInfo.getEstMax(component_);
+  refineAllowed_ = adaptInfo.isRefinementAllowed(component_);
+  coarsenAllowed_ = adaptInfo.isCoarseningAllowed(component_);
+}
+
+
+template <class Traits>
+void Marker<Traits>::finishMarking(AdaptInfo& adaptInfo)
+{
+  msg(elMarkRefine_, " elements marked for refinement");
+  msg(elMarkCoarsen_, " elements marked for coarsening");
+}
+
+
+template <class Traits>
+void Marker<Traits>::markElement(AdaptInfo& adaptInfo, const Element& elem)
+{
+  const auto& index = grid_->leafIndexSet().index(elem);
+  double lError = est_[index];
+
+  if (lError > markRLimit_ && refineAllowed_
+      && (maxRefineLevel_ == -1 || elem.level() < maxRefineLevel_)) {
+    this->mark(elem, 1);
+  } else  if (lError <= markCLimit_ && coarsenAllowed_
+              && (minRefineLevel_ == -1 || elem.level() > minRefineLevel_)) {
+    this->mark(elem, -1);
+  }
+}
+
+
+template <class Traits>
+Flag Marker<Traits>::markGrid(AdaptInfo& adaptInfo)
+{
+  if (!grid_)
+    error_exit("No grid!");
+
+  initMarking(adaptInfo);
+
+  if (!coarsenAllowed_ && !refineAllowed_)
+    return 0;
+
+  for (const auto& elem : Dune::elements(grid_->leafGridView()))
+    markElement(adaptInfo, elem);
+
+  finishMarking(adaptInfo);
+
+  Flag markFlag;
+  if (elMarkRefine_)
+    markFlag = 1;
+  if (elMarkCoarsen_)
+    markFlag |= 2;
+
+  return markFlag;
+}
+
+
+template <class Traits>
+void MSMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
+{
+  Super::initMarking(adaptInfo);
+
+  double msGammaP = std::pow(msGamma_, this->p_);
+  double msGammaCP = std::pow(msGammaC_, this->p_);
+
+  this->markRLimit_ = msGammaP * adaptInfo.getEstMax(this->component_);
+  this->markCLimit_ = msGammaCP * adaptInfo.getEstMax(this->component_);
+
+  msg("start max_est: ", adaptInfo.getEstMax(this->component_), ", mark_limits: ", this->markRLimit_, ", " , this->markCLimit_);
+}
+
+
+template <class Traits>
+void ESMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
+{
+  Super::initMarking(adaptInfo);
+
+  double esThetaP = std::pow(esTheta_, this->p_);
+  double esThetaCP = std::pow(esThetaC_, this->p_);
+  double epsP = std::pow(adaptInfo.getSpaceTolerance(this->component_), this->p_);
+
+  int nLeaves = (this->grid_->leafGridView()).size(0);
+/*#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+  Parallel::mpi::globalAdd(nLeaves);
+#endif*/
+
+  this->markRLimit_ = esThetaP * epsP / nLeaves;
+  this->markCLimit_ = esThetaCP * epsP / nLeaves;
+
+  msg("start mark_limits: ", this->markRLimit_, ", " , this->markCLimit_, "; nt = ", nLeaves);
+}
+
+
+template <class Traits>
+Flag GERSMarker<Traits>::markGrid(AdaptInfo& adaptInfo)
+{
+  Super::initMarking(adaptInfo);
+
+  if (!this->coarsenAllowed_ && !this->refineAllowed_)
+    return 0;
+
+  gersSum_ = 0.0;
+
+  double LTheta = std::pow(1.0 - gersThetaStar_, this->p_);
+  double epsP = std::pow(adaptInfo.getSpaceTolerance(this->component_), this->p_);
+
+  if (this->estSum_ < oldErrSum_) {
+    double improv = this->estSum_ / oldErrSum_;
+    double wanted = 0.8 * epsP / this->estSum_;
+    double redfac = std::min((1.0 - wanted) / (1.0 - improv), 1.0);
+    redfac = std::max(redfac, 0.0);
+
+    if (redfac < 1.0) {
+      LTheta *= redfac;
+      msg("GERS: use extrapolated theta_star = ", std::pow(LTheta, 1.0 / this->p_));
+    }
+  }
+
+  oldErrSum_ = this->estSum_;
+  double GERSGamma = 1.0;
+
+  if (this->refineAllowed_) {
+    if (LTheta > 0) {
+      do {
+        gersSum_ = 0.0;
+        GERSGamma -= gersNu_;
+        this->markRLimit_ = GERSGamma * this->estMax_;
+
+        for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
+          markElementForRefinement(adaptInfo, elem);
+
+      } while((GERSGamma > 0) && (gersSum_ < LTheta * this->estSum_));
+    }
+
+    msg("GERS refinement with gamma = ", GERSGamma);
+  }
+
+  if (this->coarsenAllowed_) {
+    GERSGamma = 0.3;
+    LTheta = gersThetaC_ * epsP;
+
+    do {
+      gersSum_ = 0.0;
+      GERSGamma -= gersNu_;
+      this->markCLimit_ = GERSGamma * this->estMax_;
+
+      for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
+        markElementForCoarsening(adaptInfo, elem);
+
+      msg("coarse loop: gamma = ", GERSGamma, ", sum = ", gersSum_, ", limit = ", LTheta);
+    } while(gersSum_ > LTheta);
+
+    msg("GERS coarsening with gamma = ", GERSGamma);
+  }
+
+  Super::finishMarking(adaptInfo);
+
+  Flag markFlag;
+  if (this->elMarkRefine_)
+    markFlag = 1;
+  if (this->elMarkCoarsen_)
+    markFlag |= 2;
+
+  return markFlag;
+}
+
+
+template <class Traits>
+void GERSMarker<Traits>::markElementForRefinement(AdaptInfo& adaptInfo, const Element& elem)
+{
+  double lError = this->est_[this->grid_->leafIndexSet().index(elem)];
+
+  if (lError > this->markRLimit_) {
+    gersSum_ += lError;
+    this->mark(elem, 1);
+  }
+}
+
+
+template <class Traits>
+void GERSMarker<Traits>::markElementForCoarsening(AdaptInfo& adaptInfo, const Element& elem)
+{
+  double lError = this->est_[this->grid_->leafIndexSet().index(elem)];
+
+  if (this->grid_->getMark(elem) <= 0) {
+    if (lError <= this->markCLimit_) {
+      gersSum_ += lError;
+      this->mark(elem, -1);
+    } else {
+      this->mark(elem, 0);
+    }
+  }
+}
+
+
+} // end namespace AMDiS

src/amdis/ProblemStat.hpp

@@ -19,6 +19,7 @@
 #include <amdis/Flag.hpp>
 #include <amdis/Initfile.hpp>
 #include <amdis/LinearAlgebra.hpp>
+#include <amdis/Marker.hpp>
 #include <amdis/Mesh.hpp>
 #include <amdis/ProblemStatBase.hpp>
 #include <amdis/ProblemStatTraits.hpp>
@@ -33,6 +34,7 @@
 
 #include <amdis/io/FileWriterInterface.hpp>
 
+#include <amdis/utility/TreeData.hpp>
 #include <amdis/utility/TreePath.hpp>
 
 namespace AMDiS
@@ -52,7 +54,7 @@ namespace AMDiS
     using Element  = typename GridView::template Codim<0>::Entity;
     using WorldVector = typename Element::Geometry::GlobalCoordinate;
 
-    /// Dimension of the mesh
+    /// Dimension of the grid
     static constexpr int dim = Grid::dimension;
 
     /// Dimension of the world
@@ -66,15 +68,15 @@ namespace AMDiS
   public:
     /**
      * \brief Constructor. Takes the name of the problem that is used to
-     * access values correpsonding to this püroblem in the parameter file.
+     * access values corresponding to this problem in the parameter file.
      **/
     explicit ProblemStat(std::string name)
       : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
       , name_(std::move(name))
     {}
 
-    /// Constructor taking additionally a reference to a mesh that is used
-    /// instead of the default created mesh, \ref ProblemStat
+    /// Constructor taking additionally a reference to a grid that is used
+    /// instead of the default created grid, \ref ProblemStat
     ProblemStat(std::string name, Grid& grid)
       : ProblemStat(std::move(name))
     {
@@ -142,14 +144,12 @@ namespace A