fvgridgeometry.hh 33.1 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
Dennis Gläser's avatar
Dennis Gläser committed
21 22
 * \ingroup CCMpfaDiscretization
 * \brief The finite volume geometry (scvs and scvfs) for cell-centered mpfa models on a grid view
23
 *        This builds up the sub control volumes and sub control volume faces
Dennis Gläser's avatar
Dennis Gläser committed
24
 *        for each element of the grid partition.
25
 */
26 27
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_FV_GRID_GEOMETRY_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_FV_GRID_GEOMETRY_HH
28 29 30

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

Dennis Gläser's avatar
Dennis Gläser committed
32 33 34
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>

Dennis Gläser's avatar
Dennis Gläser committed
35 36 37 38 39
#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>
40 41 42 43

namespace Dumux
{
/*!
Dennis Gläser's avatar
Dennis Gläser committed
44 45
 * \ingroup CCMpfaDiscretization
 * \brief The finite volume geometry (scvs and scvfs) for cell-centered mpfa models on a grid view
46
 *        This builds up the sub control volumes and sub control volume faces
Dennis Gläser's avatar
Dennis Gläser committed
47
 * \note This class is specialized for versions with and without caching the fv geometries on the grid view
48 49
 */
template<class TypeTag, bool EnableFVElementGeometryCache>
Dennis Gläser's avatar
Dennis Gläser committed
50
class CCMpfaFVGridGeometry;
51

Dennis Gläser's avatar
Dennis Gläser committed
52 53 54 55 56 57
/*!
 * \ingroup CCMpfaDiscretization
 * \brief The finite volume geometry (scvs and scvfs) for cell-centered mpfa models on a grid view
 *        This builds up the sub control volumes and sub control volume faces
 * \note For caching enabled we store the fv geometries for the whole grid view which is memory intensive but faster
 */
58
template<class TypeTag>
Dennis Gläser's avatar
Dennis Gläser committed
59
class CCMpfaFVGridGeometry<TypeTag, true> : public BaseFVGridGeometry<TypeTag>
60
{
Dennis Gläser's avatar
Dennis Gläser committed
61 62
    using ParentType = BaseFVGridGeometry<TypeTag>;
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
63
    using ElementMapper = typename GET_PROP_TYPE(TypeTag, ElementMapper);
Dennis Gläser's avatar
Dennis Gläser committed
64

65
    using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
Dennis Gläser's avatar
Dennis Gläser committed
66 67
    using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);

68 69 70 71
    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
72 73 74
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
Dennis Gläser's avatar
Dennis Gläser committed
75

Dennis Gläser's avatar
Dennis Gläser committed
76 77 78
    using Element = typename GridView::template Codim<0>::Entity;
    using Vertex = typename GridView::template Codim<dim>::Entity;
    using Intersection = typename GridView::Intersection;
79
    using CoordScalar = typename GridView::ctype;
Dennis Gläser's avatar
Dennis Gläser committed
80 81 82 83 84
    using GridIndexType = typename GridView::IndexSet::IndexType;
    using LocalIndexType = typename PrimaryInteractionVolume::Traits::LocalIndexType;

    using GridIVIndexSets = CCMpfaGridInteractionVolumeIndexSets<TypeTag>;
    using ConnectivityMap = CCMpfaConnectivityMap<TypeTag>;
85 86 87 88

    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;

public:
Dennis Gläser's avatar
Dennis Gläser committed
89
    using SecondaryIvIndicatorType = std::function<bool(const Element&, const Intersection&, bool)>;
90 91 92

    //! Constructor without indicator function for secondary interaction volumes
    //! Per default, we use the secondary IVs at branching points & boundaries
Dennis Gläser's avatar
Dennis Gläser committed
93 94 95 96 97 98 99 100 101 102 103
    CCMpfaFVGridGeometry(const GridView& gridView)
    : ParentType(gridView)
    , secondaryIvIndicator_([] (const Element& e, const Intersection& is, bool isBranching)
                               { return is.boundary() || isBranching; } )
    {}

    //! Constructor with user-defined indicator function for secondary interaction volumes
    CCMpfaFVGridGeometry(const GridView& gridView, const SecondaryIvIndicatorType& indicator)
    : ParentType(gridView)
    , secondaryIvIndicator_(indicator)
    {}
104

105 106
    //! the element mapper is the dofMapper
    //! this is convenience to have better chance to have the same main files for box/tpfa/mpfa...
Dennis Gläser's avatar
Dennis Gläser committed
107 108 109 110 111 112 113 114 115 116 117 118 119
    const ElementMapper& dofMapper() const { return this->elementMapper(); }

    //! The total number of sub control volumes
    std::size_t numScv() const { return scvs_.size(); }

    //! The total number of sub control volume faces
    std::size_t numScvf() const { return scvfs_.size(); }

    //! The total number of boundary sub control volume faces
    std::size_t numBoundaryScvf() const { return numBoundaryScvf_; }

    //! The total number of degrees of freedom
    std::size_t numDofs() const { return this->gridView().size(0); }
120

Dennis Gläser's avatar
Dennis Gläser committed
121 122 123 124 125 126 127
    //! Get an element from a global element index
    Element element(GridIndexType eIdx) const { return this->elementMap()[eIdx]; }

    //! Get an element from a sub control volume contained in it
    Element element(const SubControlVolume& scv) const { return this->elementMap()[scv.elementIndex()]; }

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

Dennis Gläser's avatar
Dennis Gläser committed
131 132
    //!Returns true if primary interaction volumes are used around a vertex (index).
    bool vertexUsesPrimaryInteractionVolume(GridIndexType vIdxGlobal) const
Dennis Gläser's avatar
Dennis Gläser committed
133
    { return primaryInteractionVolumeVertices_[vIdxGlobal]; }
134

Dennis Gläser's avatar
Dennis Gläser committed
135
    //! Returns if primary interaction volumes are used around a given vertex.
Dennis Gläser's avatar
Dennis Gläser committed
136 137
    bool vertexUsesSecondaryInteractionVolume(const Vertex& v) const
    { return secondaryInteractionVolumeVertices_[this->vertexMapper().index(v)]; }
138

Dennis Gläser's avatar
Dennis Gläser committed
139 140
    //! Returns true if primary interaction volumes are used around a given vertex index.
    bool vertexUsesSecondaryInteractionVolume(GridIndexType vIdxGlobal) const
Dennis Gläser's avatar
Dennis Gläser committed
141
    { return secondaryInteractionVolumeVertices_[vIdxGlobal]; }
142

Dennis Gläser's avatar
Dennis Gläser committed
143
    //! update all fvElementGeometries (do this again after grid adaption)
144
    void update()
145
    {
146 147
        ParentType::update();

Dennis Gläser's avatar
Dennis Gläser committed
148 149
        // stop the time required for the update
        Dune::Timer timer;
150

151
        // clear scvfs container
Dennis Gläser's avatar
Dennis Gläser committed
152 153 154 155
        scvfs_.clear();

        // determine the number of geometric entities
        const auto numVert = this->gridView().size(dim);
156
        const auto numScvs = numDofs();
157
        std::size_t numScvf = MpfaHelper::getGlobalNumScvf(this->gridView());
158

Dennis Gläser's avatar
Dennis Gläser committed
159
        // resize containers
160
        scvs_.resize(numScvs);
Dennis Gläser's avatar
Dennis Gläser committed
161
        scvfs_.reserve(numScvf);
162 163
        scvfIndicesOfScv_.resize(numScvs);

Dennis Gläser's avatar
Dennis Gläser committed
164 165 166 167
        // 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);
168 169

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

Dennis Gläser's avatar
Dennis Gläser committed
172
        // instantiate the dual grid index set (to be used for construction of interaction volumes)
Dennis Gläser's avatar
Dennis Gläser committed
173
        CCMpfaDualGridIndexSet<GridIndexType, LocalIndexType, dim> dualIdSet(this->gridView());
174 175

        // Build the SCVs and SCV faces
Dennis Gläser's avatar
Dennis Gläser committed
176
        GridIndexType scvfIdx = 0;
Dennis Gläser's avatar
Dennis Gläser committed
177 178
        numBoundaryScvf_ = 0;
        for (const auto& element : elements(this->gridView()))
179
        {
Dennis Gläser's avatar
Dennis Gläser committed
180
            const auto eIdx = this->elementMapper().index(element);
181 182 183 184 185

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

            // The local scvf index set
Dennis Gläser's avatar
Dennis Gläser committed
186
            std::vector<GridIndexType> scvfIndexSet;
187 188
            scvfIndexSet.reserve(MpfaHelper::getNumLocalScvfs(elementGeometry.type()));

Dennis Gläser's avatar
Dennis Gläser committed
189
            // for network grids there might be multiple intersections with the same geometryInInside
190
            // we indentify those by the indexInInside for now (assumes conforming grids at branching facets)
Dennis Gläser's avatar
Dennis Gläser committed
191
            std::vector<std::vector<GridIndexType>> outsideIndices;
192 193 194
            if (dim < dimWorld)
            {
                outsideIndices.resize(element.subEntities(1));
Dennis Gläser's avatar
Dennis Gläser committed
195
                for (const auto& intersection : intersections(this->gridView(), element))
196 197 198 199
                {
                    if (intersection.neighbor())
                    {
                        const auto& outside = intersection.outside();
Dennis Gläser's avatar
Dennis Gläser committed
200
                        auto nIdx = this->elementMapper().index(outside);
201 202 203 204 205 206
                        outsideIndices[intersection.indexInInside()].push_back(nIdx);
                    }
                }
            }

            // construct the sub control volume faces
Dennis Gläser's avatar
Dennis Gläser committed
207
            for (const auto& is : intersections(this->gridView(), element))
208 209 210 211 212
            {
                const auto indexInInside = is.indexInInside();
                const bool boundary = is.boundary();
                const bool neighbor = is.neighbor();

Dennis Gläser's avatar
Dennis Gläser committed
213
                // for surface grids, skip the rest if handled already
214 215 216 217 218 219 220 221 222 223
                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
224 225 226 227 228 229 230 231 232
                // 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;
233
                const bool usesSecondaryIV = secondaryIvIndicator_(element, is, isBranchingPoint);
234 235

                // make the scv faces belonging to each corner of the intersection
236
                for (std::size_t c = 0; c < numCorners; ++c)
237 238 239
                {
                    // 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
240
                    const auto vIdxGlobal = this->vertexMapper().subIndex(e, vIdxLocal, dim);
241 242 243 244 245

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

Dennis Gläser's avatar
Dennis Gläser committed
246 247
                    // if this vertex is tagged to use the secondary IVs, store info
                    if (usesSecondaryIV)
248
                    {
Dennis Gläser's avatar
Dennis Gläser committed
249 250
                        primaryInteractionVolumeVertices_[vIdxGlobal] = false;
                        secondaryInteractionVolumeVertices_[vIdxGlobal] = true;
251 252
                    }

Dennis Gläser's avatar
Dennis Gläser committed
253
                    // the quadrature point parameterizarion to be used on scvfs
254
                    static const Scalar q = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Mpfa.Q");
255

Dennis Gläser's avatar
Dennis Gläser committed
256 257 258 259 260 261
                    // make the scv face (for non-boundary scvfs on network grids, use precalculated outside indices)
                    const auto& outsideScvIndices = [&] ()
                                                    {
                                                        if (!boundary)
                                                        {
                                                            return dim == dimWorld ?
Dennis Gläser's avatar
Dennis Gläser committed
262
                                                                   std::vector<GridIndexType>({this->elementMapper().index(is.outside())}) :
Dennis Gläser's avatar
Dennis Gläser committed
263 264 265 266 267
                                                                   outsideIndices[indexInInside];
                                                        }
                                                        else
                                                        {
                                                            // compute boundary scv idx and increment counter
Dennis Gläser's avatar
Dennis Gläser committed
268
                                                            const GridIndexType bIdx = numScvs + numBoundaryScvf_;
Dennis Gläser's avatar
Dennis Gläser committed
269
                                                            numBoundaryScvf_++;
Dennis Gläser's avatar
Dennis Gläser committed
270
                                                            return std::vector<GridIndexType>(1, bIdx);
Dennis Gläser's avatar
Dennis Gläser committed
271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287
                                                        }
                                                    } ();

                    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());
288 289 290 291 292 293

                    // increment scvf counter
                    scvfIdx++;
                }

                // for network grids, clear outside indices to not make a second scvf on that facet
294 295
                if (dim < dimWorld)
                    outsideIndices[indexInInside].clear();
296 297 298 299 300 301 302 303 304 305 306 307 308
            }

            // 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());
309
            for (const auto& scvf : scvfs_)
310 311 312 313 314 315 316 317 318
            {
                if (scvf.boundary())
                    continue;

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

                flipScvfIndices_[scvf.index()].resize(numOutsideScvs);
319
                for (std::size_t i = 0; i < numOutsideScvs; ++i)
320 321 322 323 324 325
                {
                    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
326
                            MpfaHelper::vectorContainsValue(outsideScvf.outsideScvIndices(), insideScvIdx))
327 328 329 330 331 332 333 334 335 336 337
                        {
                            flipScvfIndices_[scvf.index()][i] = outsideScvfIndex;
                            // there is always only one flip face in an outside element
                            break;
                        }
                    }

                }
            }
        }

338
        // building the geometries has finished
Dennis Gläser's avatar
Dennis Gläser committed
339 340 341 342 343 344 345
        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;

346
        // build the connectivity map for an efficient assembly
Dennis Gläser's avatar
Dennis Gläser committed
347 348 349
        timer.reset();
        connectivityMap_.update(*this);
        std::cout << "Initializing of the connectivity map took " << timer.elapsed() << " seconds." << std::endl;
350 351
    }

Dennis Gläser's avatar
Dennis Gläser committed
352 353
    //! Get a sub control volume with a global scv index
    const SubControlVolume& scv(GridIndexType scvIdx) const { return scvs_[scvIdx]; }
354

Dennis Gläser's avatar
Dennis Gläser committed
355 356 357 358 359 360 361 362 363
    //! Get a sub control volume face with a global scvf index
    const SubControlVolumeFace& scvf(GridIndexType scvfIdx) const { return scvfs_[scvfIdx]; }

    //! Returns the connectivity map of which dofs
    //! have derivatives with respect to a given dof.
    const ConnectivityMap& connectivityMap() const { return connectivityMap_; }

    //! Returns the grid interaction volume seeds class.
    const GridIVIndexSets& gridInteractionVolumeIndexSets() const { return ivIndexSets_; }
364

Dennis Gläser's avatar
Dennis Gläser committed
365 366 367 368
    //! Get the scvf on the same face but from the other side
    //! Note that e.g. the normals might be different in the case of surface grids
    const SubControlVolumeFace& flipScvf(GridIndexType scvfIdx, unsigned int outsideScvfIdx = 0) const
    { return scvfs_[flipScvfIndices_[scvfIdx][outsideScvfIdx]]; }
369

Dennis Gläser's avatar
Dennis Gläser committed
370 371 372
    //! Get the sub control volume face indices of an scv by global index
    const std::vector<GridIndexType>& scvfIndicesOfScv(GridIndexType scvIdx) const
    { return scvfIndicesOfScv_[scvIdx]; }
373

Dennis Gläser's avatar
Dennis Gläser committed
374
private:
375
    // connectivity map for efficient assembly
Dennis Gläser's avatar
Dennis Gläser committed
376 377 378
    ConnectivityMap connectivityMap_;

    // the finite volume grid geometries
379 380 381
    std::vector<SubControlVolume> scvs_;
    std::vector<SubControlVolumeFace> scvfs_;

Dennis Gläser's avatar
Dennis Gläser committed
382
    // containers storing the global data
Dennis Gläser's avatar
Dennis Gläser committed
383
    std::vector<std::vector<GridIndexType>> scvfIndicesOfScv_;
Dennis Gläser's avatar
Dennis Gläser committed
384 385 386 387
    std::vector<bool> primaryInteractionVolumeVertices_;
    std::vector<bool> secondaryInteractionVolumeVertices_;
    std::size_t numBoundaryScvf_;

388
    // needed for embedded surface and network grids (dim < dimWorld)
Dennis Gläser's avatar
Dennis Gläser committed
389
    std::vector<std::vector<GridIndexType>> flipScvfIndices_;
390

Dennis Gläser's avatar
Dennis Gläser committed
391 392
    // The grid interaction volume index set
    GridIVIndexSets ivIndexSets_;
393 394

    // Indicator function on where to use the secondary IVs
Dennis Gläser's avatar
Dennis Gläser committed
395
    SecondaryIvIndicatorType secondaryIvIndicator_;
396 397
};

Dennis Gläser's avatar
Dennis Gläser committed
398 399 400 401 402 403 404
/*!
 * \ingroup CCMpfaDiscretization
 * \brief The finite volume geometry (scvs and scvfs) for cell-centered mpfa models on a grid view
 *        This builds up the sub control volumes and sub control volume faces
 * \note For caching disabled we store only some essential index maps to build up local systems on-demand in
 *       the corresponding FVElementGeometry
 */
405
template<class TypeTag>
Dennis Gläser's avatar
Dennis Gläser committed
406
class CCMpfaFVGridGeometry<TypeTag, false> : public BaseFVGridGeometry<TypeTag>
407
{
Dennis Gläser's avatar
Dennis Gläser committed
408 409
    using ParentType = BaseFVGridGeometry<TypeTag>;
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
Dennis Gläser's avatar
Dennis Gläser committed
410 411
    using ElementMapper = typename GET_PROP_TYPE(TypeTag, ElementMapper);

412
    using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
Dennis Gläser's avatar
Dennis Gläser committed
413 414
    using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);

415 416 417 418
    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
419
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
420 421
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
Dennis Gläser's avatar
Dennis Gläser committed
422

Dennis Gläser's avatar
Dennis Gläser committed
423 424 425 426
    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;
Dennis Gläser's avatar
Dennis Gläser committed
427 428 429 430 431
    using GridIndexType = typename GridView::IndexSet::IndexType;
    using LocalIndexType = typename PrimaryInteractionVolume::Traits::LocalIndexType;

    using GridIVIndexSets = CCMpfaGridInteractionVolumeIndexSets<TypeTag>;
    using ConnectivityMap = CCMpfaConnectivityMap<TypeTag>;
432 433 434 435

    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;

public:
436 437 438 439
    using SecondaryIvIndicator = std::function<bool(const Element&, const Intersection&, bool)>;

    //! Constructor without indicator function for secondary interaction volumes
    //! Per default, we use the secondary IVs at branching points & boundaries
Dennis Gläser's avatar
Dennis Gläser committed
440 441 442 443 444 445 446 447 448 449 450
    CCMpfaFVGridGeometry(const GridView& gridView)
    : ParentType(gridView)
    , secondaryIvIndicator_([] (const Element& e, const Intersection& is, bool isBranching)
                               { return is.boundary() || isBranching; } )
    {}

    //! Constructor with user-defined indicator function for secondary interaction volumes
    CCMpfaFVGridGeometry(const GridView& gridView, const SecondaryIvIndicator& indicator)
    : ParentType(gridView)
    , secondaryIvIndicator_(indicator)
    {}
451

452 453 454 455 456
    //! the element mapper is the dofMapper
    //! this is convenience to have better chance to have the same main files for box/tpfa/mpfa...
    const ElementMapper& dofMapper() const
    { return this->elementMapper(); }

Dennis Gläser's avatar
Dennis Gläser committed
457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475
    //! Returns the total number of sub control volumes.
    std::size_t numScv() const { return numScvs_; }

    //! Returns the total number of sub control volume faces.
    std::size_t numScvf() const { return numScvf_; }

    //! Returns the number of scvfs on the domain boundary.
    std::size_t numBoundaryScvf() const { return numBoundaryScvf_; }

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

    //! Gets an element from a global element index.
    Element element(GridIndexType eIdx) const { return this->elementMap()[eIdx]; }

    //! Gets an element from a sub control volume contained in it.
    Element element(const SubControlVolume& scv) const { return this->elementMap()[scv.elementIndex()]; }

    //! Returns true if primary interaction volumes are used around a given vertex.
Dennis Gläser's avatar
Dennis Gläser committed
476 477
    bool vertexUsesPrimaryInteractionVolume(const Vertex& v) const
    { return primaryInteractionVolumeVertices_[this->vertexMapper().index(v)]; }
478

Dennis Gläser's avatar
Dennis Gläser committed
479 480
    //! Returns true if primary interaction volumes are used around a given vertex (index).
    bool vertexUsesPrimaryInteractionVolume(GridIndexType vIdxGlobal) const
Dennis Gläser's avatar
Dennis Gläser committed
481
    { return primaryInteractionVolumeVertices_[vIdxGlobal]; }
482

Dennis Gläser's avatar
Dennis Gläser committed
483
    //! Returns if primary interaction volumes are used around a given vertex.
Dennis Gläser's avatar
Dennis Gläser committed
484 485
    bool vertexUsesSecondaryInteractionVolume(const Vertex& v) const
    { return secondaryInteractionVolumeVertices_[this->vertexMapper().index(v)]; }
486

Dennis Gläser's avatar
Dennis Gläser committed
487 488
    //! Returns true if primary interaction volumes are used around a given vertex (index).
    bool vertexUsesSecondaryInteractionVolume(GridIndexType vIdxGlobal) const
Dennis Gläser's avatar
Dennis Gläser committed
489
    { return secondaryInteractionVolumeVertices_[vIdxGlobal]; }
490

Dennis Gläser's avatar
Dennis Gläser committed
491 492 493 494 495 496 497
    //! 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)]; }

    //! Returns true if the vertex (index) lies on a processor boundary inside a ghost element.
    bool isGhostVertex(GridIndexType vIdxGlobal) const { return isGhostVertex_[vIdxGlobal]; }

    //! Updates all finite volume geometries of the grid. Has to be called again after grid adaption.
498
    void update()
499
    {
500 501
        ParentType::update();

Dennis Gläser's avatar
Dennis Gläser committed
502 503
        // stop the time required for the update
        Dune::Timer timer;
504

Dennis Gläser's avatar
Dennis Gläser committed
505
        // resize containers
506
        numScvs_ = numDofs();
507 508 509
        scvfIndicesOfScv_.resize(numScvs_);
        neighborVolVarIndices_.resize(numScvs_);

Dennis Gläser's avatar
Dennis Gläser committed
510 511 512 513 514
        // 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);
515

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

Dennis Gläser's avatar
Dennis Gläser committed
519
        // instantiate the dual grid index set (to be used for construction of interaction volumes)
Dennis Gläser's avatar
Dennis Gläser committed
520
        CCMpfaDualGridIndexSet<GridIndexType, LocalIndexType, dim> dualIdSet(this->gridView());
521

Dennis Gläser's avatar
Dennis Gläser committed
522
        // Build the SCVs and SCV faces
523
        numScvf_ = 0;
Dennis Gläser's avatar
Dennis Gläser committed
524 525
        numBoundaryScvf_ = 0;
        for (const auto& element : elements(this->gridView()))
526
        {
Dennis Gläser's avatar
Dennis Gläser committed
527 528 529 530
            const auto eIdx = this->elementMapper().index(element);

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

            // the element-wise index sets for finite volume geometry
Dennis Gläser's avatar
Dennis Gläser committed
533
            const auto numLocalFaces = MpfaHelper::getNumLocalScvfs(elementGeometry.type());
Dennis Gläser's avatar
Dennis Gläser committed
534 535
            std::vector<GridIndexType> scvfsIndexSet;
            std::vector< std::vector<GridIndexType> > neighborVolVarIndexSet;
536 537 538
            scvfsIndexSet.reserve(numLocalFaces);
            neighborVolVarIndexSet.reserve(numLocalFaces);

Dennis Gläser's avatar
Dennis Gläser committed
539
            // for network grids there might be multiple intersections with the same geometryInInside
540
            // we indentify those by the indexInInside for now (assumes conforming grids at branching facets)
Dennis Gläser's avatar
Dennis Gläser committed
541
            std::vector<std::vector<GridIndexType>> outsideIndices;
542 543 544
            if (dim < dimWorld)
            {
                outsideIndices.resize(element.subEntities(1));
Dennis Gläser's avatar
Dennis Gläser committed
545
                for (const auto& intersection : intersections(this->gridView(), element))
546 547 548 549
                {
                    if (intersection.neighbor())
                    {
                        const auto& outside = intersection.outside();
Dennis Gläser's avatar
Dennis Gläser committed
550
                        auto nIdx = this->elementMapper().index(outside);
551 552 553 554 555 556
                        outsideIndices[intersection.indexInInside()].push_back(nIdx);
                    }
                }
            }

            // construct the sub control volume faces
Dennis Gläser's avatar
Dennis Gläser committed
557
            for (const auto& is : intersections(this->gridView(), element))
558 559 560 561 562
            {
                const auto indexInInside = is.indexInInside();
                const bool boundary = is.boundary();
                const bool neighbor = is.neighbor();

Dennis Gläser's avatar
Dennis Gläser committed
563
                // for surface grids, skip the rest if handled already
564 565 566 567 568 569 570 571 572 573
                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
574 575
                // evaluate if vertices on this intersection use primary/secondary IVs
                const bool isBranchingPoint = dim < dimWorld ? outsideIndices[indexInInside].size() > 1 : false;
576
                const bool usesSecondaryIV = secondaryIvIndicator_(element, is, isBranchingPoint);
Dennis Gläser's avatar
Dennis Gläser committed
577 578

                // make the scv faces belonging to each corner of the intersection
579
                for (std::size_t c = 0; c < is.geometry().corners(); ++c)
580
                {
Dennis Gläser's avatar
Dennis Gläser committed
581
                    // get the global vertex index the scv face is connected to
582
                    const auto vIdxLocal = refElement.subEntity(indexInElement, 1, c, dim);
Dennis Gläser's avatar
Dennis Gläser committed
583
                    const auto vIdxGlobal = this->vertexMapper().subIndex(e, vIdxLocal, dim);
584 585

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

Dennis Gläser's avatar
Dennis Gläser committed
589 590
                    // if this vertex is tagged to use the secondary IVs, store info
                    if (usesSecondaryIV)
591
                    {
Dennis Gläser's avatar
Dennis Gläser committed
592 593
                        primaryInteractionVolumeVertices_[vIdxGlobal] = false;
                        secondaryInteractionVolumeVertices_[vIdxGlobal] = true;
594 595
                    }

Dennis Gläser's avatar
Dennis Gläser committed
596 597 598 599 600 601
                    // make the scv face (for non-boundary scvfs on network grids, use precalculated outside indices)
                    const auto& outsideScvIndices = [&] ()
                                                    {
                                                        if (!boundary)
                                                        {
                                                            return dim == dimWorld ?
Dennis Gläser's avatar
Dennis Gläser committed
602
                                                                   std::vector<GridIndexType>({this->elementMapper().index(is.outside())}) :
Dennis Gläser's avatar
Dennis Gläser committed
603 604 605 606 607
                                                                   outsideIndices[indexInInside];
                                                        }
                                                        else
                                                        {
                                                            // compute boundary scv idx and increment counter
Dennis Gläser's avatar
Dennis Gläser committed
608
                                                            const GridIndexType bIdx = numScvs_ + numBoundaryScvf_;
Dennis Gläser's avatar
Dennis Gläser committed
609
                                                            numBoundaryScvf_++;
Dennis Gläser's avatar
Dennis Gläser committed
610
                                                            return std::vector<GridIndexType>(1, bIdx);
Dennis Gläser's avatar
Dennis Gläser committed
611 612 613 614 615 616 617 618 619
                                                        }
                                                    } ();

                    // 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));
620 621 622
                }

                // for network grids, clear outside indices to not make a second scvf on that facet
623 624
                if (dim < dimWorld)
                    outsideIndices[indexInInside].clear();
625 626 627 628 629 630 631
            }

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

632
        // building the geometries has finished
Dennis Gläser's avatar
Dennis Gläser committed
633 634 635 636 637 638 639 640 641 642 643
        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;
644 645
    }

Dennis Gläser's avatar
Dennis Gläser committed
646 647
    //! Returns the sub control volume face indices of an scv by global index.
    const std::vector<GridIndexType>& scvfIndicesOfScv(GridIndexType scvIdx) const
648 649
    { return scvfIndicesOfScv_[scvIdx]; }

Dennis Gläser's avatar
Dennis Gläser committed
650 651
    //! Returns the neighboring vol var indices for each scvf contained in an scv.
    const std::vector< std::vector<GridIndexType> >& neighborVolVarIndices(GridIndexType scvIdx) const
652 653
    { return neighborVolVarIndices_[scvIdx]; }

Dennis Gläser's avatar
Dennis Gläser committed
654 655 656
    //! Returns the connectivity map of which dofs
    //! have derivatives with respect to a given dof.
    const ConnectivityMap& connectivityMap() const { return connectivityMap_; }
657

Dennis Gläser's avatar
Dennis Gläser committed
658 659
    //! Returns the grid interaction volume seeds class.
    const GridIVIndexSets& gridInteractionVolumeIndexSets() const { return ivIndexSets_; }
660

Dennis Gläser's avatar
Dennis Gläser committed
661
private:
662
    // connectivity map for efficient assembly
Dennis Gläser's avatar
Dennis Gläser committed
663
    ConnectivityMap connectivityMap_;
664

Dennis Gläser's avatar
Dennis Gläser committed
665
    // containers storing the global data
Dennis Gläser's avatar
Dennis Gläser committed
666 667
    std::vector<std::vector<GridIndexType>> scvfIndicesOfScv_;
    std::vector< std::vector< std::vector<GridIndexType> > > neighborVolVarIndices_;
Dennis Gläser's avatar
Dennis Gläser committed
668 669 670
    std::vector<bool> primaryInteractionVolumeVertices_;
    std::vector<bool> secondaryInteractionVolumeVertices_;
    std::vector<bool> isGhostVertex_;
671 672
    std::size_t numScvs_;
    std::size_t numScvf_;
Dennis Gläser's avatar
Dennis Gläser committed
673
    std::size_t numBoundaryScvf_;
674

Dennis Gläser's avatar
Dennis Gläser committed
675 676
    // The grid interaction volume index set
    GridIVIndexSets ivIndexSets_;
677 678 679

    // Indicator function on where to use the secondary IVs
    SecondaryIvIndicator secondaryIvIndicator_;
680 681
};

Dennis Gläser's avatar
Dennis Gläser committed
682
} // end namespace Dumux
683

684
#endif