Marker.inc.hpp 6.82 KB
Newer Older
1
2
3
4
5
6
7
8
// TODO: Cleanup of copied comments

#include <amdis/common/Math.hpp>

namespace AMDiS {

  using std::pow;

9
  template <class Traits>
10
  std::shared_ptr<Marker<Traits> > Marker<Traits>::createMarker(std::string name, std::string component, const EstType& est, std::shared_ptr<Grid> grid)
11
12
13
14
  {
    int strategy = 0;
    Parameters::get(name + "->strategy", strategy);

15
    std::shared_ptr<Marker<Traits> > marker;
16
17
18
19
20

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

    return marker;
  }


  
45
46
  template <class Traits>
  void Marker<Traits>::initMarking(AdaptInfo& adaptInfo)
47
48
49
  {
    elMarkRefine = 0;
    elMarkCoarsen = 0;
50
51
52
53
    estSum = pow(adaptInfo.getEstSum(component), p);
    estMax = adaptInfo.getEstMax(component);
    refineAllowed = adaptInfo.isRefinementAllowed(component);
    coarsenAllowed = adaptInfo.isCoarseningAllowed(component);
54
55
56
57
  }


  
58
59
  template <class Traits>
  void Marker<Traits>::finishMarking(AdaptInfo& adaptInfo)
60
  {
61
62
    msg(elMarkRefine, " elements marked for refinement");
    msg(elMarkCoarsen, " elements marked for coarsening");
63
64
65
66
  }


  
67
68
  template <class Traits>
  void Marker<Traits>::markElement(AdaptInfo& adaptInfo, const Element& elem)
69
70
  {
    const auto& index = grid->leafIndexSet().index(elem);
71
    double lError = est[index];
72

73
74
75
76
    if (lError > markRLimit && refineAllowed
        && (maxRefineLevel == -1 || elem.level() < maxRefineLevel)) {
      this->mark(elem, 1);
      
77
    } else {
78
79
80
81
82
      if (lError <= markCLimit && coarsenAllowed
          && (minRefineLevel == -1 || elem.level() > minRefineLevel))
        if (/*!elem->getElementData()->getElementData(COARSENABLE) ||*/
            lError /*+ elem->getCoarseningEstimation(index)*/ <= markCLimit)
          this->mark(elem, -1);
83
84
85
86
87
    }
  }


  
88
89
  template <class Traits>
  Flag Marker<Traits>::markGrid(AdaptInfo& adaptInfo)
90
  {
91
92
    if (!(grid))
      error_exit("No grid!");
93
94
95

    initMarking(adaptInfo);

96
97
    if (!coarsenAllowed &&
        !refineAllowed) {
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
      return 0;
    }
    for (const auto& elem : Dune::elements(grid->leafGridView())) {
      markElement(adaptInfo, elem);
    }

    finishMarking(adaptInfo);

    Flag markFlag;
    if (elMarkRefine)
      markFlag = 1;
    if (elMarkCoarsen)
      markFlag |= 2;

    return markFlag;
  }


116
117
  template <class Traits>
  void MSMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
118
  {
119
    Marker<Traits>::initMarking(adaptInfo);
120
121
122
123

    double MSGammaP = pow(MSGamma, this->p);
    double MSGammaCP = pow(MSGammaC, this->p);

124
125
    this->markRLimit = MSGammaP * adaptInfo.getEstMax(this->component);
    this->markCLimit = MSGammaCP * adaptInfo.getEstMax(this->component);
126

127
    msg("start max_est: ", adaptInfo.getEstMax(this->component), ", mark_limits: ", this->markRLimit, ", " , this->markCLimit);
128
129
130
  }


131
132
  template <class Traits>
  void ESMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
133
  {
134
    Marker<Traits>::initMarking(adaptInfo);
135
136
137

    double ESThetaP = pow(ESTheta, this->p);
    double ESThetaCP = pow(ESThetaC, this->p);
138
    double epsP = pow(adaptInfo.getSpaceTolerance(this->component), this->p);
139
140
141
142
143
144
145
146
147

    int nLeaves = (this->grid->leafGridView()).size(0);
/*#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    Parallel::mpi::globalAdd(nLeaves);
#endif*/

    this->markRLimit = ESThetaP * epsP / nLeaves;
    this->markCLimit = ESThetaCP * epsP / nLeaves;

148
    msg("start mark_limits: ", this->markRLimit, ", " , this->markCLimit, "; nt = ", nLeaves);
149
150
151
  }


152
153
  template <class Traits>
  Flag GERSMarker<Traits>::markGrid(AdaptInfo& adaptInfo)
154
  {
155
    Marker<Traits>::initMarking(adaptInfo);
156

157
158
    if (!this->coarsenAllowed &&
        !this->refineAllowed)
159
160
161
162
163
      return 0;

    GERSSum = 0.0;

    double LTheta = pow(1.0 - GERSThetaStar, this->p);
164
    double epsP = pow(adaptInfo.getSpaceTolerance(this->component), this->p);
165
166
167
168
169
170
171
172
173

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

      if (redfac < 1.0) {
        LTheta *= redfac;
174
        msg("GERS: use extrapolated theta_star = ", pow(LTheta, 1.0 / this->p));
175
176
177
178
179
180
      }
    }

    oldErrSum = this->estSum;
    double GERSGamma = 1.0;

181
    if (this->refineAllowed) {
182
183
184
185
186
187
188
189
190
191
192
193
      if (LTheta > 0) {
        do {
          GERSSum = 0.0;
          GERSGamma -= GERSNu;
          this->markRLimit = GERSGamma * this->estMax;

          for (const auto& elem : Dune::elements(this->grid->leafGridView()))
            markElementForRefinement(adaptInfo, elem);

        } while((GERSGamma > 0) && (GERSSum < LTheta * this->estSum));
      }

194
      msg("GERS refinement with gamma = ", GERSGamma);
195
196
    }

197
    if (this->coarsenAllowed) {
198
199
200
201
202
203
204
205
206
207
208
      GERSGamma = 0.3;
      LTheta = GERSThetaC * epsP;

      do {
        GERSSum = 0.0;
        GERSGamma -= GERSNu;
        this->markCLimit = GERSGamma * this->estMax;

        for (const auto& elem : Dune::elements(this->grid->leafGridView()))
          markElementForCoarsening(adaptInfo, elem);

209
        msg("coarse loop: gamma = ", GERSGamma, ", sum = ", GERSSum, ", limit = ", LTheta);
210
211
      } while(GERSSum > LTheta);

212
      msg("GERS coarsening with gamma = ", GERSGamma);
213
214
    }

215
    Marker<Traits>::finishMarking(adaptInfo);
216
217
218
219
220
221
222
223
224
225
226

    Flag markFlag;
    if (this->elMarkRefine)
      markFlag = 1;
    if (this->elMarkCoarsen)
      markFlag |= 2;

    return(markFlag);
  }


227
228
  template <class Traits>
  void GERSMarker<Traits>::markElementForRefinement(AdaptInfo& adaptInfo, const Element& elem)
229
  {
230
    double lError = (this->est)[(this->grid->leafIndexSet()).index(elem)];
231
232
233
234
235
236
237
238

    if (lError > this->markRLimit) {
      GERSSum += lError;
      this->mark(elem, 1);
    }
  }


239
240
  template <class Traits>
  void GERSMarker<Traits>::markElementForCoarsening(AdaptInfo& adaptInfo, const Element& elem)
241
  {
242
    double lError = (this->est)[(this->grid->leafIndexSet()).index(elem)];
243
244
245

    if ((this->grid)->getMark(elem) <= 0) {
/*      if (elem->getElementData()->getElementData(COARSENABLE))*/
246
/*        lError += elem->getCoarseningEstimation(index);*/
247
248
249
250
251
252
253
254
255
256
257
258

      if (lError <= this->markCLimit) {
        GERSSum += lError;
        this->mark(elem, -1);
      } else {
        this->mark(elem, 0);
      }
    }
  }


}