multidomain/facet/box/fvgridgeometry.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-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_FACETCOUPLING_BOX_GRID_FVGEOMETRY_HH
#define DUMUX_FACETCOUPLING_BOX_GRID_FVGEOMETRY_HH
 
#include <ranges>
#include <span>
#include <algorithm>
#include <utility>
 
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/localfunctions/lagrange/lagrangelfecache.hh>
 
#include <dumux/common/indextraits.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/extrusion.hh>
#include <dumux/discretization/basegridgeometry.hh>
#include <dumux/discretization/box/boxgeometryhelper.hh>
#include <dumux/discretization/box/subcontrolvolume.hh>
 
#include <dumux/multidomain/facet/box/fvelementgeometry.hh>
#include <dumux/multidomain/facet/box/subcontrolvolumeface.hh>
#include <dumux/multidomain/facet/vertexmapper.hh>
 
namespace Dumux {
 
namespace Detail {
template<class GV, class T>
using BoxFacetCouplingGeometryHelper_t = Dune::Std::detected_or_t<
    Dumux::BoxGeometryHelper<GV, GV::dimension, typename T::SubControlVolume, typename T::SubControlVolumeFace>,
    SpecifiesGeometryHelper,
    T
>;
} // end namespace Detail

template<class GridView>
struct BoxFacetCouplingDefaultGridGeometryTraits
{
    // use a specialized version of the box scvf
    using SubControlVolume = BoxSubControlVolume<GridView>;
    using SubControlVolumeFace = BoxFacetCouplingSubControlVolumeFace<GridView>;
 
    template<class GridGeometry, bool enableCache>
    using LocalView = BoxFacetCouplingFVElementGeometry<GridGeometry, enableCache>;
 
    // per default we use an mcmg mapper for the elements
    using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
    // the default vertex mapper is the enriched vertex dof mapper
    using VertexMapper = EnrichedVertexDofMapper<GridView>;
    // Mapper type for mapping edges
    using FacetMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
};

template<class Scalar,
         class GridView,
         bool enableGridGeometryCache = false,
         class Traits = BoxFacetCouplingDefaultGridGeometryTraits<GridView> >
class BoxFacetCouplingFVGridGeometry;

template<class Scalar, class GV, class Traits>
class BoxFacetCouplingFVGridGeometry<Scalar, GV, true, Traits>
: public BaseGridGeometry<GV, Traits>
{
    using ThisType = BoxFacetCouplingFVGridGeometry<Scalar, GV, true, Traits>;
    using ParentType = BaseGridGeometry<GV, Traits>;
    using GridIndexType = typename IndexTraits<GV>::GridIndex;
    using LocalIndexType = typename IndexTraits<GV>::LocalIndex;
 
    using Element = typename GV::template Codim<0>::Entity;
    using CoordScalar = typename GV::ctype;
    static const int dim = GV::dimension;
    static const int dimWorld = GV::dimensionworld;
 
public:
    using DiscretizationMethod = DiscretizationMethods::Box;
    static constexpr DiscretizationMethod discMethod{};

    using LocalView = typename Traits::template LocalView<ThisType, true>;
    using SubControlVolume = typename Traits::SubControlVolume;
    using SubControlVolumeFace = typename Traits::SubControlVolumeFace;
    using Extrusion = Extrusion_t<Traits>;
    using DofMapper = typename Traits::VertexMapper;
    using FeCache = Dune::LagrangeLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
    using GridView = GV;

    template<class FacetGridView, class CodimOneGridAdapter>
    BoxFacetCouplingFVGridGeometry(const GridView& gridView,
                                   const FacetGridView& facetGridView,
                                   const CodimOneGridAdapter& codimOneGridAdapter,
                                   bool verbose = false)
    : ParentType(gridView)
    , cache_(*this)
    {
        update_(facetGridView, codimOneGridAdapter, verbose);
    }

    const DofMapper& dofMapper() const
    { return this->vertexMapper(); }

    std::size_t numScv() const
    { return numScv_; }

    std::size_t numScvf() const
    { return numScvf_; }

    std::size_t numBoundaryScvf() const
    { return numBoundaryScvf_; }

    std::size_t numDofs() const
    { return this->vertexMapper().size(); }

    template<class FacetGridView, class CodimOneGridAdapter>
    void update(const GridView& gridView,
                const FacetGridView& facetGridView,
                const CodimOneGridAdapter& codimOneGridAdapter,
                bool verbose = false)
    {
        ParentType::update(gridView);
        update_(facetGridView, codimOneGridAdapter, verbose);
    }

    template<class FacetGridView, class CodimOneGridAdapter>
    void update(GridView&& gridView,
                const FacetGridView& facetGridView,
                const CodimOneGridAdapter& codimOneGridAdapter,
                bool verbose = false)
    {
        ParentType::update(std::move(gridView));
        update_(facetGridView, codimOneGridAdapter, verbose);
    }

    const FeCache& feCache() const
    { return feCache_; }

    bool dofOnBoundary(GridIndexType dofIdx) const
    { return boundaryDofIndices_[dofIdx]; }

    bool dofOnInteriorBoundary(GridIndexType dofIdx) const
    { return interiorBoundaryDofIndices_[dofIdx]; }

    bool dofOnPeriodicBoundary(GridIndexType dofIdx) const
    { return false; }

    [[deprecated("Will be removed after release 3.11. Use periodicallyMappedDofs, returning a range of dofs")]]
    GridIndexType periodicallyMappedDof(GridIndexType dofIdx) const
    { DUNE_THROW(Dune::InvalidStateException, "Periodic boundaries are not supported by the box facet coupling scheme"); }

    std::span<const GridIndexType> periodicallyMappedDofs(GridIndexType dofIdx) const
    { DUNE_THROW(Dune::InvalidStateException, "Periodic boundaries are not supported by the box facet coupling scheme"); }

    friend inline LocalView localView(const BoxFacetCouplingFVGridGeometry& gg)
    { return { gg.cache_ }; }
 
private:
    class BoxFacetCouplingGridGeometryCache
    {
        friend class BoxFacetCouplingFVGridGeometry;
    public:
        using GeometryHelper = Detail::BoxFacetCouplingGeometryHelper_t<GV, Traits>;
        using DofHelper = typename GeometryHelper::DofHelper;
 
        explicit BoxFacetCouplingGridGeometryCache(const BoxFacetCouplingFVGridGeometry& gg)
        : gridGeometry_(&gg)
        {}
 
        const BoxFacetCouplingFVGridGeometry& gridGeometry() const
        { return *gridGeometry_; }

        const std::vector<SubControlVolume>& scvs(GridIndexType eIdx) const
        { return scvs_[eIdx]; }

        const std::vector<SubControlVolumeFace>& scvfs(GridIndexType eIdx) const
        { return scvfs_[eIdx]; }
 
 
    private:
        void clear_()
        {
            scvs_.clear();
            scvfs_.clear();
        }
 
        std::vector<std::vector<SubControlVolume>> scvs_;
        std::vector<std::vector<SubControlVolumeFace>> scvfs_;
 
        const BoxFacetCouplingFVGridGeometry* gridGeometry_;
    };
 
public:
    using Cache = BoxFacetCouplingGridGeometryCache;
 
private:
    using GeometryHelper = typename Cache::GeometryHelper;
 
    template<class FacetGridView, class CodimOneGridAdapter>
    void update_(const FacetGridView& facetGridView,
                 const CodimOneGridAdapter& codimOneGridAdapter,
                 bool verbose = false)
    {
        // enrich the vertex mapper subject to the provided facet grid
        this->vertexMapper().enrich(facetGridView, codimOneGridAdapter, verbose);
 
        cache_.clear_();
 
        // resize containers
        const auto numDof = numDofs();
        const auto numElements = this->gridView().size(0);
        cache_.scvs_.resize(numElements);
        cache_.scvfs_.resize(numElements);
        boundaryDofIndices_.assign(numDof, false);
        interiorBoundaryDofIndices_.assign(numDof, false);
 
        // Build the SCV and SCV faces
        numScv_ = 0;
        numScvf_ = 0;
        numBoundaryScvf_ = 0;
        for (const auto& element : elements(this->gridView()))
        {
            auto eIdx = this->elementMapper().index(element);
 
            // keep track of number of scvs and scvfs
            numScv_ += element.subEntities(dim);
            numScvf_ += element.subEntities(dim-1);
 
            // get the element geometry
            auto elementGeometry = element.geometry();
            const auto refElement = referenceElement(elementGeometry);
 
            // instantiate the geometry helper
            GeometryHelper geometryHelper(elementGeometry);
 
            // construct the sub control volumes
            cache_.scvs_[eIdx].resize(elementGeometry.corners());
            for (LocalIndexType scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx)
                cache_.scvs_[eIdx][scvLocalIdx] = SubControlVolume(geometryHelper.getScvCorners(scvLocalIdx),
                                                                   scvLocalIdx,
                                                                   eIdx,
                                                                   this->vertexMapper().subIndex(element, scvLocalIdx, dim));
 
            // construct the sub control volume faces
            LocalIndexType scvfLocalIdx = 0;
            cache_.scvfs_[eIdx].resize(element.subEntities(dim-1));
            for (; scvfLocalIdx < element.subEntities(dim-1); ++scvfLocalIdx)
            {
                // find the global and local scv indices this scvf is belonging to
                std::array<LocalIndexType, 2> localScvIndices{{
                    static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
                    static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 1, dim))
                }};
 
                // create the sub-control volume face
                cache_.scvfs_[eIdx][scvfLocalIdx] = SubControlVolumeFace(geometryHelper,
                                                                         element,
                                                                         scvfLocalIdx,
                                                                         std::move(localScvIndices));
            }
 
            // construct the sub control volume faces on the domain/interior boundaries
            // skip handled facets (necessary for e.g. Dune::FoamGrid)
            std::vector<unsigned int> handledFacets;
            for (const auto& intersection : intersections(this->gridView(), element))
            {
                if (std::count(handledFacets.begin(), handledFacets.end(), intersection.indexInInside()))
                    continue;
 
                handledFacets.push_back(intersection.indexInInside());
 
                // determine if all corners live on the facet grid
                const auto isGeometry = intersection.geometry();
                const auto numFaceCorners = isGeometry.corners();
                const auto idxInInside = intersection.indexInInside();
                const auto boundary = intersection.boundary();
 
                std::vector<LocalIndexType> vIndicesLocal(numFaceCorners);
                for (int i = 0; i < numFaceCorners; ++i)
                    vIndicesLocal[i] = static_cast<LocalIndexType>(refElement.subEntity(idxInInside, 1, i, dim));
 
                std::vector<LocalIndexType> gridVertexIndices(numFaceCorners);
                for (int i = 0; i < numFaceCorners; ++i)
                    gridVertexIndices[i] = this->vertexMapper().vertexIndex(element, vIndicesLocal[i], dim);
 
                // if the vertices compose a facet element, the intersection is on facet grid
                const bool isOnFacet = codimOneGridAdapter.composeFacetElement(gridVertexIndices);
 
                // make sure there are no periodic boundaries
                if (boundary && intersection.neighbor())
                    DUNE_THROW(Dune::InvalidStateException, "Periodic boundaries are not supported by the box facet coupling scheme");
 
                if (isOnFacet || boundary)
                {
                    // keep track of number of faces
                    numScvf_ += numFaceCorners;
                    numBoundaryScvf_ += int(boundary)*numFaceCorners;
 
                    for (unsigned int isScvfLocalIdx = 0; isScvfLocalIdx < numFaceCorners; ++isScvfLocalIdx)
                    {
                        // find the inside scv this scvf is belonging to (localIdx = element local vertex index)
                        std::array<LocalIndexType, 2> localScvIndices{{
                            vIndicesLocal[isScvfLocalIdx], vIndicesLocal[isScvfLocalIdx]
                        }};
 
                        // create the sub-control volume face
                        cache_.scvfs_[eIdx].emplace_back(geometryHelper,
                                                         intersection,
                                                         isScvfLocalIdx,
                                                         scvfLocalIdx,
                                                         std::move(localScvIndices),
                                                         boundary,
                                                         isOnFacet);
 
                        // Mark vertices to be on domain and/or interior boundary
                        const auto dofIndex = this->vertexMapper().subIndex(element, vIndicesLocal[isScvfLocalIdx], dim);
                        if (boundary) boundaryDofIndices_[ dofIndex ] = boundary && !isOnFacet;
                        if (isOnFacet) interiorBoundaryDofIndices_[ dofIndex ] = isOnFacet;
 
                        // increment local counter
                        scvfLocalIdx++;
                    }
                }
            }
        }
    }
 
    const FeCache feCache_;
 
    // TODO do we need those?
    std::size_t numScv_;
    std::size_t numScvf_;
    std::size_t numBoundaryScvf_;
 
    // vertices on domain/interior boundaries
    std::vector<bool> boundaryDofIndices_;
    std::vector<bool> interiorBoundaryDofIndices_;
 
    Cache cache_;
};

template<class Scalar, class GV, class Traits>
class BoxFacetCouplingFVGridGeometry<Scalar, GV, false, Traits>
: public BaseGridGeometry<GV, Traits>
{
    using ThisType = BoxFacetCouplingFVGridGeometry<Scalar, GV, false, Traits>;
    using ParentType = BaseGridGeometry<GV, Traits>;
    using GridIndexType = typename IndexTraits<GV>::GridIndex;
    using LocalIndexType = typename IndexTraits<GV>::LocalIndex;
 
    static const int dim = GV::dimension;
    static const int dimWorld = GV::dimensionworld;
 
    using Element = typename GV::template Codim<0>::Entity;
    using Intersection = typename GV::Intersection;
    using CoordScalar = typename GV::ctype;
 
public:
    using DiscretizationMethod = DiscretizationMethods::Box;
    static constexpr DiscretizationMethod discMethod{};

    using LocalView = typename Traits::template LocalView<ThisType, false>;
    using SubControlVolume = typename Traits::SubControlVolume;
    using SubControlVolumeFace = typename Traits::SubControlVolumeFace;
    using Extrusion = Extrusion_t<Traits>;
    using DofMapper = typename Traits::VertexMapper;
    using FeCache = Dune::LagrangeLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
    using GridView = GV;

    template<class FacetGridView, class CodimOneGridAdapter>
    BoxFacetCouplingFVGridGeometry(const GridView& gridView,
                                   const FacetGridView& facetGridView,
                                   const CodimOneGridAdapter& codimOneGridAdapter,
                                   bool verbose = false)
    : ParentType(gridView)
    , facetMapper_(gridView, Dune::mcmgLayout(Dune::template Codim<1>()))
    , cache_(*this)
    {
        update_(facetGridView, codimOneGridAdapter, verbose);
    }

    const DofMapper& dofMapper() const
    { return this->vertexMapper(); }

    std::size_t numScv() const
    {  return numScv_; }

    std::size_t numScvf() const
    { return numScvf_; }

    std::size_t numBoundaryScvf() const
    { return numBoundaryScvf_; }

    std::size_t numDofs() const
    { return this->vertexMapper().size(); }

    template<class FacetGridView, class CodimOneGridAdapter>
    void update(const GridView& gridView,
                const FacetGridView& facetGridView,
                const CodimOneGridAdapter& codimOneGridAdapter,
                bool verbose = false)
    {
        ParentType::update(gridView);
        updateFacetMapper_();
        update_(facetGridView, codimOneGridAdapter, verbose);
    }

    template<class FacetGridView, class CodimOneGridAdapter>
    void update(GridView&& gridView,
                const FacetGridView& facetGridView,
                const CodimOneGridAdapter& codimOneGridAdapter,
                bool verbose = false)
    {
        ParentType::update(std::move(gridView));
        updateFacetMapper_();
        update_(facetGridView, codimOneGridAdapter, verbose);
    }

    const FeCache& feCache() const
    { return feCache_; }

    bool dofOnBoundary(unsigned int dofIdx) const
    { return boundaryDofIndices_[dofIdx]; }

    bool dofOnInteriorBoundary(unsigned int dofIdx) const
    { return interiorBoundaryDofIndices_[dofIdx]; }

    bool isOnInteriorBoundary(const Element& element, const Intersection& intersection) const
    { return facetIsOnInteriorBoundary_[ facetMapper_.subIndex(element, intersection.indexInInside(), 1) ]; }

    bool dofOnPeriodicBoundary(GridIndexType dofIdx) const
    { return false; }

    [[deprecated("Will be removed after release 3.11. Use periodicallyMappedDofs, returning a range of dofs")]]
    GridIndexType periodicallyMappedDof(GridIndexType dofIdx) const
    { DUNE_THROW(Dune::InvalidStateException, "Periodic boundaries are not supported by the facet coupling scheme"); }

    std::span<const GridIndexType> periodicallyMappedDofs(GridIndexType dofIdx) const
    { DUNE_THROW(Dune::InvalidStateException, "Periodic boundaries are not supported by the facet coupling scheme"); }

    friend inline LocalView localView(const BoxFacetCouplingFVGridGeometry& gg)
    { return { gg.cache_ }; }
 
private:
 
    class BoxFacetCouplingGridGeometryCache
    {
        friend class BoxFacetCouplingFVGridGeometry;
    public:
        using GeometryHelper = Detail::BoxFacetCouplingGeometryHelper_t<GV, Traits>;
        using DofHelper = typename GeometryHelper::DofHelper;
 
        explicit BoxFacetCouplingGridGeometryCache(const BoxFacetCouplingFVGridGeometry& gg)
        : gridGeometry_(&gg)
        {}
 
        const BoxFacetCouplingFVGridGeometry& gridGeometry() const
        { return *gridGeometry_; }
 
    private:
        const BoxFacetCouplingFVGridGeometry* gridGeometry_;
    };
 
public:
    using Cache = BoxFacetCouplingGridGeometryCache;
 
private:
 
    void updateFacetMapper_()
    {
        facetMapper_.update(this->gridView());
    }
 
    template<class FacetGridView, class CodimOneGridAdapter>
    void update_(const FacetGridView& facetGridView,
                 const CodimOneGridAdapter& codimOneGridAdapter,
                 bool verbose)
    {
        // enrich the vertex mapper subject to the provided facet grid
        this->vertexMapper().enrich(facetGridView, codimOneGridAdapter, verbose);
 
        // save global data on the grid's scvs and scvfs
        // TODO do we need those information?
        numScv_ = 0;
        numScvf_ = 0;
        numBoundaryScvf_ = 0;
 
        const auto numDof = numDofs();
        boundaryDofIndices_.assign(numDof, false);
        interiorBoundaryDofIndices_.assign(numDof, false);
        facetIsOnInteriorBoundary_.assign(this->gridView().size(1), false);
        for (const auto& element : elements(this->gridView()))
        {
            numScv_ += element.subEntities(dim);
            numScvf_ += element.subEntities(dim-1);
 
            const auto elementGeometry = element.geometry();
            const auto refElement = referenceElement(elementGeometry);
 
            // store the sub control volume face indices on the domain/interior boundary
            // skip handled facets (necessary for e.g. Dune::FoamGrid)
            std::vector<unsigned int> handledFacets;
            for (const auto& intersection : intersections(this->gridView(), element))
            {
                if (std::count(handledFacets.begin(), handledFacets.end(), intersection.indexInInside()))
                    continue;
 
                handledFacets.push_back(intersection.indexInInside());
 
                // determine if all corners live on the facet grid
                const auto isGeometry = intersection.geometry();
                const auto numFaceCorners = isGeometry.corners();
                const auto idxInInside = intersection.indexInInside();
                const auto boundary = intersection.boundary();
 
                std::vector<LocalIndexType> vIndicesLocal(numFaceCorners);
                for (int i = 0; i < numFaceCorners; ++i)
                    vIndicesLocal[i] = static_cast<LocalIndexType>(refElement.subEntity(idxInInside, 1, i, dim));
 
                std::vector<GridIndexType> gridVertexIndices(numFaceCorners);
                for (int i = 0; i < numFaceCorners; ++i)
                    gridVertexIndices[i] = this->vertexMapper().vertexIndex(element, vIndicesLocal[i], dim);
 
                // if all vertices are living on the facet grid, this is an interiour boundary
                const bool isOnFacet = codimOneGridAdapter.composeFacetElement(gridVertexIndices);
 
                // make sure there are no periodic boundaries
                if (boundary && intersection.neighbor())
                    DUNE_THROW(Dune::InvalidStateException, "Periodic boundaries are not supported by the box facet coupling scheme");
 
                if (isOnFacet || boundary)
                {
                    numScvf_ += numFaceCorners;
                    numBoundaryScvf_ += int(boundary)*numFaceCorners;
 
                    // Mark vertices to be on domain and/or interior boundary
                    for (int i = 0; i < numFaceCorners; ++i)
                    {
                        const auto dofIndex = this->vertexMapper().subIndex(element, vIndicesLocal[i], dim);
                        if (boundary) boundaryDofIndices_[ dofIndex ] = boundary && !isOnFacet;
                        if (isOnFacet)
                        {
                            interiorBoundaryDofIndices_[ dofIndex ] = true;
                            facetIsOnInteriorBoundary_[ facetMapper_.subIndex(element, idxInInside, 1) ] = true;
                        }
                    }
                }
            }
        }
    }
 
    const FeCache feCache_;
 
    // Information on the global number of geometries
    // TODO do we need those information?
    std::size_t numScv_;
    std::size_t numScvf_;
    std::size_t numBoundaryScvf_;
 
    // vertices on domain/interior boundaries
    std::vector<bool> boundaryDofIndices_;
    std::vector<bool> interiorBoundaryDofIndices_;
 
    // facet mapper and markers which facets lie on interior boundaries
    typename Traits::FacetMapper facetMapper_;
    std::vector<bool> facetIsOnInteriorBoundary_;
 
    Cache cache_;
};
 
} // end namespace Dumux
 
#endif