3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
gridadaptinitializationindicator.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 *****************************************************************************/
19#ifndef DUMUX_GRIDADAPTINITIALIZATIONINDICATOR_HH
20#define DUMUX_GRIDADAPTINITIALIZATIONINDICATOR_HH
21
22#include "properties.hh"
23
24#include <dune/common/dynvector.hh>
29namespace Dumux
30{
39template<class TypeTag>
40class GridAdaptInitializationIndicator
41{
42private:
43 using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
44 using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
45 using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
46 using Element = typename GridView::Traits::template Codim<0>::Entity;
47 using Intersection = typename GridView::Intersection;
48
49 using AdaptionIndicator = typename GET_PROP_TYPE(TypeTag, AdaptionIndicator);
50
51 using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
52 using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
53
54 enum
55 {
56 dim = GridView::dimension,
57 dimWorld = GridView::dimensionworld,
58 numEq = GET_PROP_VALUE(TypeTag, NumEq),
59 numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
60 };
61
62 enum
63 {
64 refineCell = 1,
65 coarsenCell = -1
66 };
67
68 using LocalPosition = Dune::FieldVector<Scalar, dim>;
69 using LocalPositionFace = Dune::FieldVector<Scalar, dim-1>;
70 using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
71
72 void virtualHierarchicSourceSearch_(PrimaryVariables &source, const Element& element)
73 {
74 int level = element.level();
75
76 if (level == maxAllowedLevel_)
77 {
78 GlobalPosition globalPos = element.geometry().center();
79 problem_.sourceAtPos(source, globalPos);
80
81 return;
82 }
83
84 unsigned int numRefine = maxAllowedLevel_ - level;
85 int numCheckCoords = pow(2, numRefine);
86
87 LocalPosition localPos(0.0);
88 GlobalPosition globalPosCheck(0.0);
89 Scalar halfInterval = (1.0/double(numCheckCoords))/2.;
90
91 PrimaryVariables sourceCheck(0.0);
92
93 using std::abs;
94 for (int i = 1; i <= numCheckCoords; i++)
95 {
96 for (int j = 1; j <= numCheckCoords; j++)
97 {
98 localPos[0] = double(i)/double(numCheckCoords) - halfInterval;
99 localPos[1] = double(j)/double(numCheckCoords) - halfInterval;
100 if (dim == 2)
101 {
102 globalPosCheck = element.geometry().global(localPos);
103 problem_.sourceAtPos(sourceCheck, globalPosCheck);
104
105 for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
106 {
107 if (abs(sourceCheck[eqIdx]) > abs(source[eqIdx]))
108 {
109 source[eqIdx] = sourceCheck[eqIdx];
110 }
111 }
112 }
113 else if (dim == 3)
114 {
115 for (int k = 1; k <= numCheckCoords; k++)
116 {
117 localPos[2] = double(k)/double(numCheckCoords) - halfInterval;
118 globalPosCheck = element.geometry().global(localPos);
119 problem_.sourceAtPos(sourceCheck, globalPosCheck);
120
121 for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
122 {
123 if (abs(sourceCheck[eqIdx]) > abs(source[eqIdx]))
124 {
125 source[eqIdx] = sourceCheck[eqIdx];
126 }
127 }
128 }
129 }
130 }
131 }
132 }
133
134 void virtualHierarchicBCSearch_(BoundaryTypes &bcTypes, PrimaryVariables &values, const Element& element, const Intersection& intersection)
135 {
136 int level = element.level();
137
138 if (level == maxAllowedLevel_)
139 {
140 GlobalPosition globalPos = intersection.geometry().center();
141 problem_.boundaryTypesAtPos(bcTypes, globalPos);
142
143 if (refineAtFluxBC_)
144 {
145 for (int i = 0; i < numEq; i++)
146 {
147 if (bcTypes.isNeumann(i))
148 {
149 PrimaryVariables fluxes;
150 problem_.neumannAtPos(fluxes, globalPos);
151
152 values += fluxes;
153 }
154 }
155 }
156 return;
157 }
158
159 unsigned int numRefine = maxAllowedLevel_ - level;
160 int numCheckCoords = pow(2, numRefine);
161
162 LocalPositionFace localPos(0.0);
163 GlobalPosition globalPosCheck(0.0);
164 Scalar halfInterval = (1.0/double(numCheckCoords))/2.;
165
166 PrimaryVariables fluxCheck(0.0);
167
168 for (int i = 1; i <= numCheckCoords; i++)
169 {
170 localPos[0] = double(i)/double(numCheckCoords) - halfInterval;
171 if (dim == 2)
172 {
173 globalPosCheck = intersection.geometry().global(localPos);
174 problem_.boundaryTypesAtPos(bcTypes, globalPosCheck);
175
176 for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
177 {
178 if (refineAtDirichletBC_ && bcTypes.isDirichlet(eqIdx))
179 {
180 return;
181 }
182 else if (refineAtFluxBC_ && bcTypes.isNeumann(eqIdx))
183 {
184 problem_.neumannAtPos(fluxCheck, globalPosCheck);
185
186 values += fluxCheck;
187 }
188 }
189 }
190 else if (dim == 3)
191 {
192 for (int k = 1; k <= numCheckCoords; k++)
193 {
194 localPos[1] = double(k)/double(numCheckCoords) - halfInterval;
195 globalPosCheck = intersection.geometry().global(localPos);
196
197 problem_.boundaryTypesAtPos(bcTypes, globalPosCheck);
198
199
200 for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
201 {
202 if (refineAtDirichletBC_ && bcTypes.isDirichlet(eqIdx))
203 {
204 return;
205 }
206 else if (refineAtFluxBC_ && bcTypes.isNeumann(eqIdx))
207 {
208 problem_.neumannAtPos(fluxCheck, globalPosCheck);
209
210 values += fluxCheck;
211 }
212
213 }
214 }
215 }
216 }
217 }
218
219
220public:
225 {
226 //First Adapt for boundary conditions and sources to get a correct pressure solution
227 if (nextMaxLevel_ == maxAllowedLevel_)
228 adaptionIndicator_.calculateIndicator();
229
230 // prepare an indicator for refinement
231 indicatorVector_.resize(problem_.variables().cellDataGlobal().size());
232
233 indicatorVector_ = coarsenCell;
234
235 if (!enableInitializationIndicator_)
236 return;
237
238 // 1) calculate Indicator -> min, maxvalues
239 // Schleife über alle Leaf-Elemente
240 using std::abs;
241 using std::max;
242 using std::min;
243 for (const auto& element : elements(problem_.gridView()))
244 {
245 int globalIdxI = problem_.variables().index(element);
246
247 int level = element.level();
248 maxLevel_ = max(level, maxLevel_);
249
250 if (level < minAllowedLevel_)
251 {
252 nextMaxLevel_ = min(max(level + 1, nextMaxLevel_), maxAllowedLevel_);
253 indicatorVector_[globalIdxI] = refineCell;
254 continue;
255 }
256
257 if (refineAtSource_)
258 {
259 PrimaryVariables source(0.0);
260 virtualHierarchicSourceSearch_(source, element);
261 for (int i = 0; i < numEq; i++)
262 {
263 if (abs(source[i]) > 1e-10)
264 {
265 nextMaxLevel_ = min(max(level + 1, nextMaxLevel_), maxAllowedLevel_);
266 indicatorVector_[globalIdxI] = refineCell;
267 break;
268 }
269 }
270 }
271
272 if (indicatorVector_[globalIdxI] != refineCell && (refineAtDirichletBC_ || refineAtFluxBC_))
273 {
274 // Berechne Verfeinerungsindikator an allen Zellen
275 for (const auto& intersection : intersections(problem_.gridView(), element))
276 {
277 if (intersection.boundary() && indicatorVector_[globalIdxI] != refineCell)
278 {
279 BoundaryTypes bcTypes;
280 PrimaryVariables values(0.0);
281
282 virtualHierarchicBCSearch_(bcTypes, values, element, intersection);
283
284
285 for (int i = 0; i < numEq; i++)
286 {
287 if (bcTypes.isDirichlet(i) && refineAtDirichletBC_)
288 {
289 nextMaxLevel_ = min(max(level + 1, nextMaxLevel_), maxAllowedLevel_);
290 indicatorVector_[globalIdxI] = refineCell;
291 break;
292 }
293 }
294 for (int j = 0; j < numPhases; j++)
295 {
296 if (abs(values[j]) > 1e-10)
297 {
298 nextMaxLevel_ = min(max(level + 1, nextMaxLevel_), maxAllowedLevel_);
299 indicatorVector_[globalIdxI] = refineCell;
300 break;
301 }
302
303 }
304 }
305 }
306 }
307 }
308 }
309
316 bool refine(const Element& element)
317 {
318 int idx = problem_.elementMapper().index(element);
319
320 if (indicatorVector_[idx] == refineCell)
321 return true;
322 else if (maxLevel_ == maxAllowedLevel_)
323 return adaptionIndicator_.refine(element);
324 else
325 return false;
326 }
327
334 bool coarsen(const Element& element)
335 {
336 int idx = problem_.elementMapper().index(element);
337
338 if (indicatorVector_[idx] == coarsenCell && maxLevel_ < maxAllowedLevel_)
339 return true;
340 else if (indicatorVector_[idx] == coarsenCell && !adaptionIndicator_.refine(element))
341 return true;
342 else
343 return false;
344 }
345
347 {
348 return maxLevel_;
349 }
350
352 void init()
353 {}
354
356 {
357 return nextMaxLevel_ == maxAllowedLevel_;
358 }
359
370 GridAdaptInitializationIndicator(Problem& problem, AdaptionIndicator& adaptionIndicator):
371 problem_(problem), adaptionIndicator_(adaptionIndicator), maxLevel_(0), nextMaxLevel_(0)
372 {
373 minAllowedLevel_ = getParam<int>("GridAdapt.MinLevel");
374 maxAllowedLevel_ = getParam<int>("GridAdapt.MaxLevel");
375 enableInitializationIndicator_ = getParam<bool>("GridAdapt.EnableInitializationIndicator");
376 refineAtDirichletBC_ = getParam<bool>("GridAdapt.RefineAtDirichletBC");
377 refineAtFluxBC_ = getParam<bool>("GridAdapt.RefineAtFluxBC");
378 refineAtSource_ = getParam<bool>("GridAdapt.RefineAtSource");
379
380 if (!refineAtDirichletBC_ && !refineAtFluxBC_ && !refineAtSource_)
381 {
382 nextMaxLevel_ = maxAllowedLevel_;
383 maxLevel_ = maxAllowedLevel_;
384 }
385 }
386
387private:
388 Problem& problem_;
389 AdaptionIndicator& adaptionIndicator_;
390 Dune::DynamicVector<int> indicatorVector_;
391 int maxLevel_;
392 int nextMaxLevel_;
393 int minAllowedLevel_;
394 int maxAllowedLevel_;
395 bool enableInitializationIndicator_;
396 bool refineAtDirichletBC_;
397 bool refineAtFluxBC_;
398 bool refineAtSource_;
399};
400}
401#endif
#define GET_PROP_VALUE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:282
#define GET_PROP_TYPE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:283
Dune::IteratorRange< typename MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper >::Intersections::const_iterator > intersections(const MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper > &glue)
Range generator to iterate with range-based for loops over all intersections as follows: for (const a...
Definition: glue.hh:62
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
Property tag NumEq
Number of equations in the system of PDEs.
Definition: porousmediumflow/sequential/properties.hh:68
Property tag NumPhases
Number of phases in the system.
Definition: porousmediumflow/sequential/properties.hh:69
void calculateIndicator()
Calculates the indicator used for refinement/coarsening for each grid cell.
Definition: gridadaptinitializationindicator.hh:224
bool initializeModel()
Definition: gridadaptinitializationindicator.hh:355
void init()
Initializes the adaption indicator class.
Definition: gridadaptinitializationindicator.hh:352
bool refine(const Element &element)
Indicator function for marking of grid cells for refinement.
Definition: gridadaptinitializationindicator.hh:316
int maxLevel()
Definition: gridadaptinitializationindicator.hh:346
GridAdaptInitializationIndicator(Problem &problem, AdaptionIndicator &adaptionIndicator)
Constructs a GridAdaptionIndicator instance.
Definition: gridadaptinitializationindicator.hh:370
bool coarsen(const Element &element)
Indicator function for marking of grid cells for coarsening.
Definition: gridadaptinitializationindicator.hh:334
Base file for properties related to sequential models.