AMDiS  2.10
The Adaptive Multi-Dimensional Simulation Toolbox
Marker.inc.hpp
1 #include <cmath>
2 
3 namespace AMDiS {
4 
5 template <class Grid>
6 void Marker<Grid>::mark(Element const& elem, int newMark)
7 {
8  int oldMark = grid_->getMark(elem);
9 
10  if (!maximumMarking_ || (newMark > oldMark)) {
11  bool marked = grid_->mark(newMark, elem);
12  if (marked) {
13  if (oldMark > 0) {
14  elMarkRefine_--;
15  } else if (oldMark < 0) {
16  elMarkCoarsen_--;
17  }
18 
19  if (newMark > 0) {
20  elMarkRefine_++;
21  } else if (newMark < 0) {
22  elMarkCoarsen_++;
23  }
24  } else {
25  msg("Marking failed");
26  }
27  }
28 }
29 
30 
31 template <class Grid>
33 {
34  this->elMarkRefine_ = 0;
35  this->elMarkCoarsen_ = 0;
36 }
37 
38 
39 template <class Grid>
41 {
42  if (info_ > 0) {
43  msg("{} elements marked for refinement", elMarkRefine_);
44  msg("{} elements marked for coarsening", elMarkCoarsen_);
45  }
46 }
47 
48 
49 template <class Grid>
51 {
52  test_exit(bool(this->grid_), "No grid!");
53 
54  initMarking(adaptInfo);
55 
56  if (!this->coarsenAllowed_ && !this->refineAllowed_)
57  return 0;
58 
59  for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
60  markElement(adaptInfo, elem);
61 
62  finishMarking(adaptInfo);
63 
64  Flag markFlag;
65  if (this->elMarkRefine_)
66  markFlag = 1;
67  if (this->elMarkCoarsen_)
68  markFlag |= 2;
69 
70  return markFlag;
71 }
72 
73 template <class Grid>
75 {
76  Super::initMarking(adaptInfo);
77  estSum_ = std::pow(adaptInfo.estSum(component_), p_);
78  estMax_ = adaptInfo.estMax(component_);
79  this->refineAllowed_ = adaptInfo.isRefinementAllowed(component_);
80  this->coarsenAllowed_ = adaptInfo.isCoarseningAllowed(component_);
81 }
82 
83 
84 template <class Grid>
85 void EstimatorMarker<Grid>::markElement(AdaptInfo& /*adaptInfo*/, const Element& elem)
86 {
87  const auto& index = this->grid_->leafIndexSet().index(elem);
88  double lError = est_[index];
89 
90  if (lError > markRLimit_ && this->refineAllowed_
91  && elem.level() < this->maxRefineLevel_) {
92  this->mark(elem, 1);
93  } else if (lError <= markCLimit_ && this->coarsenAllowed_
94  && elem.level() > this->minRefineLevel_) {
95  this->mark(elem, -1);
96  }
97 }
98 
99 
100 template <class Grid>
101 std::unique_ptr<EstimatorMarker<Grid>> EstimatorMarker<Grid>::
102 createMarker(std::string const& name, std::string const& component,
103  Estimates const& est, std::shared_ptr<Grid> const& grid)
104 {
105  int strategy = 0;
106  Parameters::get(name + "->strategy", strategy);
107 
108  switch (strategy) {
109  case 0: // no refinement/coarsening
110  break;
111  case 1:
112  return std::make_unique<GRMarker<Grid> >(name, component, est, grid);
113  break;
114  case 2:
115  return std::make_unique<MSMarker<Grid> >(name, component, est, grid);
116  break;
117  case 3:
118  return std::make_unique<ESMarker<Grid> >(name, component, est, grid);
119  break;
120  case 4:
121  return std::make_unique<GERSMarker<Grid> >(name, component, est, grid);
122  break;
123  default:
124  error_exit("invalid strategy");
125  }
126 
127  return {};
128 }
129 
130 
131 template <class Grid>
133 {
134  Super::initMarking(adaptInfo);
135 
136  double msGammaP = std::pow(msGamma_, this->p_);
137  double msGammaCP = std::pow(msGammaC_, this->p_);
138 
139  this->markRLimit_ = msGammaP * adaptInfo.estMax(this->component_);
140  this->markCLimit_ = msGammaCP * adaptInfo.estMax(this->component_);
141 
142  msg("start max_est: {}, mark_limits: {}, {}",
143  adaptInfo.estMax(this->component_), this->markRLimit_ , this->markCLimit_);
144 }
145 
146 
147 template <class Grid>
149 {
150  Super::initMarking(adaptInfo);
151 
152  double esThetaP = std::pow(esTheta_, this->p_);
153  double esThetaCP = std::pow(esThetaC_, this->p_);
154  double epsP = std::pow(adaptInfo.spaceTolerance(this->component_), this->p_);
155 
156  int nLeaves = (this->grid_->leafGridView()).size(0);
157 #if AMDIS_HAS_PARALLEL
158  Dune::Communication<>{}.sum(nLeaves, 1);
159 #endif
160 
161  this->markRLimit_ = esThetaP * epsP / nLeaves;
162  this->markCLimit_ = esThetaCP * epsP / nLeaves;
163 
164  msg("start mark_limits: {}, {}; nt = {}", this->markRLimit_, this->markCLimit_, nLeaves);
165 }
166 
167 
168 template <class Grid>
170 {
171  Super::initMarking(adaptInfo);
172 
173  if (!this->coarsenAllowed_ && !this->refineAllowed_)
174  return 0;
175 
176  gersSum_ = 0.0;
177 
178  double LTheta = std::pow(1.0 - gersThetaStar_, this->p_);
179  double epsP = std::pow(adaptInfo.spaceTolerance(this->component_), this->p_);
180 
181  if (this->estSum_ < oldErrSum_) {
182  double improv = this->estSum_ / oldErrSum_;
183  double wanted = 0.8 * epsP / this->estSum_;
184  double redfac = std::min((1.0 - wanted) / (1.0 - improv), 1.0);
185  redfac = std::max(redfac, 0.0);
186 
187  if (redfac < 1.0) {
188  LTheta *= redfac;
189  msg("GERS: use extrapolated theta_star = {}", std::pow(LTheta, 1.0 / this->p_));
190  }
191  }
192 
193  oldErrSum_ = this->estSum_;
194  double GERSGamma = 1.0;
195 
196  if (this->refineAllowed_) {
197  if (LTheta > 0) {
198  do {
199  gersSum_ = 0.0;
200  GERSGamma -= gersNu_;
201  this->markRLimit_ = GERSGamma * this->estMax_;
202 
203  for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
204  markElementForRefinement(adaptInfo, elem);
205 
206  } while((GERSGamma > 0) && (gersSum_ < LTheta * this->estSum_));
207  }
208 
209  msg("GERS refinement with gamma = {}", GERSGamma);
210  }
211 
212  if (this->coarsenAllowed_) {
213  GERSGamma = 0.3;
214  LTheta = gersThetaC_ * epsP;
215 
216  do {
217  gersSum_ = 0.0;
218  GERSGamma -= gersNu_;
219  this->markCLimit_ = GERSGamma * this->estMax_;
220 
221  for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
222  markElementForCoarsening(adaptInfo, elem);
223 
224  msg("coarse loop: gamma = {}, sum = {}, limit = {}", GERSGamma, gersSum_, LTheta);
225  } while(gersSum_ > LTheta);
226 
227  msg("GERS coarsening with gamma = {}", GERSGamma);
228  }
229 
230  Super::finishMarking(adaptInfo);
231 
232  Flag markFlag;
233  if (this->elMarkRefine_)
234  markFlag = 1;
235  if (this->elMarkCoarsen_)
236  markFlag |= 2;
237 
238  return markFlag;
239 }
240 
241 
242 template <class Grid>
243 void GERSMarker<Grid>::markElementForRefinement(AdaptInfo& /*adaptInfo*/, const Element& elem)
244 {
245  double lError = this->est_[this->grid_->leafIndexSet().index(elem)];
246 
247  if (lError > this->markRLimit_) {
248  gersSum_ += lError;
249  this->mark(elem, 1);
250  }
251 }
252 
253 
254 template <class Grid>
255 void GERSMarker<Grid>::markElementForCoarsening(AdaptInfo& /*adaptInfo*/, const Element& elem)
256 {
257  double lError = this->est_[this->grid_->leafIndexSet().index(elem)];
258 
259  if (this->grid_->getMark(elem) <= 0) {
260  if (lError <= this->markCLimit_) {
261  gersSum_ += lError;
262  this->mark(elem, -1);
263  } else {
264  this->mark(elem, 0);
265  }
266  }
267 }
268 
269 
270 template <class Grid, class PreGridFct>
272 {
273  test_exit(bool(this->grid_), "No grid!");
274 
275  Super::initMarking(adaptInfo);
276 
277  if (!this->coarsenAllowed_ && !this->refineAllowed_)
278  return 0;
279 
280  auto localFct = localFunction(gridFct_);
281 
282  for (auto const& e : Dune::elements(this->grid_->leafGridView())) {
283  localFct.bind(e);
284  int currentLevel = e.level();
285  auto refElem = Dune::referenceElement<typename Grid::ctype,Grid::dimension>(e.type());
286 
287  // evaluate in the center of the element
288  int targetLevel = int(std::round(localFct(refElem.position(0,0))));
289 
290  int m = ((((targetLevel > currentLevel) && (currentLevel < this->maxRefineLevel_))
291  || (currentLevel < this->minRefineLevel_))
292  && this->refineAllowed_)
293  - ((((targetLevel < currentLevel) && (currentLevel > this->minRefineLevel_))
294  || (currentLevel > this->maxRefineLevel_))
295  && this->coarsenAllowed_);
296  this->mark(e, m);
297  localFct.unbind();
298  }
299 
300  Super::finishMarking(adaptInfo);
301 
302  Flag markFlag;
303  if (this->elMarkRefine_)
304  markFlag = 1;
305  if (this->elMarkCoarsen_)
306  markFlag |= 2;
307 
308  return markFlag;
309 }
310 
311 } // end namespace AMDiS
void markElementForCoarsening(AdaptInfo &adaptInfo, Element const &elem)
Coarsening marking function.
Definition: Marker.inc.hpp:255
The Flag class encapsulates flags which represents simple information. Used e.g. while mesh traversal...
Definition: Flag.hpp:13
bool isRefinementAllowed(Key key) const
Returns whether coarsening is allowed or not.
Definition: AdaptInfo.hpp:582
virtual void finishMarking(AdaptInfo &adaptInfo)
Called after traversal.
Definition: Marker.inc.hpp:40
double estMax(Key key) const
Returns est_max.
Definition: AdaptInfo.hpp:357
void initMarking(AdaptInfo &adaptInfo) override
Can be used by sub classes. Called before traversal.
Definition: Marker.inc.hpp:74
double spaceTolerance(Key key) const
Returns spaceTolerance.
Definition: AdaptInfo.hpp:404
void markElement(AdaptInfo &adaptInfo, Element const &elem) override
Marks one element.
Definition: Marker.inc.hpp:85
void markElementForRefinement(AdaptInfo &adaptInfo, Element const &elem)
Refinement marking function.
Definition: Marker.inc.hpp:243
Definition: AdaptBase.hpp:6
Flag markGrid(AdaptInfo &adaptInfo) override
Marking of the whole grid.
Definition: Marker.inc.hpp:169
static std::optional< T > get(std::string const &key)
Get parameter-values from parameter-tree.
Definition: Initfile.hpp:25
bool isCoarseningAllowed(Key key) const
Returns whether coarsening is allowed or not.
Definition: AdaptInfo.hpp:570
void mark(Element const &elem, int newMark)
Definition: Marker.inc.hpp:6
virtual Flag markGrid(AdaptInfo &adaptInfo)
Marking of the whole grid.
Definition: Marker.inc.hpp:50
virtual void initMarking(AdaptInfo &adaptInfo)
Called before traversal.
Definition: Marker.inc.hpp:32
Holds adapt parameters and infos about the problem.
Definition: AdaptInfo.hpp:25
void initMarking(AdaptInfo &adaptInfo) override
Can be used by sub classes. Called before traversal.
Definition: Marker.inc.hpp:148
void initMarking(AdaptInfo &adaptInfo) override
Can be used by sub classes. Called before traversal.
Definition: Marker.inc.hpp:132
Flag markGrid(AdaptInfo &adaptInfo) override
Definition: Marker.inc.hpp:271
double estSum(Key key) const
Returns est_sum.
Definition: AdaptInfo.hpp:333
static std::unique_ptr< EstimatorMarker< Grid > > createMarker(std::string const &name, std::string const &component, Estimates const &est, std::shared_ptr< Grid > const &grid)
Creates a scalar marker depending on the strategy set in parameters.
Definition: Marker.inc.hpp:102