From ee65b18e3d92c9b96ad537d98c41da11fcebdc50 Mon Sep 17 00:00:00 2001
From: "Praetorius, Simon" <simon.praetorius@tu-dresden.de>
Date: Fri, 19 Oct 2018 18:46:05 +0200
Subject: [PATCH] GridFunctionMarker
---
examples/init/marker.dat.2d | 4 +
src/amdis/Marker.hpp | 326 ++++++++++++++++---------
src/amdis/Marker.inc.hpp | 229 +++++++++++------
src/amdis/ProblemStat.hpp | 22 +-
src/amdis/ProblemStat.inc.hpp | 19 +-
src/amdis/ProblemStatBase.hpp | 5 +-
src/amdis/StandardProblemIteration.cpp | 4 +-
test/CMakeLists.txt | 4 +
test/MarkerTest.cpp | 59 +++++
9 files changed, 468 insertions(+), 204 deletions(-)
create mode 100644 examples/init/marker.dat.2d
create mode 100644 test/MarkerTest.cpp
diff --git a/examples/init/marker.dat.2d b/examples/init/marker.dat.2d
new file mode 100644
index 00000000..7ca1df7e
--- /dev/null
+++ b/examples/init/marker.dat.2d
@@ -0,0 +1,4 @@
+dimension of world: 2
+
+test->mesh: testMesh
+testMesh->global refinements: 0
diff --git a/src/amdis/Marker.hpp b/src/amdis/Marker.hpp
index ba2fa1fd..2941248b 100644
--- a/src/amdis/Marker.hpp
+++ b/src/amdis/Marker.hpp
@@ -1,155 +1,118 @@
#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"
diff --git a/src/amdis/Marker.inc.hpp b/src/amdis/Marker.inc.hpp
index e9f08a9a..9bd74e68 100644
--- a/src/amdis/Marker.inc.hpp
+++ b/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_);
+}
+
-template <class Traits>
-void MSMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
+template <class Grid>
+void EstimatorMarker<Grid>::markElement(AdaptInfo& adaptInfo, const Element& elem)
+{
+ const auto& index = this->grid_->leafIndexSet().index(elem);
+ double lError = est_[index];
+
+ if (lError > markRLimit_ && this->refineAllowed_
+ && (this->maxRefineLevel_ == -1 || elem.level() < this->maxRefineLevel_)) {
+ this->mark(elem, 1);
+ } else if (lError <= markCLimit_ && this->coarsenAllowed_
+ && (this->minRefineLevel_ == -1 || elem.level() > this->minRefineLevel_)) {
+ this->mark(elem, -1);
+ }
+}
+
+
+template <class Grid>
+std::unique_ptr<EstimatorMarker<Grid>> EstimatorMarker<Grid>::
+createMarker(std::string const& name, std::string const& component,
+ Estimates const& est, std::shared_ptr<Grid> const& grid)
+{
+ int strategy = 0;
+ Parameters::get(name + "->strategy", strategy);
+
+ switch (strategy) {
+ case 0: // no refinement/coarsening
+ break;
+ case 1:
+ return std::make_unique<GRMarker<Grid> >(name, component, est, grid);
+ break;
+ case 2:
+ return std::make_unique<MSMarker<Grid> >(name, component, est, grid);
+ break;
+ case 3:
+ return std::make_unique<ESMarker<Grid> >(name, component, est, grid);
+ break;
+ case 4:
+ return std::make_unique<GERSMarker<Grid> >(name, component, est, grid);
+ break;
+ default:
+ error_exit("invalid strategy");
+ }
+
+ return {};
+}
+
+
+template <class Grid>
+void MSMarker<Grid>::initMarking(AdaptInfo& adaptInfo)
{
Super::initMarking(adaptInfo);
@@ -115,8 +144,8 @@ void MSMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
}
-template <class Traits>
-void ESMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
+template <class Grid>
+void ESMarker<Grid>::initMarking(AdaptInfo& adaptInfo)
{
Super::initMarking(adaptInfo);
@@ -125,9 +154,9 @@ void ESMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
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*/
+#if AMDIS_HAS_PARALLEL
+ Dune::CollectiveCommunication<>{}.sum(nLeaves, 1);
+#endif
this->markRLimit_ = esThetaP * epsP / nLeaves;
this->markCLimit_ = esThetaCP * epsP / nLeaves;
@@ -136,8 +165,8 @@ void ESMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
}
-template <class Traits>
-Flag GERSMarker<Traits>::markGrid(AdaptInfo& adaptInfo)
+template <class Grid>
+Flag GERSMarker<Grid>::markGrid(AdaptInfo& adaptInfo)
{
Super::initMarking(adaptInfo);
@@ -210,8 +239,8 @@ Flag GERSMarker<Traits>::markGrid(AdaptInfo& adaptInfo)
}
-template <class Traits>
-void GERSMarker<Traits>::markElementForRefinement(AdaptInfo& adaptInfo, const Element& elem)
+template <class Grid>
+void GERSMarker<Grid>::markElementForRefinement(AdaptInfo& adaptInfo, const Element& elem)
{
double lError = this->est_[this->grid_->leafIndexSet().index(elem)];
@@ -222,8 +251,8 @@ void GERSMarker<Traits>::markElementForRefinement(AdaptInfo& adaptInfo, const El
}
-template <class Traits>
-void GERSMarker<Traits>::markElementForCoarsening(AdaptInfo& adaptInfo, const Element& elem)
+template <class Grid>
+void GERSMarker<Grid>::markElementForCoarsening(AdaptInfo& adaptInfo, const Element& elem)
{
double lError = this->est_[this->grid_->leafIndexSet().index(elem)];
@@ -238,4 +267,48 @@ void GERSMarker<Traits>::markElementForCoarsening(AdaptInfo& adaptInfo, const El
}
+template <class Grid, class PreGridFct>
+Flag GridFunctionMarker<Grid, PreGridFct>::markGrid(AdaptInfo& adaptInfo)
+{
+ test_exit(bool(this->grid_), "No grid!");
+
+ Super::initMarking(adaptInfo);
+
+ if (!this->coarsenAllowed_ && !this->refineAllowed_)
+ return 0;
+
+ auto localFct = localFunction(gridFct_);
+
+ for (auto const& e : Dune::elements(this->grid_->leafGridView())) {
+ localFct.bind(e);
+ int currentLevel = e.level();
+ auto geo = e.geometry();
+ auto const& ref = Dune::referenceElement(geo);
+ int targetLevel = 0;
+ int codim = ref.dimension;
+
+ for (int i = 0; i < ref.size(codim); i++)
+ targetLevel = std::max(targetLevel, (int)std::round(localFct(ref.position(i, codim))));
+
+ int m = ((((targetLevel > currentLevel) && (currentLevel < this->maxRefineLevel_))
+ || (currentLevel < this->minRefineLevel_))
+ && this->refineAllowed_)
+ - ((((targetLevel < currentLevel) && (currentLevel > this->minRefineLevel_))
+ || (currentLevel > this->maxRefineLevel_))
+ && this->coarsenAllowed_);
+ this->mark(e, m);
+ localFct.unbind();
+ }
+
+ Super::finishMarking(adaptInfo);
+
+ Flag markFlag;
+ if (this->elMarkRefine_)
+ markFlag = 1;
+ if (this->elMarkCoarsen_)
+ markFlag |= 2;
+
+ return markFlag;
+}
+
} // end namespace AMDiS
diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp
index 8509e68b..74eb5c10 100644
--- a/src/amdis/ProblemStat.hpp
+++ b/src/amdis/ProblemStat.hpp
@@ -16,6 +16,7 @@
#include <amdis/CreatorInterface.hpp>
#include <amdis/CreatorMap.hpp>
#include <amdis/DirichletBC.hpp>
+//#include <amdis/Estimator.hpp>
#include <amdis/Flag.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/LinearAlgebra.hpp>
@@ -216,16 +217,28 @@ namespace AMDiS
Grid const& grid() const { return *grid_; }
/// Set the grid. Stores pointer and initializes feSpaces
- /// matrices and vectors, as well as the file-writer.
+ /// matrices and vectors, as well as markers and file-writers.
void setGrid(std::shared_ptr<Grid> const& grid)
{
grid_ = grid;
-
createGlobalBasis();
createMatricesAndVectors();
+ createMarker();
createFileWriter();
}
+ void addMarker(std::shared_ptr<Marker<Grid>> const& marker)
+ {
+ marker_.push_back(marker);
+ if (marker_.size() > 1)
+ marker_.back()->setMaximumMarking(true);
+ }
+
+ void addMarker(Marker<Grid>& marker)
+ {
+ addMarker(Dune::stackobject_to_shared_ptr(marker));
+ }
+
/// Return the gridView of the leaf-level
GridView const& gridView() { return globalBasis_->gridView(); }
@@ -363,7 +376,10 @@ namespace AMDiS
std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
/// Pointer to the adaptation markers
- std::list<std::shared_ptr<Marker<Traits> > > marker_;
+ std::list<std::shared_ptr<Marker<Grid> > > marker_;
+
+ /// Pointer to the estimators for this problem
+// std::vector<Estimator*> estimator;
/// An object of the linearSolver Interface
std::shared_ptr<LinearSolverType> linearSolver_;
diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp
index 56b904eb..e11503ba 100644
--- a/src/amdis/ProblemStat.inc.hpp
+++ b/src/amdis/ProblemStat.inc.hpp
@@ -120,6 +120,7 @@ void ProblemStat<Traits>::initialize(
template <class Traits>
void ProblemStat<Traits>::createMarker()
{
+ marker_.clear();
auto localView = globalBasis_->localView();
AMDiS::forEachNode_(localView.tree(), [&,this](auto const& node, auto treePath)
{
@@ -129,7 +130,7 @@ void ProblemStat<Traits>::createMarker()
return;
std::string tp = to_string(treePath);
- auto newMarker = Marker<Traits>::createMarker(componentName, tp, estimates_[tp], grid_);
+ auto newMarker = EstimatorMarker<Grid>::createMarker(componentName, tp, estimates_[tp], grid_);
if (newMarker)
marker_.push_back(std::move(newMarker));
@@ -145,6 +146,7 @@ void ProblemStat<Traits>::createMarker()
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
{
+ filewriter_.clear();
auto localView = globalBasis_->localView();
forEachNode_(localView.tree(), [&,this](auto const& node, auto treePath)
{
@@ -223,7 +225,8 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
template <class Traits>
-Flag ProblemStat<Traits>::markElements(AdaptInfo& adaptInfo)
+Flag ProblemStat<Traits>::
+markElements(AdaptInfo& adaptInfo)
{
Dune::Timer t;
@@ -244,12 +247,18 @@ adaptGrid(AdaptInfo& adaptInfo)
Dune::Timer t;
grid_->preAdapt();
- grid_->adapt();
- globalBasis_->update(gridView());
+ bool changed = grid_->adapt();
+
+ // update data
+ if (changed) {
+ globalBasis_->update(gridView());
+ solution_->resize(*globalBasis_);
+ }
+
grid_->postAdapt();
msg("adaptGrid needed {} seconds", t.elapsed());
- return 0;
+ return changed ? MESH_ADAPTED : Flag(0);
}
diff --git a/src/amdis/ProblemStatBase.hpp b/src/amdis/ProblemStatBase.hpp
index 20245e2b..9ba0b0b6 100644
--- a/src/amdis/ProblemStatBase.hpp
+++ b/src/amdis/ProblemStatBase.hpp
@@ -40,8 +40,7 @@ namespace AMDiS
INIT_ADAPT | INIT_FILEWRITER | INIT_INITIAL_PROBLEM |
INIT_UH_OLD | INIT_UPDATER | INIT_NONLIN_SOLVER ;
- const Flag MESH_REFINED = 1;
- const Flag MESH_COARSENED = 2;
+ const Flag MESH_ADAPTED = 1;
// forward declaration
class AdaptInfo;
@@ -70,7 +69,7 @@ namespace AMDiS
* matrices or vectors only.
*/
virtual void buildAfterAdapt(AdaptInfo& adaptInfo, Flag flag,
- bool assembleMatrix, bool assembleVector) = 0;
+ bool assembleMatrix, bool assembleVector) = 0;
/// Refinement/coarsening of the grid.
virtual Flag adaptGrid(AdaptInfo& adaptInfo) = 0;
diff --git a/src/amdis/StandardProblemIteration.cpp b/src/amdis/StandardProblemIteration.cpp
index b65dfc8f..1c7047d6 100644
--- a/src/amdis/StandardProblemIteration.cpp
+++ b/src/amdis/StandardProblemIteration.cpp
@@ -63,10 +63,8 @@ namespace AMDiS
if (toDo.isSet(MARK))
markFlag = problem_.markElements(adaptInfo);
- else
- markFlag = 3;
- if (toDo.isSet(ADAPT) && (markFlag.isSet(MESH_REFINED) || markFlag.isSet(MESH_COARSENED)))
+ if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_ADAPTED))
flag |= problem_.adaptGrid(adaptInfo);
if (toDo.isSet(BUILD))
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index afb4a3cc..a18847f2 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -24,6 +24,10 @@ dune_add_test(SOURCES FiniteElementTypeTest.cpp
dune_add_test(SOURCES FilesystemTest.cpp
LINK_LIBRARIES amdis)
+dune_add_test(SOURCES MarkerTest.cpp
+ LINK_LIBRARIES amdis
+ CMD_ARGS "${CMAKE_SOURCE_DIR}/examples/init/marker.dat.2d")
+
dune_add_test(SOURCES MultiTypeVectorTest.cpp
LINK_LIBRARIES amdis)
diff --git a/test/MarkerTest.cpp b/test/MarkerTest.cpp
new file mode 100644
index 00000000..f34e0f75
--- /dev/null
+++ b/test/MarkerTest.cpp
@@ -0,0 +1,59 @@
+#include <amdis/AMDiS.hpp>
+#include <amdis/ProblemStat.hpp>
+
+#include "Tests.hpp"
+
+using namespace AMDiS;
+
+static const int d = 2; // problem dimension
+
+using Grid = Dune::UGGrid<d>;
+using Basis = LagrangeBasis<Grid::LeafGridView, 1>;
+using Problem = ProblemStat<Basis>;
+using DomainType = typename Dune::FieldVector<double,d>;
+
+int main(int argc, char** argv)
+{
+ AMDiS::init(argc, argv);
+
+ DomainType lowerLeft; lowerLeft = 0.0; // lower left grid corner
+ DomainType upperRight; upperRight = 1.0; // upper right grid corner
+ std::array<unsigned int, d> s; s.fill(1); // number of elements on each axis
+
+ std::shared_ptr<Grid> grid =
+ Dune::StructuredGridFactory<Grid>::createCubeGrid(lowerLeft, upperRight, s);
+
+ Problem prob("test", *grid);
+ prob.initialize(INIT_ALL);
+
+ AdaptInfo adaptInfo("adapt");
+
+ // refine the grid close to the points (0, 0) and (1, 1)
+ auto markerFunc = [](auto const& x) -> double
+ {
+ static Dune::FieldVector<double, d> x0{0.0, 0.0};
+ static Dune::FieldVector<double, d> x1{1.0, 1.0};
+ return { 1.0 / (0.1 + std::min(distance(x0, x), distance(x1, x))) };
+ };
+
+#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
+ GridFunctionMarker marker("mymarker", grid, markerFunc);
+#else
+ auto marker = makeGridFunctionMarker("mymarker", grid, markerFunc);
+#endif
+
+ prob.addMarker(marker);
+
+ for (int i = 0; i < 10; ++i) {
+ prob.markElements(adaptInfo);
+ if (!prob.adaptGrid(adaptInfo))
+ break;
+ }
+
+ prob.getSolution().interpolate(markerFunc);
+ prob.writeFiles(adaptInfo);
+
+ AMDiS::finalize();
+
+ return report_errors();
+}
--
GitLab