adaptionhelper.hh 15 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
Thomas Fetzer's avatar
Thomas Fetzer committed
19
20
#ifndef DUMUX_ADAPTIONHELPER_HH
#define DUMUX_ADAPTIONHELPER_HH
21

Thomas Fetzer's avatar
Thomas Fetzer committed
22
23
#include <dune/common/version.hh>
#include <dune/grid/common/gridenums.hh>
24
#include <dune/grid/utility/persistentcontainer.hh>
25
#include <dune/localfunctions/lagrange/pqkfactory.hh>
Thomas Fetzer's avatar
Thomas Fetzer committed
26
27

#include <dumux/common/propertysystem.hh>
28
29
30
31
32
33
34
35
36

/**
 * @file
 * @brief  Base class holding the variables for implicit models.
 */

namespace Dumux
{

Thomas Fetzer's avatar
Thomas Fetzer committed
37
38
39
40
41
42
43
44
45
namespace Properties
{
NEW_PROP_TAG(GridView);
NEW_PROP_TAG(ImplicitIsBox);
NEW_PROP_TAG(PrimaryVariables);
NEW_PROP_TAG(Problem);
NEW_PROP_TAG(Scalar);
}

46
template<class TypeTag>
Thomas Fetzer's avatar
Thomas Fetzer committed
47
class AdaptionHelper
48
49
50
51
52
{
private:
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
53
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
54
55
56
57
58
59
60
61
62
63
64

    struct AdaptedValues
    {
    	PrimaryVariables u;
        int count;
        AdaptedValues()
        {
            count = 0;
        }
    };

65
66
67
68
69
70
71
72
73
    enum {
        // Grid and world dimension
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld
    };

    enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
    enum { dofCodim = isBox ? dim : 0 };

74
75
    typedef typename GridView::Grid Grid;
    typedef typename Grid::LevelGridView LevelGridView;
76
77
    typedef typename LevelGridView::template Codim<dofCodim>::Iterator LevelIterator;
    typedef typename LevelGridView::template Codim<0>::Iterator ElementLevelIterator;
78
79
    typedef typename GridView::Traits::template Codim<0>::EntityPointer ElementPointer;
    typedef typename GridView::Traits::template Codim<0>::Entity Element;
80
81
    typedef typename GridView::Traits::template Codim<dofCodim>::Entity DofEntity;
    typedef typename GridView::Traits::template Codim<dofCodim>::EntityPointer DofPointer;
82
83
    typedef Dune::PersistentContainer<Grid, AdaptedValues> PersistentContainer;

84
85
86
87
88
89
90
    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
    typedef typename GridView::ctype CoordScalar;
    typedef Dune::FieldVector<CoordScalar, dim> LocalPosition;

    typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> LocalFiniteElementCache;
    typedef typename LocalFiniteElementCache::FiniteElementType LocalFiniteElement;

91
92
93
private:
    const GridView gridView_;
    const Grid& grid_;
Thomas Fetzer's avatar
Thomas Fetzer committed
94
    PersistentContainer adaptionMap_;
95
96
97
98
99
100
101
102
103

public:
    //! Constructs an adaptive helper object
    /**
     * In addition to providing a storage object for cell-centered Methods, this class provides
     * mapping functionality to adapt the grid.
     *
     *  @param gridView a DUNE gridview object corresponding to diffusion and transport equation
     */
Thomas Fetzer's avatar
Thomas Fetzer committed
104
105
    AdaptionHelper(const GridView& gridView) :
    	gridView_(gridView), grid_(gridView.grid()), adaptionMap_(grid_, dofCodim)
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
    {}


    /*!
     * Store primary variables
     *
     * To reconstruct the solution in father elements, problem properties might
     * need to be accessed.
     * From upper level on downwards, the old solution is stored into an container
     * object, before the grid is adapted. Father elements hold averaged information
     * from the son cells for the case of the sons being coarsened.
     *
     * @param problem The current problem
     */
    void storePrimVars(Problem& problem)
    {
Thomas Fetzer's avatar
Thomas Fetzer committed
122
        adaptionMap_.resize();
123
124
125
126
127
128
129

        // loop over all levels of the grid
        for (int level = grid_.maxLevel(); level >= 0; level--)
        {
            //get grid view on level grid
            LevelGridView levelView = grid_.levelGridView(level);

130
131
132
133
134
            if(!isBox)
            {
				for (ElementLevelIterator eIt = levelView.template begin<0>(); eIt != levelView.template end<0>(); ++eIt)
				{
					//get your map entry
Thomas Fetzer's avatar
Thomas Fetzer committed
135
					AdaptedValues &adaptedValues = adaptionMap_[*eIt];
136
137
138
139
140
141
142

					// put your value in the map
					if (eIt->isLeaf())
					{
						// get index
						int indexI = this->elementIndex(problem, *eIt);

Thomas Fetzer's avatar
Thomas Fetzer committed
143
						storeAdaptionValues(adaptedValues, problem.model().curSol()[indexI]);
144
145
146
147
148
149
150

						adaptedValues.count = 1;
					}
					//Average in father
					if (eIt->level() > 0)
					{
						ElementPointer epFather = eIt->father();
Thomas Fetzer's avatar
Thomas Fetzer committed
151
						AdaptedValues& adaptedValuesFather = adaptionMap_[*epFather];
152
						adaptedValuesFather.count += 1;
Thomas Fetzer's avatar
Thomas Fetzer committed
153
						storeAdaptionValues(adaptedValues, adaptedValuesFather);
154
155
156
157
158
					}

				}
        	}
            else
159
            {
160
161
162
				for (LevelIterator dofIt = levelView.template begin<dofCodim>(); dofIt != levelView.template end<dofCodim>(); ++dofIt)
				{
					//get your map entry
Thomas Fetzer's avatar
Thomas Fetzer committed
163
					AdaptedValues &adaptedValues = adaptionMap_[*dofIt];
164

165
166
					// put your value in the map
					int indexI = this->dofIndex(problem, *dofIt);
167

Thomas Fetzer's avatar
Thomas Fetzer committed
168
					storeAdaptionValues(adaptedValues, problem.model().curSol()[indexI]);
169

170
171
172
					adaptedValues.count = 1;

				}
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
            }
        }
    }

    /*!
     * Reconstruct missing primary variables (where elements are created/deleted)
     *
     * To reconstruct the solution in father elements, problem properties might
     * need to be accessed.
     * Starting from the lowest level, the old solution is mapped on the new grid:
     * Where coarsened, new cells get information from old father element.
     * Where refined, a new solution is reconstructed from the old father cell,
     * and then a new son is created. That is then stored into the general data
     * structure (CellData).
     *
     * @param problem The current problem
     */
    void reconstructPrimVars(Problem& problem)
    {
Thomas Fetzer's avatar
Thomas Fetzer committed
192
        adaptionMap_.resize();
193
194
195
196
197

        for (int level = 0; level <= grid_.maxLevel(); level++)
        {
            LevelGridView levelView = grid_.levelGridView(level);

198
            for (ElementLevelIterator eIt = levelView.template begin<0>(); eIt != levelView.template end<0>(); ++eIt)
199
200
201
202
203
204
205
206
207
208
            {
                // only treat non-ghosts, ghost data is communicated afterwards
                if (eIt->partitionType() == Dune::GhostEntity)
                    continue;

                if (!eIt->isNew())
                {
                    //entry is in map, write in leaf
                    if (eIt->isLeaf())
                    {
209
210
                    	if(!isBox)
                    	{
Thomas Fetzer's avatar
Thomas Fetzer committed
211
							AdaptedValues &adaptedValues = adaptionMap_[*eIt];
212
213
							int newIdxI = this->elementIndex(problem, *eIt);

Thomas Fetzer's avatar
Thomas Fetzer committed
214
							setAdaptionValues(adaptedValues, problem.model().curSol()[newIdxI]);
215
216
217
218
219
220
221
222
223
224
225
226
						}
                    	else
                    	{
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
                    		unsigned int numSubEntities = eIt->subEntities(dofCodim);
#else
                    		int numSubEntities = eIt->template count <dofCodim>();
#endif

                        	for(unsigned int i = 0; i < numSubEntities; i++)
                        	{
                        		DofPointer subEntity = eIt->template subEntity <dofCodim>(i);
Thomas Fetzer's avatar
Thomas Fetzer committed
227
    							AdaptedValues &adaptedValues = adaptionMap_[*subEntity];
228
    							int newIdxI = this->dofIndex(problem, *subEntity);
229

Thomas Fetzer's avatar
Thomas Fetzer committed
230
    							setAdaptionValues(adaptedValues, problem.model().curSol()[newIdxI]);
231
232
233

                        	}
                    	}
234
235
236
237
238
                    }
                }
                else
                {
                    // value is not in map, interpolate from father element
239
                    if (eIt->level() > 0 && eIt->hasFather())
240
241
242
                    {
                        ElementPointer epFather = eIt->father();

243
                        if(!isBox)
244
                        {
Thomas Fetzer's avatar
Thomas Fetzer committed
245
246
247
							// create new entry: reconstruct from adaptionMap_[*father] to a new
							// adaptionMap_[*son]
							reconstructAdaptionValues(adaptionMap_, *epFather, *eIt, problem);
248
249

							// access new son
Thomas Fetzer's avatar
Thomas Fetzer committed
250
							AdaptedValues& adaptedValues = adaptionMap_[*eIt];
251
252
253
254
255
256
257
258
							adaptedValues.count = 1;

							// if we are on leaf, store reconstructed values of son in CellData object
							if (eIt->isLeaf())
							{
								// acess new CellData object
								int newIdxI = this->elementIndex(problem, *eIt);

Thomas Fetzer's avatar
Thomas Fetzer committed
259
								setAdaptionValues(adaptedValues, problem.model().curSol()[newIdxI]);
260
261
262
263
264
265
266
267
268
269
270
271
272
273
							}
                        }
                        else
                        {
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
                    		unsigned int numSubEntities= eIt->subEntities(dofCodim);
#else
                    		int numSubEntities= eIt->template count <dofCodim>();
#endif
                    		const auto geometryI = eIt->geometry();

                        	for(unsigned int i = 0; i < numSubEntities; i++)
                        	{
                        		DofPointer subEntity = eIt->template subEntity <dofCodim>(i);
Thomas Fetzer's avatar
Thomas Fetzer committed
274
    							AdaptedValues &adaptedValues = adaptionMap_[*subEntity];
275
276
277
278
279
280
281
282
283
284
285
286
287

    							if(adaptedValues.count == 0){
									LocalPosition dofCenterPos = geometryI.local(subEntity->geometry().center());
									const LocalFiniteElementCache feCache;
									Dune::GeometryType geomType = epFather->geometry().type();

									const LocalFiniteElement &localFiniteElement = feCache.get(geomType);
									std::vector<Dune::FieldVector<Scalar, 1> > shapeVal;
									localFiniteElement.localBasis().evaluateFunction(dofCenterPos, shapeVal);
									PrimaryVariables u(0);
									for (int j = 0; j < shapeVal.size(); ++j)
									{
										DofPointer subEntityFather = epFather->template subEntity <dofCodim>(j);
Thomas Fetzer's avatar
Thomas Fetzer committed
288
										AdaptedValues & adaptedValuesFather = adaptionMap_[*subEntityFather];
289
290
291
292
293
294
295
296
297
298
										u.axpy(shapeVal[j], adaptedValuesFather.u);
									}

									adaptedValues.u = u;
									adaptedValues.count = 1;
    							}

    							if (eIt->isLeaf())
    							{
        							int newIdxI = this->dofIndex(problem, *subEntity);
Thomas Fetzer's avatar
Thomas Fetzer committed
299
    								setAdaptionValues(adaptedValues, problem.model().curSol()[newIdxI]);
300
301
302
    							}

                        	}
303
304
                        }
                    }
305

306
307
308
309
310
                }
            }

        }
        // reset entries in restrictionmap
Thomas Fetzer's avatar
Thomas Fetzer committed
311
312
313
        adaptionMap_.resize( typename PersistentContainer::Value() );
        adaptionMap_.shrinkToFit();
        adaptionMap_.fill( typename PersistentContainer::Value() );
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329

//#if HAVE_MPI
//        // communicate ghost data
//        typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
//        typedef typename SolutionTypes::ElementMapper ElementMapper;
//        typedef VectorExchange<ElementMapper, std::vector<CellData> > DataHandle;
//        DataHandle dataHandle(problem.elementMapper(), this->cellDataGlobal());
//        problem.gridView().template communicate<DataHandle>(dataHandle,
//                                                            Dune::InteriorBorder_All_Interface,
//                                                            Dune::ForwardCommunication);
//#endif
    }

    //! Stores values to be adapted in an adaptedValues container
    /**
     * Stores values to be adapted from the current CellData objects into
Thomas Fetzer's avatar
Thomas Fetzer committed
330
     * the adaption container in order to be mapped on a new grid.
331
332
     *
     * \param adaptedValues Container for model-specific values to be adapted
333
     * \param u The variables to be stored
334
     */
Thomas Fetzer's avatar
Thomas Fetzer committed
335
    static void storeAdaptionValues(AdaptedValues& adaptedValues, const PrimaryVariables& u)
336
337
338
339
340
341
342
343
344
345
346
347
    {
        adaptedValues.u = u;
    }
    //! Stores sons entries into father element for averaging
    /**
     * Sum up the adaptedValues (sons values) into father element. We store from leaf
     * upwards, so sons are stored first, then cells on the next leaf (=fathers)
     * can be averaged.
     *
     * \param adaptedValues Container for model-specific values to be adapted
     * \param adaptedValuesFather Values to be adapted of father cell
     */
Thomas Fetzer's avatar
Thomas Fetzer committed
348
    static void storeAdaptionValues(AdaptedValues& adaptedValues,
Martin Schneider's avatar
Martin Schneider committed
349
                                    AdaptedValues& adaptedValuesFather)
350
    {
351
352
353
354
355
356
357
358
359
    	if(!isBox)
    	{
			adaptedValuesFather.u += adaptedValues.u;
			adaptedValuesFather.u /= adaptedValues.count;
    	}
    	else
    	{
    		adaptedValuesFather.u = adaptedValues.u;
    	}
360
361
362
363
364
365
366
367
    }
    //! Set adapted values in CellData
    /**
     * This methods stores reconstructed values into the cellData object, by
     * this setting a newly mapped solution to the storage container of the
     * decoupled models.
     *
     * \param adaptedValues Container for model-specific values to be adapted
368
     * \param u The variables to be stored
369
     */
Thomas Fetzer's avatar
Thomas Fetzer committed
370
    static void setAdaptionValues(AdaptedValues& adaptedValues, PrimaryVariables& u)
371
372
373
374
375
376
377
378
379
380
381
    {
    	PrimaryVariables uNew = adaptedValues.u;
    	uNew /= adaptedValues.count;

    	u = uNew;
    }

    //! Reconstructs sons entries from data of father cell
    /**
     * Reconstructs a new solution from a father cell into a newly
     * generated son cell. New cell is stored into the global
Thomas Fetzer's avatar
Thomas Fetzer committed
382
     * adaptionMap.
383
     *
Thomas Fetzer's avatar
Thomas Fetzer committed
384
     * \param adaptionMap Global map storing all values to be adapted
385
386
387
388
     * \param father Entity Pointer to the father cell
     * \param son Entity Pointer to the newly created son cell
     * \param problem The problem
     */
Thomas Fetzer's avatar
Thomas Fetzer committed
389
    static void reconstructAdaptionValues(Dune::PersistentContainer<Grid, AdaptedValues>& adaptionMap,
390
    		const Element& father, const Element& son, const Problem& problem)
391
    {
Thomas Fetzer's avatar
Thomas Fetzer committed
392
393
		AdaptedValues& adaptedValues = adaptionMap[son];
		AdaptedValues& adaptedValuesFather = adaptionMap[father];
394

395
396
397
		adaptedValues.u = adaptedValuesFather.u;
		adaptedValues.u /= adaptedValuesFather.count;
    }
398

399
400
401
402
403
404
405
    int dofIndex(const Problem& problem, const DofEntity& entity) const
    {
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
                return problem.model().dofMapper().index(entity);
#else
                return problem.model().dofMapper().map(entity);
#endif
406
407
    }

408
    int elementIndex(const Problem& problem, const Element& element) const
409
410
    {
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
411
                return problem.elementMapper().index(element);
412
#else
413
                return problem.elementMapper().map(element);
414
415
416
417
418
419
#endif
    }

};
}
#endif