GridFunctionMarker

examples/init/marker.dat.2d

+dimension of world: 2
+
+test->mesh: testMesh
+testMesh->global refinements: 0

src/amdis/Marker.hpp

 #pragma once
 
+#include <limits>
 #include <string>
+#include <utility>
 
+#include <dune/common/std/optional.hh>
 #include <dune/grid/common/grid.hh>
 
+#include <amdis/common/ConceptsBase.hpp>
+
+#include <amdis/gridfunctions/GridFunctionConcepts.hpp>
+
 #include <amdis/AdaptInfo.hpp>
 #include <amdis/Flag.hpp>
 #include <amdis/Initfile.hpp>
 
-namespace AMDiS {
-
-
+namespace AMDiS
+{
   /**
    * \ingroup Adaption
    *
    * \brief
-   * Base class for all scalar markers.
+   * Base class for all markers.
    */
-  template <class Traits>
+  template <class Grid>
   class Marker
   {
   protected:
-    using GlobalBasis = typename Traits::GlobalBasis;
-    using GridView = typename GlobalBasis::GridView;
-    using Grid     = typename GridView::Grid;
-    using IndexSet = typename GridView::IndexSet;
+    using GridView = typename Grid::LeafGridView;
     using Element  = typename GridView::template Codim<0>::Entity;
-    using Estimates = std::vector<double>;
 
   public:
     /// Constructor.
-    Marker() {}
+    Marker() = default;
 
     /// Constructor.
-    Marker(std::string const& name, std::string const& component, Estimates const& est,
-           std::shared_ptr<Grid> const& grid)
+    Marker(std::string const& name, 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() {}
+    /// Destructor.
+    virtual ~Marker() = default;
 
     /// 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");
-        }
-      }
-    }
+    void mark(Element const& elem, int newMark);
 
-    /// Can be used by sub classes. Called before traversal.
+    /// Called before traversal.
     virtual void initMarking(AdaptInfo& adaptInfo);
 
-    /// Can be used by sub classes. Called after traversal.
+    /// Called after traversal.
     virtual void finishMarking(AdaptInfo& adaptInfo);
 
     /// Marks one element.
-    virtual void markElement(AdaptInfo& adaptInfo, Element const& elem);
+    virtual void markElement(AdaptInfo& adaptInfo, Element const& elem) = 0;
 
-    /// Marking of the mesh.
+    /// Marking of the whole grid.
     virtual Flag markGrid(AdaptInfo& adaptInfo);
 
-    /// Sets \ref maximumMarking.
-    void setMaximumMarking(bool maxMark)
-    {
-      maximumMarking_ = maxMark;
-    }
-
-    int getElMarkRefine()
+    /// Returns \ref elMarkRefine_ of the Marker
+    int getElMarkRefine() const
     {
       return elMarkRefine_;
     }
 
-    int getElMarkCoarsen()
+    /// Returns \ref elMarkCoarsen_ of the Marker
+    int getElMarkCoarsen() const
     {
       return elMarkCoarsen_;
     }
 
-    /// Returns \ref name of the Marker
-    std::string const& getName() const
+    /// 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::unique_ptr<Marker<Traits> > createMarker(std::string const& name, std::string const& component, Estimates const& est, std::shared_ptr<Grid> const& grid);
+    /// Sets \ref maximumMarking_.
+    void setMaximumMarking(bool maxMark)
+    {
+      maximumMarking_ = maxMark;
+    }
+
+    /// Sets \ref refineAllowed_.
+    void allowRefinement(bool allow)
+    {
+      refineAllowed_ = allow;
+    }
+
+    /// Sets \ref coarsenAllowed_.
+    void allowCoarsening(bool allow)
+    {
+      coarsenAllowed_ = allow;
+    }
 
   protected:
-    /// Name of the scalar marker.
+    /// Name of the 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;
+    bool maximumMarking_ = false;
 
     /// Info level.
-    int info_ = 10;
+    int info_ = 0;
 
     /// Counter for elements marked for refinement
     int elMarkRefine_ = 0;
@@ -163,9 +126,76 @@ namespace AMDiS {
     /// Minimal level of all elements.
     int minRefineLevel_ = -1;
 
+    /// Allow elements to be marked for refinement
     bool refineAllowed_ = true;
 
-    bool coarsenAllowed_ = false;
+    /// Allow elements to be marked for coarsening
+    bool coarsenAllowed_ = true;
+  };
+
+
+  /**
+   * \ingroup Adaption
+   *
+   * \brief
+   * Base class for all estimator-based markers.
+   */
+  template <class Grid>
+  class EstimatorMarker
+      : public Marker<Grid>
+  {
+  protected:
+    using Super = Marker<Grid>;
+    using Element   = typename Super::Element;
+    using Estimates = std::vector<double>;
+
+  public:
+    /// Constructor.
+    EstimatorMarker() = default;
+
+    /// Constructor.
+    EstimatorMarker(std::string name, std::string component, Estimates const& est,
+                    std::shared_ptr<Grid> const& grid)
+      : Super{std::move(name), grid}
+      , component_(std::move(component))
+      , est_(est)
+    {
+      Parameters::get(this->name_ + "->p", p_);
+    }
+
+    /// Can be used by sub classes. Called before traversal.
+    virtual void initMarking(AdaptInfo& adaptInfo) override;
+
+    /// Marks one element.
+    virtual void markElement(AdaptInfo& adaptInfo, Element const& elem) override;
+
+    /// Creates a scalar marker depending on the strategy set in parameters.
+    static std::unique_ptr<EstimatorMarker<Grid> > createMarker(std::string const& name,
+      std::string const& component, Estimates const& est, std::shared_ptr<Grid> const& grid);
+
+  protected:
+    /// Problem component to work on
+    std::string component_;
+
+    /// Pointer to the local error estimates
+    Estimates est_;
+
+    /// Estimation sum
+    double estSum_ = 0.0;
+
+    /// Estmation maximum
+    double estMax_ = 0.0;
+
+    /// 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_ = 2.0;
   };
 
 
@@ -175,24 +205,22 @@ namespace AMDiS {
    * \brief
    * Global refinement.
    */
-  template <class Traits>
+  template <class Grid>
   class GRMarker
-      : public Marker<Traits>
+      : public EstimatorMarker<Grid>
   {
-    using Super = Marker<Traits>;
-    using Grid      = typename Super::Grid;
+    using Super = EstimatorMarker<Grid>;
     using Element   = typename Super::Element;
     using Estimates = typename Super::Estimates;
 
   public:
     /// Constructor.
-    GRMarker(std::string const& name, std::string component, Estimates const& est,
+    GRMarker(std::string const& name, std::string const& component, Estimates const& est,
              std::shared_ptr<Grid> const& grid)
-      : Marker<Traits>(name, std::move(component), est, grid)
+      : Super{name, component, est, grid}
     {}
 
-    /// Implementation of Marker::markElement().
-    virtual void markElement(AdaptInfo& adaptInfo, Element const& elem)
+    virtual void markElement(AdaptInfo& adaptInfo, Element const& elem) override
     {
       if (this->refineAllowed_)
         this->mark(elem, 1);
@@ -207,12 +235,11 @@ namespace AMDiS {
    * Maximum strategy.
    */
 
-  template <class Traits>
+  template <class Grid>
   class MSMarker
-      : public Marker<Traits>
+      : public EstimatorMarker<Grid>
   {
-    using Super = Marker<Traits>;
-    using Grid      = typename Super::Grid;
+    using Super = EstimatorMarker<Grid>;
     using Estimates = typename Super::Estimates;
 
   public:
@@ -225,8 +252,7 @@ namespace AMDiS {
       Parameters::get(name + "->MSGammaC", msGammaC_);
     }
 
-    /// Implementation of Marker::initMarking().
-    void initMarking(AdaptInfo& adaptInfo);
+    virtual void initMarking(AdaptInfo& adaptInfo) override;
 
   protected:
     /// Marking parameter.
@@ -244,12 +270,11 @@ namespace AMDiS {
    * Equidistribution strategy
    */
 
-  template <class Traits>
+  template <class Grid>
   class ESMarker
-      : public Marker<Traits>
+      : public EstimatorMarker<Grid>
   {
-    using Super = Marker<Traits>;
-    using Grid      = typename Super::Grid;
+    using Super = EstimatorMarker<Grid>;
     using Estimates = typename Super::Estimates;
 
   public:
@@ -262,8 +287,7 @@ namespace AMDiS {
       Parameters::get(name + "->ESThetaC", esThetaC_);
     }
 
-    /// Implementation of Marker::initMarking().
-    virtual void initMarking(AdaptInfo& adaptInfo);
+    virtual void initMarking(AdaptInfo& adaptInfo) override;
 
   protected:
     /// Marking parameter.
@@ -281,12 +305,11 @@ namespace AMDiS {
    * Guaranteed error reduction strategy
    */
 
-  template <class Traits>
+  template <class Grid>
   class GERSMarker
-      : public Marker<Traits>
+      : public EstimatorMarker<Grid>
   {
-    using Super = Marker<Traits>;
-    using Grid      = typename Super::Grid;
+    using Super = EstimatorMarker<Grid>;
     using Element   = typename Super::Element;
     using Estimates = typename Super::Estimates;
 
@@ -301,8 +324,7 @@ namespace AMDiS {
       Parameters::get(name + "->GERSThetaC", gersThetaC_);
     }
 
-    /// Implementation of Marker::markGrid().
-    virtual Flag markGrid(AdaptInfo& adaptInfo);
+    virtual Flag markGrid(AdaptInfo& adaptInfo) override;
 
   protected:
     /// Refinement marking function.
@@ -328,6 +350,86 @@ namespace AMDiS {
     double gersThetaC_ = 0.1;
   };
 
-}
+
+  /**
+   * \ingroup Adaption
+   *
+   * \brief
+   * Marker based on an indicator given as grid-function.
+   *
+   * On each element the grid-function is evaluated in the barycenter. The returned
+   * values must be an integer (or convertible to an integer) indicating the
+   * desired refinement level of this element. The element is marked for refinement
+   * if the current level is < the desired level or coarsened, if >.
+   *
+   * \tparam Grid     An Implementation of the \ref Dune::Grid interface
+   * \tparam GridFct  A grid-function with `Range` convertible to `int`
+   */
+  template <class Grid, class GridFct>
+  class GridFunctionMarker
+      : public Marker<Grid>
+  {
+    using Super = Marker<Grid>;
+    using Element = typename Super::Element;
+
+    template <class GF>
+    using IsGridFunction = decltype(localFunction(std::declval<GF>()));
+
+    static_assert(Dune::Std::is_detected<IsGridFunction,GridFct>::value, "");
+
+  public:
+    /// Constructor.
+    template <class GF>
+    GridFunctionMarker(std::string const& name, std::shared_ptr<Grid> const& grid, GF&& gf,
+                       Dune::Std::optional<int> maxRef = { /*max*/ },
+                       Dune::Std::optional<int> minRef = { /*0*/ })
+      : Super{name, grid}
+      , gridFct_{makeGridFunction(std::forward<GF>(gf), grid->leafGridView())}
+    {
+      
+      this->maxRefineLevel_ = maxRef ? maxRef.value() : 
+        (this->maxRefineLevel_ == -1 ? std::numeric_limits<int>::max() : this->maxRefineLevel_);
+      this->minRefineLevel_ = minRef ? minRef.value() : 
+        (this->minRefineLevel_ == -1 ? 0 : this->minRefineLevel_);
+    }
+
+    /// \brief Implementation of \ref Marker::markElement. Does nothing since marking is
+    /// done in \ref markGrid().
+    virtual void markElement(AdaptInfo& adaptInfo, Element const& elem) final {}
+
+    /// Mark element for refinement, if grid-function \ref gridFct_ evaluates to
+    /// a value larger than the current level and is marked for coarsening of
+    /// the result is less than the current value.
+    virtual Flag markGrid(AdaptInfo& adaptInfo) override;
+
+  private:
+    /// Indicator function for adaptation
+    GridFct gridFct_;
+  };
+
+#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
+  // Deduction guide for GridFunctionMarker class
+  template <class Grid, class PreGridFct>
+  GridFunctionMarker(std::string const& name, std::shared_ptr<Grid> const& grid,
+                     PreGridFct&& preGridFct,
+                     Dune::Std::optional<int> maxRef = {},
+                     Dune::Std::optional<int> minRef = {})
+    -> GridFunctionMarker<Grid,
+          std::decay_t<decltype(makeGridFunction(std::forward<PreGridFct>(preGridFct), grid->leafGridView()))>>;
+#endif
+
+  // Generator function for GridFunctionMarker class
+  template <class Grid, class PreGridFct>
+  auto makeGridFunctionMarker(std::string const& name, std::shared_ptr<Grid> const& grid,
+                              PreGridFct&& preGridFct,
+                              Dune::Std::optional<int> maxRef = { /*max*/ },
+                              Dune::Std::optional<int> minRef = { /*0*/ })
+  {
+    auto gridFct = makeGridFunction(std::forward<PreGridFct>(preGridFct), grid->leafGridView());
+    using GridFct = decltype(gridFct);
+    return GridFunctionMarker<Grid,GridFct>{name, grid, gridFct, maxRef, minRef};
+  }
+
+} // end namespace AMDiS
 
 #include "Marker.inc.hpp"

src/amdis/Marker.inc.hpp

@@ -2,105 +2,134 @@
 
 namespace AMDiS {
 
-template <class Traits>
-std::unique_ptr<Marker<Traits> > Marker<Traits>::
-createMarker(std::string const& name, std::string const& component,
-             Estimates const& est, std::shared_ptr<Grid> const& grid)
+template <class Grid>
+void Marker<Grid>::mark(Element const& elem, int newMark)
 {
-  int strategy = 0;
-  Parameters::get(name + "->strategy", strategy);
-
-  switch (strategy) {
-  case 0:  // no refinement/coarsening
-    break;
-  case 1:
-    msg("Creating global refinement (GR) marker");
-    return std::make_unique<GRMarker<Traits> >(name, component, est, grid);
-    break;
-  case 2:
-    msg("Creating maximum strategy (MS) marker");
-    return std::make_unique<MSMarker<Traits> >(name, component, est, grid);
-    break;
-  case 3:
-    msg("Creating equidistribution strategy (ES) marker");
-    return std::make_unique<ESMarker<Traits> >(name, component, est, grid);
-    break;
-  case 4:
-    msg("Creating guaranteed error reduction strategy (GERS) marker");
-    return std::make_unique<GERSMarker<Traits> >(name, component, est, grid);
-    break;
-  default:
-    error_exit("invalid strategy");
+  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");
+    }
   }
-
-  return {};
-}
-
-
-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)
+template <class Grid>
+void Marker<Grid>::initMarking(AdaptInfo& adaptInfo)
 {
-  msg("{} elements marked for refinement", elMarkRefine_);
-  msg("{} elements marked for coarsening", elMarkCoarsen_);
+  this->elMarkRefine_ = 0;
+  this->elMarkCoarsen_ = 0;
 }
 
 
-template <class Traits>
-void Marker<Traits>::markElement(AdaptInfo& adaptInfo, const Element& elem)
+template <class Grid>
+void Marker<Grid>::finishMarking(AdaptInfo& adaptInfo)
 {
-  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);
+  if (info_ > 0) {
+    msg("{} elements marked for refinement", elMarkRefine_);
+    msg("{} elements marked for coarsening", elMarkCoarsen_);
   }
 }
 
 
-template <class Traits>
-Flag Marker<Traits>::markGrid(AdaptInfo& adaptInfo)
+template <class Grid>
+Flag Marker<Grid>::markGrid(AdaptInfo& adaptInfo)
 {
-  if (!grid_)
-    error_exit("No grid!");
+  test_exit(bool(this->grid_), "No grid!");
 
   initMarking(adaptInfo);
 
-  if (!coarsenAllowed_ && !refineAllowed_)
+  if (!this->coarsenAllowed_ && !this->refineAllowed_)
     return 0;
 
-  for (const auto& elem : Dune::elements(grid_->leafGridView()))
+  for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
     markElement(adaptInfo, elem);
 
   finishMarking(adaptInfo);
 
   Flag markFlag;
-  if (elMarkRefine_)
+  if (this->elMarkRefine_)
     markFlag = 1;
-  if (elMarkCoarsen_)
+  if (this->elMarkCoarsen_)
     markFlag |= 2;
 
   return markFlag;
 }
 
+template <class Grid>
+void EstimatorMarker<Grid>::initMarking(AdaptInfo& adaptInfo)
+{
+  Super::initMarking(adaptInfo);
+  estSum_ = std::pow(adaptInfo.getEstSum(component_), p_);
+  estMax_ = adaptInfo.getEstMax(component_);
+  this->refineAllowed_ = adaptInfo.isRefinementAllowed(component_);
+  this->coarsenAllowed_ = adaptInfo.isCoarseningAllowed(component_);