fvgridgeometry.hh 33 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// -*- 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/>.   *
 *****************************************************************************/
/*!
 * \file
21
 * \brief Base class for the finite volume geometry vector for mpfa models
22
23
24
 *        This builds up the sub control volumes and sub control volume faces
 *        for each element.
 */
25
26
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_FV_GRID_GEOMETRY_BASE_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_FV_GRID_GEOMETRY_BASE_HH
27
28
29

#include <dune/geometry/multilineargeometry.hh>
#include <dune/geometry/referenceelements.hh>
30

31
#include <dumux/common/elementmap.hh>
Dennis Gläser's avatar
Dennis Gläser committed
32
#include <dumux/common/boundingboxtree.hh>
33

Dennis Gläser's avatar
Dennis Gläser committed
34
35
36
37
38
#include <dumux/discretization/basefvgridgeometry.hh>
#include <dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh>
#include <dumux/discretization/cellcentered/mpfa/connectivitymap.hh>
#include <dumux/discretization/cellcentered/mpfa/dualgridindexset.hh>
#include <dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh>
39
40
41
42
43
44
45
46
47
48

namespace Dumux
{
/*!
 * \ingroup ImplicitModel
 * \brief Base class for the finite volume geometry vector for mpfa models
 *        This builds up the sub control volumes and sub control volume faces
 *        for each element.
 */
template<class TypeTag, bool EnableFVElementGeometryCache>
49
class CCMpfaFVGridGeometry
50
{};
51

Dennis Gläser's avatar
Dennis Gläser committed
52
// specialization in case the finite volume grid geometries are stored
53
template<class TypeTag>
Dennis Gläser's avatar
Dennis Gläser committed
54
class CCMpfaFVGridGeometry<TypeTag, true> : public BaseFVGridGeometry<TypeTag>
55
{
Dennis Gläser's avatar
Dennis Gläser committed
56
57
    using ParentType = BaseFVGridGeometry<TypeTag>;
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
58
59
60
61
    using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
Dennis Gläser's avatar
Dennis Gläser committed
62
63
    using GridIVIndexSets = CCMpfaGridInteractionVolumeIndexSets<TypeTag>;
    using ConnectivityMap = CCMpfaConnectivityMap<TypeTag>;
64

Dennis Gläser's avatar
Dennis Gläser committed
65
66
67
68
69
70
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
    using Element = typename GridView::template Codim<0>::Entity;
    using Vertex = typename GridView::template Codim<dim>::Entity;
    using Intersection = typename GridView::Intersection;
71
72
73
    using CoordScalar = typename GridView::ctype;
    using IndexType = typename GridView::IndexSet::IndexType;

Dennis Gläser's avatar
Dennis Gläser committed
74
    using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
75
76
    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;

Dennis Gläser's avatar
Dennis Gläser committed
77
78
79
    //! The local class needs access to the scv, scvfs and the fv element geometry
    //! as they are globally cached
    friend typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
80
81
82

public:
    //! Constructor
83
    CCMpfaFVGridGeometry(const GridView gridView)
Dennis Gläser's avatar
Dennis Gläser committed
84
85
    : ParentType(gridView), elementMap_(gridView)
    {}
86

Dennis Gläser's avatar
Dennis Gläser committed
87
88
89
    /*!
     * \brief Returns the total number of sub control volumes.
     */
90
91
92
    std::size_t numScv() const
    { return scvs_.size(); }

Dennis Gläser's avatar
Dennis Gläser committed
93
94
95
    /*!
     * \brief Returns the total number of sub control volume faces.
     */
96
97
98
    std::size_t numScvf() const
    { return scvfs_.size(); }

Dennis Gläser's avatar
Dennis Gläser committed
99
100
101
102
103
    /*!
     * \brief Returns the number of scvfs on the domain boundary.
     */
    std::size_t numBoundaryScvf() const
    { return numBoundaryScvf_; }
104

105
106
107
108
109
110
111

    /*!
     * \brief Returns the total number of degrees of freedom.
     */
    std::size_t numDofs() const
    { return this->gridView().size(0); }

Dennis Gläser's avatar
Dennis Gläser committed
112
113
114
    /*!
     * \brief Gets an element from a sub control volume contained in it.
     */
115
116
117
    Element element(const SubControlVolume& scv) const
    { return elementMap_.element(scv.elementIndex()); }

Dennis Gläser's avatar
Dennis Gläser committed
118
119
120
    /*!
     * \brief Gets an element from a global element index.
     */
121
122
123
    Element element(IndexType eIdx) const
    { return elementMap_.element(eIdx); }

Dennis Gläser's avatar
Dennis Gläser committed
124
125
126
127
128
129
    /*!
     * \brief Returns true if primary interaction volumes are used around a given vertex,
     *        false otherwise.
     */
    bool vertexUsesPrimaryInteractionVolume(const Vertex& v) const
    { return primaryInteractionVolumeVertices_[this->vertexMapper().index(v)]; }
130

Dennis Gläser's avatar
Dennis Gläser committed
131
132
133
134
135
136
    /*!
     * \brief Returns true if primary interaction volumes are used around a given vertex index,
     *        false otherwise.
     */
    bool vertexUsesPrimaryInteractionVolume(IndexType vIdxGlobal) const
    { return primaryInteractionVolumeVertices_[vIdxGlobal]; }
137

Dennis Gläser's avatar
Dennis Gläser committed
138
139
140
141
142
    /*!
     * \brief Returns if primary interaction volumes are used around a given vertex.
     */
    bool vertexUsesSecondaryInteractionVolume(const Vertex& v) const
    { return secondaryInteractionVolumeVertices_[this->vertexMapper().index(v)]; }
143

Dennis Gläser's avatar
Dennis Gläser committed
144
145
146
147
148
149
    /*!
     * \brief Returns true if primary interaction volumes are used around a given vertex index,
     *        false otherwise.
     */
    bool vertexUsesSecondaryInteractionVolume(IndexType vIdxGlobal) const
    { return secondaryInteractionVolumeVertices_[vIdxGlobal]; }
150

Dennis Gläser's avatar
Dennis Gläser committed
151
152
153
154
155
156
157
158
159
160
161
    /*!
     * \brief Updates all finite volume geometries of the grid. This has to be called again
     *        after grid adaption. A function can be passed to this method specifying where the secondary
     *        interaction volume type should be used. Per default we use it on boundaries
     *        and branching points.
     *
     * \param useSecondaryIV Indicator function at which vertices to apply the secondary interaction volume
     */
    void update( std::function<bool(const Element&, const Intersection&, bool)> useSecondaryIV
                              = [] (const Element& e, const Intersection& is, bool isBranching)
                                   { return is.boundary() || isBranching; } )
162
    {
Dennis Gläser's avatar
Dennis Gläser committed
163
164
        // stop the time required for the update
        Dune::Timer timer;
165

Dennis Gläser's avatar
Dennis Gläser committed
166
167
168
169
170
171
172
173
        // clear containers (necessary after grid refinement)
        scvs_.clear();
        scvfs_.clear();
        scvfIndicesOfScv_.clear();
        elementMap_.clear();

        // determine the number of geometric entities
        const auto numVert = this->gridView().size(dim);
174
        const auto numScvs = numDofs();
Dennis Gläser's avatar
Dennis Gläser committed
175
        std::size_t numScvf = MpfaHelper::estimateNumScvf(this->gridView());
176

Dennis Gläser's avatar
Dennis Gläser committed
177
        // resize containers
178
        scvs_.resize(numScvs);
Dennis Gläser's avatar
Dennis Gläser committed
179
        scvfs_.reserve(numScvf);
180
181
182
        scvfIndicesOfScv_.resize(numScvs);
        elementMap_.resize(numScvs);

Dennis Gläser's avatar
Dennis Gläser committed
183
184
185
186
        // Some methods require to use a second type of interaction volume, e.g.
        // around vertices on the boundary or branching points (surface grids)
        primaryInteractionVolumeVertices_.resize(numVert, true);
        secondaryInteractionVolumeVertices_.resize(numVert, false);
187
188

        // find vertices on processor boundaries
Dennis Gläser's avatar
Dennis Gläser committed
189
        const auto isGhostVertex = MpfaHelper::findGhostVertices(this->gridView(), this->vertexMapper());
190

Dennis Gläser's avatar
Dennis Gläser committed
191
192
        // instantiate the dual grid index set (to be used for construction of interaction volumes)
        CCMpfaDualGridIndexSet<TypeTag> dualIdSet(this->gridView());
193
194
195

        // Build the SCVs and SCV faces
        IndexType scvfIdx = 0;
Dennis Gläser's avatar
Dennis Gläser committed
196
197
        numBoundaryScvf_ = 0;
        for (const auto& element : elements(this->gridView()))
198
        {
Dennis Gläser's avatar
Dennis Gläser committed
199
            const auto eIdx = this->elementMapper().index(element);
200
201
202
203
204
205
206
207
208
209
210

            // the element geometry
            auto elementGeometry = element.geometry();

            // fill the element map with seeds
            elementMap_[eIdx] = element.seed();

            // The local scvf index set
            std::vector<IndexType> scvfIndexSet;
            scvfIndexSet.reserve(MpfaHelper::getNumLocalScvfs(elementGeometry.type()));

Dennis Gläser's avatar
Dennis Gläser committed
211
            // for network grids there might be multiple intersections with the same geometryInInside
212
213
214
215
216
            // we indentify those by the indexInInside for now (assumes conforming grids at branching facets)
            std::vector<std::vector<IndexType>> outsideIndices;
            if (dim < dimWorld)
            {
                outsideIndices.resize(element.subEntities(1));
Dennis Gläser's avatar
Dennis Gläser committed
217
                for (const auto& intersection : intersections(this->gridView(), element))
218
219
220
221
                {
                    if (intersection.neighbor())
                    {
                        const auto& outside = intersection.outside();
Dennis Gläser's avatar
Dennis Gläser committed
222
                        auto nIdx = this->elementMapper().index(outside);
223
224
225
226
227
228
                        outsideIndices[intersection.indexInInside()].push_back(nIdx);
                    }
                }
            }

            // construct the sub control volume faces
Dennis Gläser's avatar
Dennis Gläser committed
229
            for (const auto& is : intersections(this->gridView(), element))
230
231
232
233
234
            {
                const auto indexInInside = is.indexInInside();
                const bool boundary = is.boundary();
                const bool neighbor = is.neighbor();

Dennis Gläser's avatar
Dennis Gläser committed
235
                // for surface grids, skip the rest if handled already
236
237
238
239
240
241
242
243
244
245
                if (dim < dimWorld && neighbor && outsideIndices[indexInInside].empty())
                    continue;

                // if outside level > inside level, use the outside element in the following
                const bool useNeighbor = neighbor && is.outside().level() > element.level();
                const auto& e = useNeighbor ? is.outside() : element;
                const auto indexInElement = useNeighbor ? is.indexInOutside() : indexInInside;
                const auto eg = e.geometry();
                const auto& refElement = ReferenceElements::general(eg.type());

Dennis Gläser's avatar
Dennis Gläser committed
246
247
248
249
250
251
252
253
254
255
                // Set up a container with all relevant positions for scvf corner computation
                const auto numCorners = is.geometry().corners();
                const auto isPositions = MpfaHelper::computeScvfCornersOnIntersection(eg,
                                                                                      refElement,
                                                                                      indexInElement,
                                                                                      numCorners);

                // evaluate if vertices on this intersection use primary/secondary IVs
                const bool isBranchingPoint = dim < dimWorld ? outsideIndices[indexInInside].size() > 1 : false;
                const bool usesSecondaryIV = useSecondaryIV(element, is, isBranchingPoint);
256
257

                // make the scv faces belonging to each corner of the intersection
Dennis Gläser's avatar
Dennis Gläser committed
258
                for (int c = 0; c < numCorners; ++c)
259
260
261
                {
                    // get the global vertex index the scv face is connected to
                    const auto vIdxLocal = refElement.subEntity(indexInElement, 1, c, dim);
Dennis Gläser's avatar
Dennis Gläser committed
262
                    const auto vIdxGlobal = this->vertexMapper().subIndex(e, vIdxLocal, dim);
263
264
265
266
267

                    // do not build scvfs connected to a processor boundary
                    if (isGhostVertex[vIdxGlobal])
                        continue;

Dennis Gläser's avatar
Dennis Gläser committed
268
269
                    // if this vertex is tagged to use the secondary IVs, store info
                    if (usesSecondaryIV)
270
                    {
Dennis Gläser's avatar
Dennis Gläser committed
271
272
                        primaryInteractionVolumeVertices_[vIdxGlobal] = false;
                        secondaryInteractionVolumeVertices_[vIdxGlobal] = true;
273
274
                    }

Dennis Gläser's avatar
Dennis Gläser committed
275
276
                    // the quadrature point to be used on the scvfs
                    static const Scalar q = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Mpfa, Q);
277

Dennis Gläser's avatar
Dennis Gläser committed
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
                    // make the scv face (for non-boundary scvfs on network grids, use precalculated outside indices)
                    const auto& outsideScvIndices = [&] ()
                                                    {
                                                        if (!boundary)
                                                        {
                                                            return dim == dimWorld ?
                                                                   std::vector<IndexType>({this->elementMapper().index(is.outside())}) :
                                                                   outsideIndices[indexInInside];
                                                        }
                                                        else
                                                        {
                                                            // compute boundary scv idx and increment counter
                                                            const IndexType bIdx = numScvs + numBoundaryScvf_;
                                                            numBoundaryScvf_++;
                                                            return std::vector<IndexType>(1, bIdx);
                                                        }
                                                    } ();

                    scvfIndexSet.push_back(scvfIdx);
                    scvfs_.emplace_back(MpfaHelper(),
                                        MpfaHelper::getScvfCorners(isPositions, numCorners, c),
                                        is.centerUnitOuterNormal(),
                                        vIdxGlobal,
                                        vIdxLocal,
                                        scvfIdx,
                                        eIdx,
                                        outsideScvIndices,
                                        q,
                                        boundary);

                    // insert the scvf data into the dual grid index set
                    dualIdSet[vIdxGlobal].insert(scvfs_.back());
310
311
312
313
314
315

                    // increment scvf counter
                    scvfIdx++;
                }

                // for network grids, clear outside indices to not make a second scvf on that facet
Dennis Gläser's avatar
Dennis Gläser committed
316
                if (dim < dimWorld) outsideIndices[indexInInside].clear();
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
            }

            // create sub control volume for this element
            scvs_[eIdx] = SubControlVolume(std::move(elementGeometry), eIdx);

            // Save the scvf indices belonging to this scv to build up fv element geometries fast
            scvfIndicesOfScv_[eIdx] = scvfIndexSet;
        }

        // Make the flip index set for network and surface grids
        if (dim < dimWorld)
        {
            flipScvfIndices_.resize(scvfs_.size());
            for (auto&& scvf : scvfs_)
            {
                if (scvf.boundary())
                    continue;

                const auto numOutsideScvs = scvf.numOutsideScvs();
                const auto vIdxGlobal = scvf.vertexIndex();
                const auto insideScvIdx = scvf.insideScvIdx();

                flipScvfIndices_[scvf.index()].resize(numOutsideScvs);
                for (unsigned int i = 0; i < numOutsideScvs; ++i)
                {
                    const auto outsideScvIdx = scvf.outsideScvIdx(i);
                    for (auto outsideScvfIndex : scvfIndicesOfScv_[outsideScvIdx])
                    {
                        const auto& outsideScvf = this->scvf(outsideScvfIndex);
                        if (outsideScvf.vertexIndex() == vIdxGlobal &&
Dennis Gläser's avatar
Dennis Gläser committed
347
                            MpfaHelper::vectorContainsValue(outsideScvf.outsideScvIndices(), insideScvIdx))
348
349
350
351
352
353
354
355
356
357
358
                        {
                            flipScvfIndices_[scvf.index()][i] = outsideScvfIndex;
                            // there is always only one flip face in an outside element
                            break;
                        }
                    }

                }
            }
        }

359
        // building the geometries has finished
Dennis Gläser's avatar
Dennis Gläser committed
360
361
362
363
364
365
366
367
368
369
370
        std::cout << "Initializing of the grid finite volume geometry took " << timer.elapsed() << " seconds." << std::endl;

        // Initialize the grid interaction volume seeds
        timer.reset();
        ivIndexSets_.update(*this, std::move(dualIdSet));
        std::cout << "Initializing of the grid interaction volume index sets took " << timer.elapsed() << " seconds." << std::endl;

        // build the connectivity map for an effecient assembly
        timer.reset();
        connectivityMap_.update(*this);
        std::cout << "Initializing of the connectivity map took " << timer.elapsed() << " seconds." << std::endl;
371
372
373
    }

    /*!
Dennis Gläser's avatar
Dennis Gläser committed
374
     * \brief Returns a sub control volume for a given global scv index.
375
376
377
378
     */
    const SubControlVolume& scv(IndexType scvIdx) const
    { return scvs_[scvIdx]; }

Dennis Gläser's avatar
Dennis Gläser committed
379
380
381
    /*!
     * \brief Returns a sub control volume face for a given global scvf index.
     */
382
383
384
    const SubControlVolumeFace& scvf(IndexType scvfIdx) const
    { return scvfs_[scvfIdx]; }

Dennis Gläser's avatar
Dennis Gläser committed
385
386
387
388
    /*!
     * \brief Returns the "flipped" scvf, i.e. the correspongin scvf in an outside element.
     *        The second argument is optional and only comes into play on network/surface grids.
     */
389
390
391
    const SubControlVolumeFace& flipScvf(IndexType scvfIdx, unsigned int outsideScvfIdx = 0) const
    { return scvfs_[flipScvfIndices_[scvfIdx][outsideScvfIdx]]; }

Dennis Gläser's avatar
Dennis Gläser committed
392
393
394
    /*!
     * \brief Returns the sub control volume face indices of an scv by global index.
     */
395
396
397
    const std::vector<IndexType>& scvfIndicesOfScv(IndexType scvIdx) const
    { return scvfIndicesOfScv_[scvIdx]; }

Dennis Gläser's avatar
Dennis Gläser committed
398
399
400
401
402
403
    /*!
     * \brief Returns the connectivity map of which dofs have derivatives with respect
     *        to a given dof.
     */
    const ConnectivityMap& connectivityMap() const
    { return connectivityMap_; }
404

Dennis Gläser's avatar
Dennis Gläser committed
405
406
407
408
409
    /*!
     * \brief Returns the grit interaction volume seeds class.
     */
    const GridIVIndexSets& gridInteractionVolumeIndexSets() const
    { return ivIndexSets_; }
410

Dennis Gläser's avatar
Dennis Gläser committed
411
412
private:
    // mappers
413
    Dumux::ElementMap<GridView> elementMap_;
Dennis Gläser's avatar
Dennis Gläser committed
414
415
416
    ConnectivityMap connectivityMap_;

    // the finite volume grid geometries
417
418
419
    std::vector<SubControlVolume> scvs_;
    std::vector<SubControlVolumeFace> scvfs_;

Dennis Gläser's avatar
Dennis Gläser committed
420
    // containers storing the global data
421
    std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
Dennis Gläser's avatar
Dennis Gläser committed
422
423
424
425
    std::vector<bool> primaryInteractionVolumeVertices_;
    std::vector<bool> secondaryInteractionVolumeVertices_;
    std::size_t numBoundaryScvf_;

426
427
428
    // needed for embedded surface and network grids (dim < dimWorld)
    std::vector<std::vector<IndexType>> flipScvfIndices_;

Dennis Gläser's avatar
Dennis Gläser committed
429
430
    // The grid interaction volume index set
    GridIVIndexSets ivIndexSets_;
431
432
433
434
};

// specialization in case the FVElementGeometries are not stored
template<class TypeTag>
Dennis Gläser's avatar
Dennis Gläser committed
435
class CCMpfaFVGridGeometry<TypeTag, false> : public BaseFVGridGeometry<TypeTag>
436
{
Dennis Gläser's avatar
Dennis Gläser committed
437
438
    using ParentType = BaseFVGridGeometry<TypeTag>;
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
439
440
441
442
    using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
Dennis Gläser's avatar
Dennis Gläser committed
443
444
    using GridIVIndexSets = CCMpfaGridInteractionVolumeIndexSets<TypeTag>;
    using ConnectivityMap = CCMpfaConnectivityMap<TypeTag>;
445

Dennis Gläser's avatar
Dennis Gläser committed
446
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
447
448
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
Dennis Gläser's avatar
Dennis Gläser committed
449
450
451
452
453
    using Element = typename GridView::template Codim<0>::Entity;
    using Vertex = typename GridView::template Codim<dim>::Entity;
    using Intersection = typename GridView::Intersection;
    using CoordScalar = typename GridView::ctype;
    using IndexType = typename GridView::IndexSet::IndexType;
454

Dennis Gläser's avatar
Dennis Gläser committed
455
    using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
456
457
    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;

Dennis Gläser's avatar
Dennis Gläser committed
458
459
    //! Todo is this necessary?
    friend typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
460
461
462

public:
    //! Constructor
463
    CCMpfaFVGridGeometry(const GridView gridView)
Dennis Gläser's avatar
Dennis Gläser committed
464
465
    : ParentType(gridView), elementMap_(gridView)
    {}
466

Dennis Gläser's avatar
Dennis Gläser committed
467
468
469
    /*!
     * \brief Returns the total number of sub control volumes.
     */
470
471
472
    std::size_t numScv() const
    { return numScvs_; }

Dennis Gläser's avatar
Dennis Gläser committed
473
474
475
    /*!
     * \brief Returns the total number of sub control volume faces.
     */
476
477
478
    std::size_t numScvf() const
    { return numScvf_; }

Dennis Gläser's avatar
Dennis Gläser committed
479
480
481
482
483
    /*!
     * \brief Returns the number of scvfs on the domain boundary.
     */
    std::size_t numBoundaryScvf() const
    { return numBoundaryScvf_; }
484

485
486
487
488
489
490
    /*!
     * \brief Returns the total number of degrees of freedom.
     */
    std::size_t numDofs() const
    { return this->gridView().size(0); }

Dennis Gläser's avatar
Dennis Gläser committed
491
492
493
    /*!
     * \brief Gets an element from a sub control volume contained in it.
     */
494
495
496
    Element element(const SubControlVolume& scv) const
    { return elementMap_.element(scv.elementIndex()); }

Dennis Gläser's avatar
Dennis Gläser committed
497
498
499
    /*!
     * \brief Gets an element from a global element index.
     */
500
501
502
    Element element(IndexType eIdx) const
    { return elementMap_.element(eIdx); }

Dennis Gläser's avatar
Dennis Gläser committed
503
504
505
506
507
508
    /*!
     * \brief Returns true if primary interaction volumes are used around a given vertex,
     *        false otherwise.
     */
    bool vertexUsesPrimaryInteractionVolume(const Vertex& v) const
    { return primaryInteractionVolumeVertices_[this->vertexMapper().index(v)]; }
509

Dennis Gläser's avatar
Dennis Gläser committed
510
511
512
513
514
515
    /*!
     * \brief Returns true if primary interaction volumes are used around a given vertex index,
     *        false otherwise.
     */
    bool vertexUsesPrimaryInteractionVolume(IndexType vIdxGlobal) const
    { return primaryInteractionVolumeVertices_[vIdxGlobal]; }
516

Dennis Gläser's avatar
Dennis Gläser committed
517
518
519
520
521
    /*!
     * \brief Returns if primary interaction volumes are used around a given vertex.
     */
    bool vertexUsesSecondaryInteractionVolume(const Vertex& v) const
    { return secondaryInteractionVolumeVertices_[this->vertexMapper().index(v)]; }
522

Dennis Gläser's avatar
Dennis Gläser committed
523
524
525
526
527
528
    /*!
     * \brief Returns true if primary interaction volumes are used around a given vertex index,
     *        false otherwise.
     */
    bool vertexUsesSecondaryInteractionVolume(IndexType vIdxGlobal) const
    { return secondaryInteractionVolumeVertices_[vIdxGlobal]; }
529

Dennis Gläser's avatar
Dennis Gläser committed
530
531
532
533
534
    /*!
     * \brief Returns true if a given vertex lies on a processor boundary inside a ghost element.
     */
    bool isGhostVertex(const Vertex& v) const
    { return isGhostVertex_[this->vertexMapper().index(v)]; }
535

Dennis Gläser's avatar
Dennis Gläser committed
536
537
538
539
540
541
    /*!
     * \brief Returns true if the vertex corresponding to a given vertex index lies on a
     *        processor boundary inside a ghost element.
     */
    bool isGhostVertex(IndexType vIdxGlobal) const
    { return isGhostVertex_[vIdxGlobal]; }
542

Dennis Gläser's avatar
Dennis Gläser committed
543
544
545
546
547
548
549
550
551
552
553
    /*!
     * \brief Updates all finite volume geometries of the grid. This has to be called again
     *        after grid adaption. A function can be passed to this method specifying where the secondary
     *        interaction volume type should be used. Per default we use it on boundaries
     *        and branching points.
     *
     * \param useSecondaryIV Indicator function at which vertices to apply the secondary interaction volume
     */
    void update( std::function<bool(const Element&, const Intersection&, bool)> useSecondaryIV
                              = [] (const Element& e, const Intersection& is, bool isBranching)
                                   { return is.boundary() || isBranching; } )
554
    {
Dennis Gläser's avatar
Dennis Gläser committed
555
556
        // stop the time required for the update
        Dune::Timer timer;
557

Dennis Gläser's avatar
Dennis Gläser committed
558
559
560
561
        // clear containers (necessary after grid refinement)
        scvfIndicesOfScv_.clear();
        neighborVolVarIndices_.clear();
        elementMap_.clear();
562

Dennis Gläser's avatar
Dennis Gläser committed
563
        // resize containers
564
        numScvs_ = numDofs();
565
566
        scvfIndicesOfScv_.resize(numScvs_);
        neighborVolVarIndices_.resize(numScvs_);
Dennis Gläser's avatar
Dennis Gläser committed
567
        elementMap_.resize(numScvs_);
568

Dennis Gläser's avatar
Dennis Gläser committed
569
570
571
572
573
        // Some methods require to use a second type of interaction volume, e.g.
        // around vertices on the boundary or branching points (surface grids)
        const auto numVert = this->gridView().size(dim);
        primaryInteractionVolumeVertices_.resize(numVert, true);
        secondaryInteractionVolumeVertices_.resize(numVert, false);
574

Dennis Gläser's avatar
Dennis Gläser committed
575
576
        // find vertices on processor boundaries HERE!!
        isGhostVertex_ = MpfaHelper::findGhostVertices(this->gridView(), this->vertexMapper());
577

Dennis Gläser's avatar
Dennis Gläser committed
578
579
        // instantiate the dual grid index set (to be used for construction of interaction volumes)
        CCMpfaDualGridIndexSet<TypeTag> dualIdSet(this->gridView());
580

Dennis Gläser's avatar
Dennis Gläser committed
581
        // Build the SCVs and SCV faces
582
        numScvf_ = 0;
Dennis Gläser's avatar
Dennis Gläser committed
583
584
        numBoundaryScvf_ = 0;
        for (const auto& element : elements(this->gridView()))
585
        {
Dennis Gläser's avatar
Dennis Gläser committed
586
587
588
589
            const auto eIdx = this->elementMapper().index(element);

            // the element geometry
            auto elementGeometry = element.geometry();
590
591
592
593
594

            // fill the element map with seeds
            elementMap_[eIdx] = element.seed();

            // the element-wise index sets for finite volume geometry
Dennis Gläser's avatar
Dennis Gläser committed
595
            const auto numLocalFaces = MpfaHelper::getNumLocalScvfs(elementGeometry.type());
596
597
598
599
600
            std::vector<IndexType> scvfsIndexSet;
            std::vector< std::vector<IndexType> > neighborVolVarIndexSet;
            scvfsIndexSet.reserve(numLocalFaces);
            neighborVolVarIndexSet.reserve(numLocalFaces);

Dennis Gläser's avatar
Dennis Gläser committed
601
            // for network grids there might be multiple intersections with the same geometryInInside
602
603
604
605
606
            // we indentify those by the indexInInside for now (assumes conforming grids at branching facets)
            std::vector<std::vector<IndexType>> outsideIndices;
            if (dim < dimWorld)
            {
                outsideIndices.resize(element.subEntities(1));
Dennis Gläser's avatar
Dennis Gläser committed
607
                for (const auto& intersection : intersections(this->gridView(), element))
608
609
610
611
                {
                    if (intersection.neighbor())
                    {
                        const auto& outside = intersection.outside();
Dennis Gläser's avatar
Dennis Gläser committed
612
                        auto nIdx = this->elementMapper().index(outside);
613
614
615
616
617
618
619
                        outsideIndices[intersection.indexInInside()].push_back(nIdx);
                    }
                }
            }

            unsigned int localFaceIdx = 0;
            // construct the sub control volume faces
Dennis Gläser's avatar
Dennis Gläser committed
620
            for (const auto& is : intersections(this->gridView(), element))
621
622
623
624
625
            {
                const auto indexInInside = is.indexInInside();
                const bool boundary = is.boundary();
                const bool neighbor = is.neighbor();

Dennis Gläser's avatar
Dennis Gläser committed
626
                // for surface grids, skip the rest if handled already
627
628
629
630
631
632
633
634
635
636
                if (dim < dimWorld && neighbor && outsideIndices[indexInInside].empty())
                    continue;

                // if outside level > inside level, use the outside element in the following
                const bool useNeighbor = neighbor && is.outside().level() > element.level();
                const auto& e = useNeighbor ? is.outside() : element;
                const auto indexInElement = useNeighbor ? is.indexInOutside() : indexInInside;
                const auto eg = e.geometry();
                const auto& refElement = ReferenceElements::general(eg.type());

Dennis Gläser's avatar
Dennis Gläser committed
637
638
639
640
641
642
                // evaluate if vertices on this intersection use primary/secondary IVs
                const bool isBranchingPoint = dim < dimWorld ? outsideIndices[indexInInside].size() > 1 : false;
                const bool usesSecondaryIV = useSecondaryIV(element, is, isBranchingPoint);

                // make the scv faces belonging to each corner of the intersection
                for (int c = 0; c < is.geometry().corners(); ++c)
643
                {
Dennis Gläser's avatar
Dennis Gläser committed
644
                    // get the global vertex index the scv face is connected to
645
                    const auto vIdxLocal = refElement.subEntity(indexInElement, 1, c, dim);
Dennis Gläser's avatar
Dennis Gläser committed
646
                    const auto vIdxGlobal = this->vertexMapper().subIndex(e, vIdxLocal, dim);
647
648

                    // do not build scvfs connected to a processor boundary
Dennis Gläser's avatar
Dennis Gläser committed
649
                    if (isGhostVertex_[vIdxGlobal])
650
651
                        continue;

Dennis Gläser's avatar
Dennis Gläser committed
652
653
                    // if this vertex is tagged to use the secondary IVs, store info
                    if (usesSecondaryIV)
654
                    {
Dennis Gläser's avatar
Dennis Gläser committed
655
656
                        primaryInteractionVolumeVertices_[vIdxGlobal] = false;
                        secondaryInteractionVolumeVertices_[vIdxGlobal] = true;
657
658
                    }

Dennis Gläser's avatar
Dennis Gläser committed
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
                    // make the scv face (for non-boundary scvfs on network grids, use precalculated outside indices)
                    const auto& outsideScvIndices = [&] ()
                                                    {
                                                        if (!boundary)
                                                        {
                                                            return dim == dimWorld ?
                                                                   std::vector<IndexType>({this->elementMapper().index(is.outside())}) :
                                                                   outsideIndices[indexInInside];
                                                        }
                                                        else
                                                        {
                                                            // compute boundary scv idx and increment counter
                                                            const IndexType bIdx = numScvs_ + numBoundaryScvf_;
                                                            numBoundaryScvf_++;
                                                            return std::vector<IndexType>(1, bIdx);
                                                        }
                                                    } ();

                    // insert the scvf data into the dual grid index set
                    dualIdSet[vIdxGlobal].insert(boundary, numScvf_, eIdx, outsideScvIndices);

                    // store information on the scv face
                    scvfsIndexSet.push_back(numScvf_++);
                    neighborVolVarIndexSet.emplace_back(std::move(outsideScvIndices));
683
684
685
686
687
688

                    // increment counter
                    localFaceIdx++;
                }

                // for network grids, clear outside indices to not make a second scvf on that facet
Dennis Gläser's avatar
Dennis Gläser committed
689
                if (dim < dimWorld) outsideIndices[indexInInside].clear();
690
691
692
693
694
695
696
            }

            // store the sets of indices in the data container
            scvfIndicesOfScv_[eIdx] = scvfsIndexSet;
            neighborVolVarIndices_[eIdx] = neighborVolVarIndexSet;
        }

697
        // building the geometries has finished
Dennis Gläser's avatar
Dennis Gläser committed
698
699
700
701
702
703
704
705
706
707
708
        std::cout << "Initializing of the grid finite volume geometry took " << timer.elapsed() << " seconds." << std::endl;

        // Initialize the grid interaction volume seeds
        timer.reset();
        ivIndexSets_.update(*this, std::move(dualIdSet));
        std::cout << "Initializing of the grid interaction volume index sets took " << timer.elapsed() << " seconds." << std::endl;

        // build the connectivity map for an effecient assembly
        timer.reset();
        connectivityMap_.update(*this);
        std::cout << "Initializing of the connectivity map took " << timer.elapsed() << " seconds." << std::endl;
709
710
    }

Dennis Gläser's avatar
Dennis Gläser committed
711
712
713
    /*!
     * \brief Returns the sub control volume face indices of an scv by global index.
     */
714
715
716
    const std::vector<IndexType>& scvfIndicesOfScv(IndexType scvIdx) const
    { return scvfIndicesOfScv_[scvIdx]; }

Dennis Gläser's avatar
Dennis Gläser committed
717
718
719
    /*!
     * \brief Returns the neighboring vol var indices for each scvf contained in an scv.
     */
720
721
722
723
    const std::vector< std::vector<IndexType> >& neighborVolVarIndices(IndexType scvIdx) const
    { return neighborVolVarIndices_[scvIdx]; }

    /*!
Dennis Gläser's avatar
Dennis Gläser committed
724
725
     * \brief Returns the connectivity map of which dofs have derivatives with respect
     *        to a given dof.
726
     */
Dennis Gläser's avatar
Dennis Gläser committed
727
728
    const ConnectivityMap& connectivityMap() const
    { return connectivityMap_; }
729

Dennis Gläser's avatar
Dennis Gläser committed
730
731
732
733
734
    /*!
     * \brief Returns the grit interaction volume seeds class.
     */
    const GridIVIndexSets& gridInteractionVolumeIndexSets() const
    { return ivIndexSets_; }
735

Dennis Gläser's avatar
Dennis Gläser committed
736
737
738
739
private:
    // mappers
    Dumux::ElementMap<GridView> elementMap_;
    ConnectivityMap connectivityMap_;
740

Dennis Gläser's avatar
Dennis Gläser committed
741
742
743
744
745
746
    // containers storing the global data
    std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
    std::vector< std::vector< std::vector<IndexType> > > neighborVolVarIndices_;
    std::vector<bool> primaryInteractionVolumeVertices_;
    std::vector<bool> secondaryInteractionVolumeVertices_;
    std::vector<bool> isGhostVertex_;
747
748
    std::size_t numScvs_;
    std::size_t numScvf_;
Dennis Gläser's avatar
Dennis Gläser committed
749
    std::size_t numBoundaryScvf_;
750

Dennis Gläser's avatar
Dennis Gläser committed
751
752
    // The grid interaction volume index set
    GridIVIndexSets ivIndexSets_;
753
754
755
};

} // end namespace
756

757
#endif