discretization/cellcentered/mpfa/fvelementgeometry.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_DISCRETIZATION_CCMPFA_FV_ELEMENT_GEOMETRY_HH
#define DUMUX_DISCRETIZATION_CCMPFA_FV_ELEMENT_GEOMETRY_HH
 
#include <optional>
#include <utility>
 
#include <dune/common/exceptions.hh>
#include <dune/common/iteratorrange.hh>
#include <dune/geometry/type.hh>
 
#include <dumux/common/parameters.hh>
#include <dumux/common/indextraits.hh>
#include <dumux/discretization/scvandscvfiterators.hh>
 
namespace Dumux {
 
#ifndef DOXYGEN
namespace Detail::Mpfa {
 
template<typename GridGeometry, typename SubControlVolumeFace>
typename SubControlVolumeFace::Traits::Geometry makeScvfGeometry(const GridGeometry& gridGeometry,
                                                                 const SubControlVolumeFace& scvf)
{
    static constexpr int dim = GridGeometry::GridView::dimension;
 
    const auto& facetInfo = scvf.facetInfo();
    const auto element = gridGeometry.element(facetInfo.elementIndex);
    const auto elemGeo = element.geometry();
    const auto refElement = referenceElement(elemGeo);
    for (const auto& is : intersections(gridGeometry.gridView(), element))
    {
        if (is.indexInInside() == facetInfo.facetIndex)
        {
            const auto numCorners = is.geometry().corners();
            const auto isPositions = GridGeometry::MpfaHelper::computeScvfCornersOnIntersection(
                elemGeo, refElement, facetInfo.facetIndex, numCorners
            );
            return {
                Dune::GeometryTypes::cube(dim-1),
                GridGeometry::MpfaHelper::getScvfCorners(
                    isPositions, numCorners, facetInfo.facetCornerIndex
                )
            };
        }
    }
    DUNE_THROW(Dune::InvalidStateException, "Could not construct scvf geometry");
}
 
template<typename GridGeometry, typename SubControlVolumeFace>
auto getVertexCorner(const GridGeometry& gridGeometry, const SubControlVolumeFace& scvf)
{
    static constexpr int dim = GridGeometry::GridView::dimension;
 
    const auto& facetInfo = scvf.facetInfo();
    const auto element = gridGeometry.element(facetInfo.elementIndex);
    const auto elemGeo = element.geometry();
    const auto refElement = referenceElement(elemGeo);
    return elemGeo.global(refElement.position(
        refElement.subEntity(facetInfo.facetIndex, 1, facetInfo.facetCornerIndex, dim),
        dim
    ));
}
 
template<typename GridGeometry, typename SubControlVolumeFace>
auto getFacetCorner(const GridGeometry& gridGeometry, const SubControlVolumeFace& scvf)
{
    const auto& facetInfo = scvf.facetInfo();
    const auto element = gridGeometry.element(facetInfo.elementIndex);
    const auto elemGeo = element.geometry();
    const auto refElement = referenceElement(elemGeo);
    return elemGeo.global(refElement.position(facetInfo.facetIndex, 1));
}
 
} // namespace Detail::Mpfa
#endif // DOXYGEN
 
template<class GG, bool enableGridGeometryCache>
class CCMpfaFVElementGeometry;
 
template<class GG>
class CCMpfaFVElementGeometry<GG, true>
{
    using ThisType = CCMpfaFVElementGeometry<GG, true>;
    using GridView = typename GG::GridView;
    using GridIndexType = typename IndexTraits<GridView>::GridIndex;
 
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
 
public:
    using Element = typename GridView::template Codim<0>::Entity;
    using SubControlVolume = typename GG::SubControlVolume;
    using SubControlVolumeFace = typename GG::SubControlVolumeFace;
    using GridGeometry = GG;
    static constexpr std::size_t maxNumElementScvs = 1;
    static constexpr std::size_t maxNumElementScvfs = dim == 3 ? 24 : 8;
 
    CCMpfaFVElementGeometry(const GridGeometry& gridGeometry)
    : gridGeometryPtr_(&gridGeometry) {}
 
    const SubControlVolume& scv(GridIndexType scvIdx) const
    {
        return gridGeometry().scv(scvIdx);
    }
 
    const SubControlVolumeFace& scvf(GridIndexType scvfIdx) const
    {
        return gridGeometry().scvf(scvfIdx);
    }
 
    const SubControlVolumeFace& flipScvf(GridIndexType scvfIdx, unsigned int outsideScvIdx = 0) const
    {
        return gridGeometry().flipScvf(scvfIdx, outsideScvIdx);
    }
 
    friend inline Dune::IteratorRange< ScvIterator<SubControlVolume, std::array<GridIndexType, 1>, ThisType> >
    scvs(const CCMpfaFVElementGeometry& fvGeometry)
    {
        using ScvIterator = Dumux::ScvIterator<SubControlVolume, std::array<GridIndexType, 1>, ThisType>;
        return Dune::IteratorRange<ScvIterator>(ScvIterator(fvGeometry.scvIndices_.begin(), fvGeometry),
                                                ScvIterator(fvGeometry.scvIndices_.end(), fvGeometry));
    }
 
    friend inline Dune::IteratorRange< ScvfIterator<SubControlVolumeFace, std::vector<GridIndexType>, ThisType> >
    scvfs(const CCMpfaFVElementGeometry& fvGeometry)
    {
        const auto& g = fvGeometry.gridGeometry();
        const auto scvIdx = fvGeometry.scvIndices_[0];
        using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<GridIndexType>, ThisType>;
        return Dune::IteratorRange<ScvfIterator>(ScvfIterator(g.scvfIndicesOfScv(scvIdx).begin(), fvGeometry),
                                                 ScvfIterator(g.scvfIndicesOfScv(scvIdx).end(), fvGeometry));
    }
 
    std::size_t numScv() const
    {
        return scvIndices_.size();
    }
 
    std::size_t numScvf() const
    {
        return gridGeometry().scvfIndicesOfScv(scvIndices_[0]).size();
    }
 
    CCMpfaFVElementGeometry bind(const Element& element) &&
    {
        this->bindElement(element);
        return std::move(*this);
    }
 
    void bind(const Element& element) &
    {
        this->bindElement(element);
    }
 
    CCMpfaFVElementGeometry bindElement(const Element& element) &&
    {
        this->bindElement(element);
        return std::move(*this);
    }
 
    void bindElement(const Element& element) &
    {
        element_ = element;
        scvIndices_[0] = gridGeometry().elementMapper().index(element);
    }
 
    bool isBound() const
    { return static_cast<bool>(element_); }
 
    const Element& element() const
    { return *element_; }
 
    const GridGeometry& gridGeometry() const
    { return *gridGeometryPtr_; }
 
    bool hasBoundaryScvf() const
    { return gridGeometry().hasBoundaryScvf(scvIndices_[0]); }
 
    typename Element::Geometry geometry(const SubControlVolume& scv) const
    { return gridGeometryPtr_->element(scv.dofIndex()).geometry(); }
 
    typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const
    { return Detail::Mpfa::makeScvfGeometry(gridGeometry(), scvf); }
 
    typename SubControlVolumeFace::Traits::GlobalPosition vertexCorner(const SubControlVolumeFace& scvf) const
    { return Detail::Mpfa::getVertexCorner(gridGeometry(), scvf); }
 
    typename SubControlVolumeFace::Traits::GlobalPosition facetCorner(const SubControlVolumeFace& scvf) const
    { return Detail::Mpfa::getFacetCorner(gridGeometry(), scvf); }
 
private:
 
    std::optional<Element> element_;
    std::array<GridIndexType, 1> scvIndices_;
    const GridGeometry* gridGeometryPtr_;
};
 
template<class GG>
class CCMpfaFVElementGeometry<GG, false>
{
    using ThisType = CCMpfaFVElementGeometry<GG, false>;
    using GridView = typename GG::GridView;
    using GridIndexType = typename IndexTraits<GridView>::GridIndex;
    using MpfaHelper = typename GG::MpfaHelper;
 
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
    using CoordScalar = typename GridView::ctype;
 
public:
    using Element = typename GridView::template Codim<0>::Entity;
    using SubControlVolume = typename GG::SubControlVolume;
    using SubControlVolumeFace = typename GG::SubControlVolumeFace;
    using GridGeometry = GG;
    static constexpr std::size_t maxNumElementScvs = 1;
    static constexpr std::size_t maxNumElementScvfs = dim == 3 ? 24 : 8;
 
    CCMpfaFVElementGeometry(const GridGeometry& gridGeometry)
    : gridGeometryPtr_(&gridGeometry) {}
 
    const SubControlVolume& scv(GridIndexType scvIdx) const
    {
        if (scvIdx == scvIndices_[0])
            return scvs_[0];
        else
            return neighborScvs_[findLocalIndex(scvIdx, neighborScvIndices_)];
    }
 
    const SubControlVolumeFace& scvf(GridIndexType scvfIdx) const
    {
        auto it = std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx);
        if (it != scvfIndices_.end())
            return scvfs_[std::distance(scvfIndices_.begin(), it)];
        else
            return neighborScvfs_[findLocalIndex(scvfIdx, neighborScvfIndices_)];
    }
 
    const SubControlVolumeFace& flipScvf(GridIndexType scvfIdx, unsigned int outsideScvIdx = 0) const
    {
        return scvf( gridGeometry().flipScvfIdx(scvfIdx, outsideScvIdx) );
    }
 
    friend inline Dune::IteratorRange<typename std::array<SubControlVolume, 1>::const_iterator>
    scvs(const ThisType& g)
    {
        using IteratorType = typename std::array<SubControlVolume, 1>::const_iterator;
        return Dune::IteratorRange<IteratorType>(g.scvs_.begin(), g.scvs_.end());
    }
 
    friend inline Dune::IteratorRange<typename std::vector<SubControlVolumeFace>::const_iterator>
    scvfs(const ThisType& g)
    {
        using IteratorType = typename std::vector<SubControlVolumeFace>::const_iterator;
        return Dune::IteratorRange<IteratorType>(g.scvfs_.begin(), g.scvfs_.end());
    }
 
    std::size_t numScv() const
    { return scvs_.size(); }
 
    std::size_t numScvf() const
    { return scvfs_.size(); }
 
    CCMpfaFVElementGeometry bindElement(const Element& element) &&
    {
        this->bindElement_(element);
        return std::move(*this);
    }
 
    void bindElement(const Element& element) &
    {
        this->bindElement_(element);
    }
 
    CCMpfaFVElementGeometry bind(const Element& element) &&
    {
        this->bind_(element);
        return std::move(*this);
    }
 
    void bind(const Element& element) &
    {
        this->bind_(element);
    }
 
    bool isBound() const
    { return static_cast<bool>(element_); }
 
    const Element& element() const
    { return *element_; }
 
    const GridGeometry& gridGeometry() const
    { return *gridGeometryPtr_; }
 
    bool hasBoundaryScvf() const
    { return hasBoundaryScvf_; }
 
    typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const
    { return gridGeometryPtr_->element(scv.dofIndex()).geometry(); }
 
    typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const
    { return Detail::Mpfa::makeScvfGeometry(gridGeometry(), scvf); }
 
    typename SubControlVolumeFace::Traits::GlobalPosition vertexCorner(const SubControlVolumeFace& scvf) const
    { return Detail::Mpfa::getVertexCorner(gridGeometry(), scvf); }
 
    typename SubControlVolumeFace::Traits::GlobalPosition facetCorner(const SubControlVolumeFace& scvf) const
    { return Detail::Mpfa::getFacetCorner(gridGeometry(), scvf); }
 
private:
 
    void bind_(const Element& element)
    {
        // make inside geometries
        bindElement_(element);
 
        // get some references for convenience
        const auto globalI = gridGeometry().elementMapper().index(element);
        const auto& assemblyMapI = gridGeometry().connectivityMap()[globalI];
 
        // reserve memory
        const auto numNeighbors = assemblyMapI.size();
        const auto numNeighborScvfs = numNeighborScvfs_(assemblyMapI);
        neighborScvs_.reserve(numNeighbors);
        neighborScvIndices_.reserve(numNeighbors);
        neighborScvfIndices_.reserve(numNeighborScvfs);
        neighborScvfs_.reserve(numNeighborScvfs);
 
        // make neighbor geometries
        // use the assembly map to determine which faces are necessary
        for (const auto& dataJ : assemblyMapI)
            makeNeighborGeometries(gridGeometry().element(dataJ.globalJ),
                                   dataJ.globalJ,
                                   dataJ.scvfsJ,
                                   dataJ.additionalScvfs);
 
        // //! TODO Check if user added additional DOF dependencies, i.e. the residual of DOF globalI depends
        // //! on additional DOFs not included in the discretization schemes' occupation pattern
        // const auto& additionalDofDependencies = problem.getAdditionalDofDependencies(globalI);
        // if (!additionalDofDependencies.empty())
        // {
        //     const auto newNumNeighbors = neighborScvs_.size() + additionalDofDependencies.size();
        //     neighborScvs_.reserve(newNumNeighbors);
        //     neighborScvIndices_.reserve(newNumNeighbors);
        //     for (auto globalJ : additionalDofDependencies)
        //     {
        //         neighborScvs_.emplace_back(gridGeometry().element(globalJ).geometry(), globalJ);
        //         neighborScvIndices_.emplace_back(globalJ);
        //     }
        // }
    }
 
    void bindElement_(const Element& element)
    {
        clear();
        element_ = element;
        makeElementGeometries(element);
    }
 
    template<class DataJContainer>
    std::size_t numNeighborScvfs_(const DataJContainer& dataJContainer)
    {
        std::size_t numNeighborScvfs = 0;
        for (const auto& dataJ : dataJContainer)
            numNeighborScvfs += dataJ.scvfsJ.size() + dataJ.additionalScvfs.size();
        return numNeighborScvfs;
    }
 
    void makeElementGeometries(const Element& element)
    {
        // make the scv
        const auto eIdx = gridGeometry().elementMapper().index(element);
        scvs_[0] = SubControlVolume(element.geometry(), eIdx);
        scvIndices_[0] = eIdx;
 
        // get data on the scv faces
        const auto& scvFaceIndices = gridGeometry().scvfIndicesOfScv(eIdx);
        const auto& neighborVolVarIndices = gridGeometry().neighborVolVarIndices(eIdx);
 
        // the quadrature point parameterizaion to be used on scvfs
        static const auto q = getParam<CoordScalar>("MPFA.Q");
 
        // reserve memory for the scv faces
        const auto numLocalScvf = scvFaceIndices.size();
        scvfIndices_.reserve(numLocalScvf);
        scvfs_.reserve(numLocalScvf);
 
        // for network grids we only want to do one scvf per half facet
        // this approach assumes conforming grids at branching facets
        std::vector<bool> finishedFacets;
        if (dim < dimWorld)
            finishedFacets.resize(element.subEntities(1), false);
 
        int scvfCounter = 0;
        for (const auto& is : intersections(gridGeometry().gridView(), element))
        {
            // if we are dealing with a lower dimensional network
            // only make a new scvf if we haven't handled it yet
            if (dim < dimWorld)
            {
                const auto indexInInside = is.indexInInside();
                if (finishedFacets[indexInInside])
                    continue;
                else
                    finishedFacets[indexInInside] = true;
            }
 
            // if outside level > inside level, use the outside element in the following
            bool useNeighbor = is.neighbor() && is.outside().level() > element.level();
            const auto& e = useNeighbor ? is.outside() : element;
            const auto indexInElement = useNeighbor ? is.indexInOutside() : is.indexInInside();
            const auto eg = e.geometry();
            const auto refElement = referenceElement(eg);
 
            // 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);
 
            // make the scv faces belonging to each corner of the intersection
            for (int c = 0; c < numCorners; ++c)
            {
                // get the global vertex index the scv face is connected to
                auto vIdxLocal = refElement.subEntity(indexInElement, 1, c, dim);
                auto vIdxGlobal = gridGeometry().vertexMapper().subIndex(e, vIdxLocal, dim);
 
                // do not build scvfs connected to a processor boundary
                if (gridGeometry().isGhostVertex(vIdxGlobal))
                    continue;
 
                hasBoundaryScvf_ = (hasBoundaryScvf_ || is.boundary());
                typename SubControlVolumeFace::FacetInfo facetInfo{
                    gridGeometry().elementMapper().index(e),
                    useNeighbor ? is.indexInOutside() : is.indexInInside(),
                    c
                };
                scvfs_.emplace_back(MpfaHelper(),
                                    MpfaHelper::getScvfCorners(isPositions, numCorners, c),
                                    is,
                                    std::move(facetInfo),
                                    vIdxGlobal,
                                    vIdxLocal,
                                    scvFaceIndices[scvfCounter],
                                    eIdx,
                                    neighborVolVarIndices[scvfCounter],
                                    q,
                                    is.boundary());
 
                scvfIndices_.emplace_back(scvFaceIndices[scvfCounter]);
                scvfCounter++;
            }
        }
    }
 
    template<typename IndexVector>
    void makeNeighborGeometries(const Element& element,
                                GridIndexType eIdxGlobal,
                                const IndexVector& scvfIndices,
                                const IndexVector& additionalScvfs)
    {
        // create the neighbor scv if it doesn't exist yet
        neighborScvs_.emplace_back(element.geometry(), eIdxGlobal);
        neighborScvIndices_.push_back(eIdxGlobal);
 
        // get data on the scv faces
        const auto& scvFaceIndices = gridGeometry().scvfIndicesOfScv(eIdxGlobal);
        const auto& neighborVolVarIndices = gridGeometry().neighborVolVarIndices(eIdxGlobal);
 
        // the quadrature point parameterizaion to be used on scvfs
        static const auto q = getParam<CoordScalar>("MPFA.Q");
 
        // for network grids we only want to do one scvf per half facet
        // this approach assumes conforming grids at branching facets
        std::vector<bool> finishedFacets;
        if (dim < dimWorld)
            finishedFacets.resize(element.subEntities(1), false);
 
        int scvfCounter = 0;
        for (const auto& is : intersections(gridGeometry().gridView(), element))
        {
            // if we are dealing with a lower dimensional network
            // only make a new scvf if we haven't handled it yet
            if (dim < dimWorld)
            {
                auto indexInInside = is.indexInInside();
                if(finishedFacets[indexInInside])
                    continue;
                else
                    finishedFacets[indexInInside] = true;
            }
 
            // if outside level > inside level, use the outside element in the following
            bool useNeighbor = is.neighbor() && is.outside().level() > element.level();
            const auto& e = useNeighbor ? is.outside() : element;
            const auto indexInElement = useNeighbor ? is.indexInOutside() : is.indexInInside();
            const auto eg = e.geometry();
            const auto refElement = referenceElement(eg);
 
            // 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);
 
            // make the scv faces belonging to each corner of the intersection
            for (int c = 0; c < numCorners; ++c)
            {
                // get the global vertex index the scv face is connected to
                auto vIdxLocal = refElement.subEntity(indexInElement, 1, c, dim);
                auto vIdxGlobal = gridGeometry().vertexMapper().subIndex(e, vIdxLocal, dim);
 
                // do not build scvfs connected to a processor boundary
                if (gridGeometry().isGhostVertex(vIdxGlobal))
                    continue;
 
                // only build the scvf if it is in the list of necessary indices
                if (!MpfaHelper::vectorContainsValue(scvfIndices, scvFaceIndices[scvfCounter])
                    && !MpfaHelper::vectorContainsValue(additionalScvfs, scvFaceIndices[scvfCounter]))
                {
                    // increment counter either way
                    scvfCounter++;
                    continue;
                }
 
                // build scvf
                typename SubControlVolumeFace::FacetInfo facetInfo{
                    gridGeometry().elementMapper().index(e),
                    useNeighbor ? is.indexInOutside() : is.indexInInside(),
                    c
                };
                neighborScvfs_.emplace_back(MpfaHelper(),
                                            MpfaHelper::getScvfCorners(isPositions, numCorners, c),
                                            is,
                                            std::move(facetInfo),
                                            vIdxGlobal,
                                            vIdxLocal,
                                            scvFaceIndices[scvfCounter],
                                            eIdxGlobal,
                                            neighborVolVarIndices[scvfCounter],
                                            q,
                                            is.boundary());
 
                neighborScvfIndices_.emplace_back(scvFaceIndices[scvfCounter]);
 
                // increment counter
                scvfCounter++;
            }
        }
    }
 
    unsigned int findLocalIndex(const GridIndexType idx,
                                const std::vector<GridIndexType>& indices) const
    {
        auto it = std::find(indices.begin(), indices.end(), idx);
        assert(it != indices.end() && "Could not find the scv/scvf! Make sure to properly bind this class!");
        return std::distance(indices.begin(), it);
    }
 
    void clear()
    {
        scvfIndices_.clear();
        scvfs_.clear();
 
        neighborScvIndices_.clear();
        neighborScvfIndices_.clear();
        neighborScvs_.clear();
        neighborScvfs_.clear();
 
        hasBoundaryScvf_ = false;
    }
 
    const GridGeometry* gridGeometryPtr_;
    std::optional<Element> element_;
 
    // local storage after binding an element
    std::array<GridIndexType, 1> scvIndices_;
    std::vector<GridIndexType> scvfIndices_;
    std::array<SubControlVolume, 1> scvs_;
    std::vector<SubControlVolumeFace> scvfs_;
 
    std::vector<GridIndexType> neighborScvIndices_;
    std::vector<GridIndexType> neighborScvfIndices_;
    std::vector<SubControlVolume> neighborScvs_;
    std::vector<SubControlVolumeFace> neighborScvfs_;
 
    bool hasBoundaryScvf_ = false;
};
 
} // end namespace
 
#endif