Marker.inc.hpp 7.84 KB
Newer Older
1
#include <cmath>
2
3
4

namespace AMDiS {

Praetorius, Simon's avatar
Praetorius, Simon committed
5
6
template <class Grid>
void Marker<Grid>::mark(Element const& elem, int newMark)
7
{
Praetorius, Simon's avatar
Praetorius, Simon committed
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
  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");
    }
27
  }
28
}
29
30


Praetorius, Simon's avatar
Praetorius, Simon committed
31
32
template <class Grid>
void Marker<Grid>::initMarking(AdaptInfo& adaptInfo)
33
{
Praetorius, Simon's avatar
Praetorius, Simon committed
34
35
  this->elMarkRefine_ = 0;
  this->elMarkCoarsen_ = 0;
36
}
37
38


Praetorius, Simon's avatar
Praetorius, Simon committed
39
40
template <class Grid>
void Marker<Grid>::finishMarking(AdaptInfo& adaptInfo)
41
{
Praetorius, Simon's avatar
Praetorius, Simon committed
42
43
44
  if (info_ > 0) {
    msg("{} elements marked for refinement", elMarkRefine_);
    msg("{} elements marked for coarsening", elMarkCoarsen_);
45
  }
46
}
47
48


Praetorius, Simon's avatar
Praetorius, Simon committed
49
50
template <class Grid>
Flag Marker<Grid>::markGrid(AdaptInfo& adaptInfo)
51
{
Praetorius, Simon's avatar
Praetorius, Simon committed
52
  test_exit(bool(this->grid_), "No grid!");
53

54
  initMarking(adaptInfo);
55

Praetorius, Simon's avatar
Praetorius, Simon committed
56
  if (!this->coarsenAllowed_ && !this->refineAllowed_)
57
    return 0;
58

Praetorius, Simon's avatar
Praetorius, Simon committed
59
  for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
60
    markElement(adaptInfo, elem);
61

62
  finishMarking(adaptInfo);
63

64
  Flag markFlag;
Praetorius, Simon's avatar
Praetorius, Simon committed
65
  if (this->elMarkRefine_)
66
    markFlag = 1;
Praetorius, Simon's avatar
Praetorius, Simon committed
67
  if (this->elMarkCoarsen_)
68
    markFlag |= 2;
69

70
71
  return markFlag;
}
72

Praetorius, Simon's avatar
Praetorius, Simon committed
73
74
75
76
template <class Grid>
void EstimatorMarker<Grid>::initMarking(AdaptInfo& adaptInfo)
{
  Super::initMarking(adaptInfo);
Praetorius, Simon's avatar
Praetorius, Simon committed
77
78
  estSum_ = std::pow(adaptInfo.estSum(component_), p_);
  estMax_ = adaptInfo.estMax(component_);
Praetorius, Simon's avatar
Praetorius, Simon committed
79
80
81
82
  this->refineAllowed_ = adaptInfo.isRefinementAllowed(component_);
  this->coarsenAllowed_ = adaptInfo.isCoarseningAllowed(component_);
}

83

Praetorius, Simon's avatar
Praetorius, Simon committed
84
85
86
87
88
89
90
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_
91
      && elem.level() < this->maxRefineLevel_) {
Praetorius, Simon's avatar
Praetorius, Simon committed
92
93
    this->mark(elem, 1);
  } else  if (lError <= markCLimit_ && this->coarsenAllowed_
94
              && elem.level() > this->minRefineLevel_) {
Praetorius, Simon's avatar
Praetorius, Simon committed
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
    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)
133
134
{
  Super::initMarking(adaptInfo);
135

136
137
  double msGammaP = std::pow(msGamma_, this->p_);
  double msGammaCP = std::pow(msGammaC_, this->p_);
138

Praetorius, Simon's avatar
Praetorius, Simon committed
139
140
  this->markRLimit_ = msGammaP * adaptInfo.estMax(this->component_);
  this->markCLimit_ = msGammaCP * adaptInfo.estMax(this->component_);
141

142
  msg("start max_est: {}, mark_limits: {}, {}",
Praetorius, Simon's avatar
Praetorius, Simon committed
143
    adaptInfo.estMax(this->component_), this->markRLimit_ , this->markCLimit_);
144
}
145
146


Praetorius, Simon's avatar
Praetorius, Simon committed
147
148
template <class Grid>
void ESMarker<Grid>::initMarking(AdaptInfo& adaptInfo)
149
150
{
  Super::initMarking(adaptInfo);
151

152
153
  double esThetaP = std::pow(esTheta_, this->p_);
  double esThetaCP = std::pow(esThetaC_, this->p_);
Praetorius, Simon's avatar
Praetorius, Simon committed
154
  double epsP = std::pow(adaptInfo.spaceTolerance(this->component_), this->p_);
155

156
  int nLeaves = (this->grid_->leafGridView()).size(0);
Praetorius, Simon's avatar
Praetorius, Simon committed
157
158
159
#if AMDIS_HAS_PARALLEL
  Dune::CollectiveCommunication<>{}.sum(nLeaves, 1);
#endif
160

161
162
  this->markRLimit_ = esThetaP * epsP / nLeaves;
  this->markCLimit_ = esThetaCP * epsP / nLeaves;
163

164
  msg("start mark_limits: {}, {}; nt = {}", this->markRLimit_, this->markCLimit_, nLeaves);
165
}
166
167


Praetorius, Simon's avatar
Praetorius, Simon committed
168
169
template <class Grid>
Flag GERSMarker<Grid>::markGrid(AdaptInfo& adaptInfo)
170
171
{
  Super::initMarking(adaptInfo);
172

173
174
  if (!this->coarsenAllowed_ && !this->refineAllowed_)
    return 0;
175

176
  gersSum_ = 0.0;
177

178
  double LTheta = std::pow(1.0 - gersThetaStar_, this->p_);
Praetorius, Simon's avatar
Praetorius, Simon committed
179
  double epsP = std::pow(adaptInfo.spaceTolerance(this->component_), this->p_);
180

181
182
183
184
185
  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);
186

187
188
    if (redfac < 1.0) {
      LTheta *= redfac;
189
      msg("GERS: use extrapolated theta_star = {}", std::pow(LTheta, 1.0 / this->p_));
190
    }
191
  }
192

193
194
  oldErrSum_ = this->estSum_;
  double GERSGamma = 1.0;
195

196
197
198
199
200
201
  if (this->refineAllowed_) {
    if (LTheta > 0) {
      do {
        gersSum_ = 0.0;
        GERSGamma -= gersNu_;
        this->markRLimit_ = GERSGamma * this->estMax_;
202

203
204
        for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
          markElementForRefinement(adaptInfo, elem);
205

206
      } while((GERSGamma > 0) && (gersSum_ < LTheta * this->estSum_));
207
208
    }

209
    msg("GERS refinement with gamma = {}", GERSGamma);
210
  }
211

212
213
214
  if (this->coarsenAllowed_) {
    GERSGamma = 0.3;
    LTheta = gersThetaC_ * epsP;
215

216
217
218
219
    do {
      gersSum_ = 0.0;
      GERSGamma -= gersNu_;
      this->markCLimit_ = GERSGamma * this->estMax_;
220

221
222
      for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
        markElementForCoarsening(adaptInfo, elem);
223

224
      msg("coarse loop: gamma = {}, sum = {}, limit = {}", GERSGamma, gersSum_, LTheta);
225
226
    } while(gersSum_ > LTheta);

227
    msg("GERS coarsening with gamma = {}", GERSGamma);
228
  }
229

230
  Super::finishMarking(adaptInfo);
231

232
233
234
235
236
  Flag markFlag;
  if (this->elMarkRefine_)
    markFlag = 1;
  if (this->elMarkCoarsen_)
    markFlag |= 2;
237

238
239
  return markFlag;
}
240
241


Praetorius, Simon's avatar
Praetorius, Simon committed
242
243
template <class Grid>
void GERSMarker<Grid>::markElementForRefinement(AdaptInfo& adaptInfo, const Element& elem)
244
245
{
  double lError = this->est_[this->grid_->leafIndexSet().index(elem)];
246

247
248
249
  if (lError > this->markRLimit_) {
    gersSum_ += lError;
    this->mark(elem, 1);
250
  }
251
}
252
253


Praetorius, Simon's avatar
Praetorius, Simon committed
254
255
template <class Grid>
void GERSMarker<Grid>::markElementForCoarsening(AdaptInfo& adaptInfo, const Element& elem)
256
257
{
  double lError = this->est_[this->grid_->leafIndexSet().index(elem)];
258

259
260
261
262
263
264
  if (this->grid_->getMark(elem) <= 0) {
    if (lError <= this->markCLimit_) {
      gersSum_ += lError;
      this->mark(elem, -1);
    } else {
      this->mark(elem, 0);
265
266
    }
  }
267
}
268
269


Praetorius, Simon's avatar
Praetorius, Simon committed
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
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 codim = ref.dimension;

289
290
    // evaluate in the center of the element
    int targetLevel = int(std::round(localFct(ref.position(0,0))));
Praetorius, Simon's avatar
Praetorius, Simon committed
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312

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

313
} // end namespace AMDiS