3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
gridadapt.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
23#ifndef DUMUX_GRIDADAPT_HH
24#define DUMUX_GRIDADAPT_HH
25
26#include <unordered_map>
27#include <dune/grid/common/partitionset.hh>
28
29#include "properties.hh"
31
32namespace Dumux
33{
34
45template<class TypeTag, bool adaptive>
47{
51
52 using Grid = typename GridView::Grid;
53 using LeafGridView = typename Grid::LeafGridView;
54 using Element = typename Grid::template Codim<0>::Entity;
55
58 using AdaptionInitializationIndicator = GetPropType<TypeTag, Properties::AdaptionInitializationIndicator>;
59
60 using IdSet = typename Grid::Traits::LocalIdSet;
61 using IdType = typename IdSet::IdType;
62
63public:
68 GridAdapt (Problem& problem)
69 : problem_(problem), adaptionIndicator_(problem), marked_(0), coarsened_(0)
70 {
71 levelMin_ = getParam<int>("GridAdapt.MinLevel");
72 levelMax_ = getParam<int>("GridAdapt.MaxLevel");
73 adaptationInterval_ = getParam<int>("GridAdapt.AdaptionInterval", 1);
74
75 if (levelMin_ < 0)
76 Dune::dgrave << __FILE__<< ":" <<__LINE__
77 << " : Dune cannot coarsen to gridlevels smaller 0! "<< std::endl;
78 }
79
87 void init()
88 {
89 adaptionIndicator_.init();
90
91 if (!getParam<bool>("GridAdapt.EnableInitializationIndicator"))
92 return;
93
94 AdaptionInitializationIndicator adaptionInitIndicator(problem_, adaptionIndicator_);
95
96 int maxIter = 2*levelMax_;
97 int iter = 0;
98 while (iter <= maxIter)
99 {
100 adaptGrid(adaptionInitIndicator);
101
102 if (!wasAdapted())
103 {
104 break;
105 }
106
107 int shouldInitialize = adaptionInitIndicator.initializeModel();
108 if (problem_.grid().comm().max(shouldInitialize))
109 problem_.model().initialize();
110
111 iter++;
112 }
113 }
114
130 {
131 adaptGrid(adaptionIndicator_) ;
132 }
133
150 template<class Indicator>
151 void adaptGrid(Indicator& indicator)
152 {
153 // reset internal counter for marked elements
154 marked_ = coarsened_ = 0;
155
156 // check for adaption interval: Adapt only at certain time step indices
157 if (problem_.timeManager().timeStepIndex() % adaptationInterval_ != 0)
158 return;
159
160 /**** 1) determine refining parameter if standard is used ***/
161 // if not, the indicatorVector and refinement Bounds have to
162 // specified by the problem through setIndicator()
163 indicator.calculateIndicator();
164
165 /**** 2) mark elements according to indicator *********/
166 markElements(indicator);
167
168 // abort if nothing in grid is marked
169 int sumMarked = problem_.grid().comm().sum(marked_);
170 int sumCoarsened = problem_.grid().comm().sum(coarsened_);
171 if (sumMarked == 0 && sumCoarsened == 0)
172 return;
173 else
174 Dune::dinfo << marked_ << " cells have been marked_ to be refined, "
175 << coarsened_ << " to be coarsened." << std::endl;
176
177 /**** 2b) Do pre-adaption step *****/
178 problem_.grid().preAdapt();
179 problem_.preAdapt();
180
181 /**** 3) Put primary variables in a map *********/
182 problem_.variables().storePrimVars(problem_);
183
184 /**** 4) Adapt Grid and size of variable vectors *****/
185 problem_.grid().adapt();
186
187 // forceRefineRatio(1);
188
189 // update mapper to new cell indices
190 problem_.variables().elementMapper().update();
191
192 // adapt size of vectors
193 problem_.variables().adaptVariableSize(problem_.variables().elementMapper().size());
194
195 /**** 5) (Re-)construct primary variables to new grid **/
196 problem_.variables().reconstructPrimVars(problem_);
197
198 // delete markers in grid
199 problem_.grid().postAdapt();
200
201 // call method in problem for potential output etc.
202 problem_.postAdapt();
203
204 return;
205 }
206
211 template<class Indicator>
212 void markElements(Indicator& indicator)
213 {
214 using CoarsenMarkerType = std::unordered_map<IdType, IdType>;
215 CoarsenMarkerType coarsenMarker;
216 const IdSet& idSet(problem_.grid().localIdSet());
217
218 for (const auto& element : elements(problem_.gridView()))
219 {
220 // only mark non-ghost elements
221 if (element.partitionType() == Dune::GhostEntity)
222 continue;
223
224 // refine?
225 if (indicator.refine(element) && element.level() < levelMax_)
226 {
227 problem_.grid().mark( 1, element);
228 ++marked_;
229
230 // this also refines the neighbor elements
231 checkNeighborsRefine_(element);
232 }
233 if (indicator.coarsen(element) && element.hasFather())
234 {
235 IdType idx = idSet.id(element.father());
236 auto it = coarsenMarker.find(idx);
237 if (it != coarsenMarker.end())
238 {
239 ++it->second;
240 }
241 else
242 {
243 coarsenMarker[idx] = 1;
244 }
245 }
246 }
247 // coarsen
248 for (const auto& element : elements(problem_.gridView()))
249 {
250 // only mark non-ghost elements
251 if (element.partitionType() == Dune::GhostEntity)
252 continue;
253
254 if (indicator.coarsen(element) && element.level() > levelMin_)
255 {
256 IdType idx = idSet.id(element.father());
257 auto it = coarsenMarker.find(idx);
258 if (it != coarsenMarker.end())
259 {
260 if (problem_.grid().getMark(element) == 0
261 && it->second == element.geometry().corners())
262 {
263 // check if coarsening is possible
264 bool coarsenPossible = true;
265 for(const auto& intersection : intersections(problem_.gridView(), element))
266 {
267 if(intersection.neighbor())
268 {
269 auto outside = intersection.outside();
270 if ((problem_.grid().getMark(outside) > 0)
271 || outside.level() > element.level())
272 {
273 coarsenPossible = false;
274 }
275 }
276 }
277
278 if(coarsenPossible)
279 {
280 problem_.grid().mark( -1, element );
281 ++coarsened_;
282 }
283 }
284 }
285 }
286 }
287 }
288
293 {
294 int sumMarked = problem_.grid().comm().sum(marked_);
295 int sumCoarsened = problem_.grid().comm().sum(coarsened_);
296
297 return (sumMarked != 0 || sumCoarsened != 0);
298 }
299
306 void setLevels(int levMin, int levMax)
307 {
308 if (levMin < 0)
309 Dune::dgrave << __FILE__<< ":" <<__LINE__
310 << " : Dune cannot coarsen to gridlevels smaller 0! "<< std::endl;
311 levelMin_ = levMin;
312 levelMax_ = levMax;
313 }
314
322 int getMaxLevel() const
323 {
324 return levelMax_;
325 }
333 int getMinLevel() const
334 {
335 return levelMin_;
336 }
337
338 AdaptionIndicator& adaptionIndicator()
339 {
340 return adaptionIndicator_;
341 }
342
343 AdaptionIndicator& adaptionIndicator() const
344 {
345 return adaptionIndicator_;
346 }
347
348private:
361 bool checkNeighborsRefine_(const Element &entity, int level = 1)
362 {
363 // this also refines the neighbor elements
364 for(const auto& intersection : intersections(problem_.gridView(), entity))
365 {
366 if(!intersection.neighbor())
367 continue;
368
369 auto outside = intersection.outside();
370
371 // only mark non-ghost elements
372 if (outside.partitionType() == Dune::GhostEntity)
373 continue;
374
375 if ((outside.level() < levelMax_)
376 && (outside.level() < entity.level()))
377 {
378 problem_.grid().mark(1, outside);
379 ++marked_;
380
381 if(level != levelMax_)
382 checkNeighborsRefine_(outside, ++level);
383 }
384 }
385 return true;
386 }
387
388
398 void forceRefineRatio(int maxLevelDelta = 1)
399 {
400 LeafGridView leafGridView = problem_.gridView();
401 // delete all existing marks
402 problem_.grid().postAdapt();
403 bool done;
404 do
405 {
406 // run through all cells
407 done=true;
408 for (const auto& element : elements(problem_.gridView()))
409 {
410 // only mark non-ghost elements
411 if (element.partitionType() == Dune::GhostEntity)
412 continue;
413
414 // run through all neighbor-cells (intersections)
415 for (const auto& intersection : intersections(leafGridView, element))
416 {
417 if(!intersection.neighbor())
418 continue;
419
420 if (element.level() + maxLevelDelta < intersection.outside().level())
421 {
422 problem_.grid().mark( 1, element );
423 done=false;
424 }
425 }
426 }
427 if (done==false)
428 {
429 // adapt the grid
430 problem_.grid().adapt();
431 // delete marks
432 problem_.grid().postAdapt();
433 }
434 }
435 while (done!=true);
436 }
437
438 // private Variables
439 Problem& problem_;
440 AdaptionIndicator adaptionIndicator_;
441
442 int marked_;
443 int coarsened_;
444
445 int levelMin_;
446 int levelMax_;
447
448 int adaptationInterval_;
449};
450
458template<class TypeTag>
459class GridAdapt<TypeTag, false>
460{
464 using ScalarSolutionType = typename SolutionTypes::ScalarSolution;
465
466public:
467 void init()
468 {}
470 {}
472 {
473 return false;
474 }
475 void setLevels(int, int)
476 {}
477 void setTolerance(int, int)
478 {}
479 void setIndicator(const ScalarSolutionType&,
480 const Scalar&, const Scalar&)
481 {}
482 GridAdapt (Problem& problem)
483 {}
484};
485
486}
487#endif
Defines a type tag and some fundamental properties for linear solvers.
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type GetProp
get the type of a property (equivalent to old macro GET_PROP(...))
Definition: propertysystem.hh:140
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
Standard Module for h-adaptive simulations.
Definition: gridadapt.hh:47
void adaptGrid(Indicator &indicator)
Method to adapt the grid with individual indicator vector.
Definition: gridadapt.hh:151
GridAdapt(Problem &problem)
Definition: gridadapt.hh:68
int getMaxLevel() const
Returns maximum refinement level.
Definition: gridadapt.hh:322
AdaptionIndicator & adaptionIndicator()
Definition: gridadapt.hh:338
AdaptionIndicator & adaptionIndicator() const
Definition: gridadapt.hh:343
void adaptGrid()
Standard method to adapt the grid.
Definition: gridadapt.hh:129
void setLevels(int levMin, int levMax)
Definition: gridadapt.hh:306
bool wasAdapted()
Returns true if grid cells have been marked for adaptation.
Definition: gridadapt.hh:292
void markElements(Indicator &indicator)
Definition: gridadapt.hh:212
void init()
Initalization method of the h-adaptive module.
Definition: gridadapt.hh:87
int getMinLevel() const
Returns minimum refinement level.
Definition: gridadapt.hh:333
GridAdapt(Problem &problem)
Definition: gridadapt.hh:482
void adaptGrid()
Definition: gridadapt.hh:469
void setLevels(int, int)
Definition: gridadapt.hh:475
void setIndicator(const ScalarSolutionType &, const Scalar &, const Scalar &)
Definition: gridadapt.hh:479
void init()
Definition: gridadapt.hh:467
bool wasAdapted()
Definition: gridadapt.hh:471
void setTolerance(int, int)
Definition: gridadapt.hh:477
Base file for properties related to sequential models.