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

namespace AMDiS {

5
template <class Traits>
6
7
8
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)
9
10
11
12
13
14
15
16
17
{
  int strategy = 0;
  Parameters::get(name + "->strategy", strategy);

  switch (strategy) {
  case 0:  // no refinement/coarsening
    break;
  case 1:
    msg("Creating global refinement (GR) marker");
18
    return std::make_unique<GRMarker<Traits> >(name, component, est, grid);
19
20
21
    break;
  case 2:
    msg("Creating maximum strategy (MS) marker");
22
    return std::make_unique<MSMarker<Traits> >(name, component, est, grid);
23
24
25
    break;
  case 3:
    msg("Creating equidistribution strategy (ES) marker");
26
    return std::make_unique<ESMarker<Traits> >(name, component, est, grid);
27
28
29
    break;
  case 4:
    msg("Creating guaranteed error reduction strategy (GERS) marker");
30
    return std::make_unique<GERSMarker<Traits> >(name, component, est, grid);
31
32
33
    break;
  default:
    error_exit("invalid strategy");
34
35
  }

36
  return {};
37
}
38
39


40
41
42
43
44
45
46
47
48
49
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_);
}
50
51


52
53
54
template <class Traits>
void Marker<Traits>::finishMarking(AdaptInfo& adaptInfo)
{
55
56
  msg("{} elements marked for refinement", elMarkRefine_);
  msg("{} elements marked for coarsening", elMarkCoarsen_);
57
}
58
59


60
61
62
63
64
65
66
67
68
69
70
71
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);
72
  }
73
}
74
75


76
77
78
79
80
template <class Traits>
Flag Marker<Traits>::markGrid(AdaptInfo& adaptInfo)
{
  if (!grid_)
    error_exit("No grid!");
81

82
  initMarking(adaptInfo);
83

84
85
  if (!coarsenAllowed_ && !refineAllowed_)
    return 0;
86

87
88
  for (const auto& elem : Dune::elements(grid_->leafGridView()))
    markElement(adaptInfo, elem);
89

90
  finishMarking(adaptInfo);
91

92
93
94
95
96
  Flag markFlag;
  if (elMarkRefine_)
    markFlag = 1;
  if (elMarkCoarsen_)
    markFlag |= 2;
97

98
99
  return markFlag;
}
100
101


102
103
104
105
template <class Traits>
void MSMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
{
  Super::initMarking(adaptInfo);
106

107
108
  double msGammaP = std::pow(msGamma_, this->p_);
  double msGammaCP = std::pow(msGammaC_, this->p_);
109

110
111
112
  this->markRLimit_ = msGammaP * adaptInfo.getEstMax(this->component_);
  this->markCLimit_ = msGammaCP * adaptInfo.getEstMax(this->component_);

113
114
  msg("start max_est: {}, mark_limits: {}, {}",
    adaptInfo.getEstMax(this->component_), this->markRLimit_ , this->markCLimit_);
115
}
116
117


118
119
120
121
template <class Traits>
void ESMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
{
  Super::initMarking(adaptInfo);
122

123
124
125
  double esThetaP = std::pow(esTheta_, this->p_);
  double esThetaCP = std::pow(esThetaC_, this->p_);
  double epsP = std::pow(adaptInfo.getSpaceTolerance(this->component_), this->p_);
126

127
  int nLeaves = (this->grid_->leafGridView()).size(0);
128
/*#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
129
  Parallel::mpi::globalAdd(nLeaves);
130
131
#endif*/

132
133
  this->markRLimit_ = esThetaP * epsP / nLeaves;
  this->markCLimit_ = esThetaCP * epsP / nLeaves;
134

135
  msg("start mark_limits: {}, {}; nt = {}", this->markRLimit_, this->markCLimit_, nLeaves);
136
}
137
138


139
140
141
142
template <class Traits>
Flag GERSMarker<Traits>::markGrid(AdaptInfo& adaptInfo)
{
  Super::initMarking(adaptInfo);
143

144
145
  if (!this->coarsenAllowed_ && !this->refineAllowed_)
    return 0;
146

147
  gersSum_ = 0.0;
148

149
150
  double LTheta = std::pow(1.0 - gersThetaStar_, this->p_);
  double epsP = std::pow(adaptInfo.getSpaceTolerance(this->component_), this->p_);
151

152
153
154
155
156
  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);
157

158
159
    if (redfac < 1.0) {
      LTheta *= redfac;
160
      msg("GERS: use extrapolated theta_star = {}", std::pow(LTheta, 1.0 / this->p_));
161
    }
162
  }
163

164
165
  oldErrSum_ = this->estSum_;
  double GERSGamma = 1.0;
166

167
168
169
170
171
172
  if (this->refineAllowed_) {
    if (LTheta > 0) {
      do {
        gersSum_ = 0.0;
        GERSGamma -= gersNu_;
        this->markRLimit_ = GERSGamma * this->estMax_;
173

174
175
        for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
          markElementForRefinement(adaptInfo, elem);
176

177
      } while((GERSGamma > 0) && (gersSum_ < LTheta * this->estSum_));
178
179
    }

180
    msg("GERS refinement with gamma = {}", GERSGamma);
181
  }
182

183
184
185
  if (this->coarsenAllowed_) {
    GERSGamma = 0.3;
    LTheta = gersThetaC_ * epsP;
186

187
188
189
190
    do {
      gersSum_ = 0.0;
      GERSGamma -= gersNu_;
      this->markCLimit_ = GERSGamma * this->estMax_;
191

192
193
      for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
        markElementForCoarsening(adaptInfo, elem);
194

195
      msg("coarse loop: gamma = {}, sum = {}, limit = {}", GERSGamma, gersSum_, LTheta);
196
197
    } while(gersSum_ > LTheta);

198
    msg("GERS coarsening with gamma = {}", GERSGamma);
199
  }
200

201
  Super::finishMarking(adaptInfo);
202

203
204
205
206
207
  Flag markFlag;
  if (this->elMarkRefine_)
    markFlag = 1;
  if (this->elMarkCoarsen_)
    markFlag |= 2;
208

209
210
  return markFlag;
}
211
212


213
214
215
216
template <class Traits>
void GERSMarker<Traits>::markElementForRefinement(AdaptInfo& adaptInfo, const Element& elem)
{
  double lError = this->est_[this->grid_->leafIndexSet().index(elem)];
217

218
219
220
  if (lError > this->markRLimit_) {
    gersSum_ += lError;
    this->mark(elem, 1);
221
  }
222
}
223
224


225
226
227
228
template <class Traits>
void GERSMarker<Traits>::markElementForCoarsening(AdaptInfo& adaptInfo, const Element& elem)
{
  double lError = this->est_[this->grid_->leafIndexSet().index(elem)];
229

230
231
232
233
234
235
  if (this->grid_->getMark(elem) <= 0) {
    if (lError <= this->markCLimit_) {
      gersSum_ += lError;
      this->mark(elem, -1);
    } else {
      this->mark(elem, 0);
236
237
    }
  }
238
}
239
240


241
} // end namespace AMDiS