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

namespace AMDiS {

5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
template <class Traits>
std::shared_ptr<Marker<Traits> > Marker<Traits>::createMarker(std::string name, std::string component, Estimates const& est, std::shared_ptr<Grid> const& grid)
{
  int strategy = 0;
  Parameters::get(name + "->strategy", strategy);

  std::shared_ptr<Marker<Traits> > marker;

  switch (strategy) {
  case 0:  // no refinement/coarsening
    break;
  case 1:
    msg("Creating global refinement (GR) marker");
    marker = std::make_shared<GRMarker<Traits> >(name, component, est, grid);
    break;
  case 2:
    msg("Creating maximum strategy (MS) marker");
    marker = std::make_shared<MSMarker<Traits> >(name, component, est, grid);
    break;
  case 3:
    msg("Creating equidistribution strategy (ES) marker");
    marker = std::make_shared<ESMarker<Traits> >(name, component, est, grid);
    break;
  case 4:
    msg("Creating guaranteed error reduction strategy (GERS) marker");
    marker = std::make_shared<GERSMarker<Traits> >(name, component, est, grid);
    break;
  default:
    error_exit("invalid strategy");
34
35
  }

36
37
  return marker;
}
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
55
56
57
template <class Traits>
void Marker<Traits>::finishMarking(AdaptInfo& adaptInfo)
{
  msg(elMarkRefine_, " elements marked for refinement");
  msg(elMarkCoarsen_, " elements marked for coarsening");
}
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
113
114
  this->markRLimit_ = msGammaP * adaptInfo.getEstMax(this->component_);
  this->markCLimit_ = msGammaCP * adaptInfo.getEstMax(this->component_);

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


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

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

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

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

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


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

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

146
  gersSum_ = 0.0;
147

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

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

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

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

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

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

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

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

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

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

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

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

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

200
  Super::finishMarking(adaptInfo);
201

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

208
209
  return markFlag;
}
210
211


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

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


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

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


240
} // end namespace AMDiS