discretization/facecentered/diamond/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_DISCRETIZATION_FACECENTERED_DIAMOND_FV_GRID_GEOMETRY
#define DUMUX_DISCRETIZATION_FACECENTERED_DIAMOND_FV_GRID_GEOMETRY
 
#include <memory>
#include <ranges>
#include <unordered_map>
 
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/geometry/type.hh>
 
#include <dumux/common/defaultmappertraits.hh>
#include <dumux/common/indextraits.hh>
#include <dumux/common/math.hh>
#include <dumux/geometry/volume.hh>
#include <dumux/geometry/center.hh>
#include <dumux/discretization/basegridgeometry.hh>
#include <dumux/discretization/checkoverlapsize.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/extrusion.hh>
#include <dumux/discretization/nonconformingfecache.hh>
 
#include <dumux/discretization/facecentered/diamond/subcontrolvolume.hh>
#include <dumux/discretization/facecentered/diamond/subcontrolvolumeface.hh>
#include <dumux/discretization/facecentered/diamond/fvelementgeometry.hh>
#include <dumux/discretization/facecentered/diamond/geometryhelper.hh>
#include <dumux/discretization/boundaryface.hh>
 
#include <dumux/io/grid/periodicgridtraits.hh>
 
namespace Dumux {
 
namespace Detail {
template<class GV, class T>
using FaceCenteredDiamondGeometryHelper_t = Dune::Std::detected_or_t<
    Dumux::DiamondGeometryHelper<GV, typename T::SubControlVolume, typename T::SubControlVolumeFace>,
    SpecifiesGeometryHelper,
    T
>;
} // end namespace Detail

template<class GridView, class ScvRule = Dumux::QuadratureRules::MidpointQuadrature, class ScvfRule = Dumux::QuadratureRules::MidpointQuadrature>
using FaceCenteredDiamondQuadratureTraits = CVFE::DefaultQuadratureTraits<GridView, ScvRule, ScvfRule>;

template<class GridView, class QuadratureTraits = FaceCenteredDiamondQuadratureTraits<GridView>>
struct FaceCenteredDiamondDefaultGridGeometryTraits : public DefaultMapperTraits<GridView>, public QuadratureTraits
{
    using SubControlVolume = FaceCenteredDiamondSubControlVolume<GridView>;
    using SubControlVolumeFace = FaceCenteredDiamondSubControlVolumeFace<GridView>;
    using DofMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
 
    template<class GridGeometry, bool enableCache>
    using LocalView = FaceCenteredDiamondFVElementGeometry<GridGeometry, enableCache>;
};

template<class GV,
         bool enableCaching = true,
         class Traits = FaceCenteredDiamondDefaultGridGeometryTraits<GV>>
class FaceCenteredDiamondFVGridGeometry
: public BaseGridGeometry<GV, Traits>
{
    using ThisType = FaceCenteredDiamondFVGridGeometry<GV, enableCaching, Traits>;
    using ParentType = BaseGridGeometry<GV, Traits>;
    using GridIndexType = typename IndexTraits<GV>::GridIndex;
    using LocalIndexType = typename IndexTraits<GV>::SmallLocalIndex;
    using Element = typename GV::template Codim<0>::Entity;
 
    using Scalar = typename GV::ctype;
 
    static const int dim = GV::dimension;
    static const int dimWorld = GV::dimensionworld;
 
    static_assert(dim > 1, "Only implemented for dim > 1");
 
public:
    using DiscretizationMethod = DiscretizationMethods::FCDiamond;
    static constexpr DiscretizationMethod discMethod = DiscretizationMethod{};
    static constexpr bool cachingEnabled = true;
 
    static constexpr bool enableHybridCVFE = false;

    using LocalView = typename Traits::template LocalView<ThisType, true>;
    using SubControlVolume = typename Traits::SubControlVolume;
    using SubControlVolumeFace = typename Traits::SubControlVolumeFace;
    using GridView = GV;
    using DofMapper = typename Traits::DofMapper;
    using Extrusion = Extrusion_t<Traits>;
    using FeCache = NonconformingFECache<Scalar, Scalar, dim>;
    using SupportsPeriodicity = typename PeriodicGridTraits<typename GV::Grid>::SupportsPeriodicity;
    using ScvQuadratureRule = typename Traits::ScvQuadratureRule;
    using ScvfQuadratureRule = typename Traits::ScvfQuadratureRule;
    using BoundaryFace = Experimental::BoundaryFace<GV>;

    FaceCenteredDiamondFVGridGeometry(const GridView& gridView, const std::string& paramGroup = "")
    : ParentType(gridView)
    , dofMapper_(gridView, Dune::mcmgLayout(Dune::Codim<1>{}))
    , cache_(*this)
    , periodicGridTraits_(this->gridView().grid())
    {
        update_();
    }

    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->gridView().size(1); }

    void update(const GridView& gridView)
    {
        ParentType::update(gridView);
        update_();
    }

    void update(GridView&& gridView)
    {
        ParentType::update(std::move(gridView));
        update_();
    }

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

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

    const DofMapper& dofMapper() const
    { return dofMapper_; }

    bool dofOnPeriodicBoundary(GridIndexType dofIdx) const
    { return periodicFaceMap_.count(dofIdx); }

    [[deprecated("Will be removed after release 3.11. Use periodicallyMappedDofs, returning a range of dofs")]]
    GridIndexType periodicallyMappedDof(GridIndexType dofIdx) const
    { return periodicFaceMap_.at(dofIdx)[0]; }

    const std::ranges::range auto periodicallyMappedDofs(GridIndexType dofIdx) const
    { return std::views::all(periodicFaceMap_.at(dofIdx)); }

    const std::unordered_map<GridIndexType, std::array<GridIndexType, 1>>& periodicDofMap() const
    { return periodicFaceMap_; }

    friend inline LocalView localView(const FaceCenteredDiamondFVGridGeometry& gg)
    { return { gg.cache_ }; }
 
private:
 
    class FCDiamondGridGeometryCache
    {
        friend class FaceCenteredDiamondFVGridGeometry;
    public:
        using GeometryHelper = Detail::FaceCenteredDiamondGeometryHelper_t<GV, Traits>;
        using DofHelper = typename GeometryHelper::DofHelper;
 
        explicit FCDiamondGridGeometryCache(const FaceCenteredDiamondFVGridGeometry& gg)
        : gridGeometry_(&gg)
        {}
 
        const FaceCenteredDiamondFVGridGeometry& 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]; }

        bool hasBoundaryScvf(GridIndexType eIdx) const
        { return hasBoundaryScvf_[eIdx]; }

        const auto& boundaryFaces(GridIndexType eIdx) const
        { return boundaryFaces_[eIdx]; }

        const auto& boundaryFaceScvfRanges(GridIndexType eIdx) const
        { return boundaryFaceScvfRanges_[eIdx]; }
 
    private:
        void clear_()
        {
            scvs_.clear();
            scvfs_.clear();
            hasBoundaryScvf_.clear();
            boundaryFaces_.clear();
            boundaryFaceScvfRanges_.clear();
        }
 
        std::vector<std::vector<SubControlVolume>> scvs_;
        std::vector<std::vector<SubControlVolumeFace>> scvfs_;
        std::vector<bool> hasBoundaryScvf_;
        std::vector<Dune::ReservedVector<BoundaryFace, 2*dim>> boundaryFaces_;
        std::vector<Dune::ReservedVector<std::array<LocalIndexType, 2>, 2*dim>> boundaryFaceScvfRanges_;
 
        const FaceCenteredDiamondFVGridGeometry* gridGeometry_;
    };
 
public:
    using Cache = FCDiamondGridGeometryCache;
private:
    using GeometryHelper = typename Cache::GeometryHelper;

    void update_()
    {
        // clear containers (necessary after grid refinement)
        cache_.clear_();
        dofMapper_.update(this->gridView());
 
        // determine size of containers
        const auto numElements = this->gridView().size(0);
        cache_.scvs_.resize(numElements);
        cache_.scvfs_.resize(numElements);
        cache_.hasBoundaryScvf_.resize(numElements, false);
        cache_.boundaryFaces_.resize(numElements);
        cache_.boundaryFaceScvfRanges_.resize(numElements);
 
        boundaryDofIndices_.assign(numDofs(), false);
 
        numScv_ = 0;
        numScvf_ = 0;
        numBoundaryScvf_ = 0;
 
        // Build the scvs and scv faces
        for (const auto& element : elements(this->gridView()))
        {
            const auto eIdx = this->elementMapper().index(element);
 
            const auto geometry = element.geometry();
            GeometryHelper geometryHelper(geometry);
 
            // build the scvs
            cache_.scvs_[eIdx].reserve(geometryHelper.numScv());
            numScv_ += geometryHelper.numScv();
            for (LocalIndexType localScvIdx = 0; localScvIdx < geometryHelper.numScv(); ++localScvIdx)
            {
                const auto dofIndex = dofMapper().subIndex(element, localScvIdx, 1);
                const auto& corners = geometryHelper.getScvCorners(localScvIdx);
                const auto volume = Dumux::convexPolytopeVolume<dim>(
                    SubControlVolume::Traits::geometryType(geometry.type()),
                    [&](unsigned int i){ return corners[i]; }
                );
 
                cache_.scvs_[eIdx].emplace_back(
                    volume,
                    geometryHelper.facetCenter(localScvIdx),
                    Dumux::center(corners),
                    localScvIdx,
                    eIdx,
                    dofIndex
                );
            }
 
            // build interior scvfs
            LocalIndexType localScvfIdx = 0;
            cache_.scvfs_[eIdx].reserve(geometryHelper.numInteriorScvf());
            numScvf_ += geometryHelper.numInteriorScvf();
            for (; localScvfIdx < geometryHelper.numInteriorScvf(); ++localScvfIdx)
            {
                const auto& corners = geometryHelper.getScvfCorners(localScvfIdx);
                const auto& scvPair = geometryHelper.getInsideOutsideScvForScvf(localScvfIdx);
                const auto area = Dumux::convexPolytopeVolume<dim-1>(
                    SubControlVolumeFace::Traits::interiorGeometryType(geometry.type()),
                    [&](unsigned int i){ return corners[i]; }
                );
 
                cache_.scvfs_[eIdx].emplace_back(
                    Dumux::center(corners),
                    area,
                    geometryHelper.normal(corners, scvPair),
                    scvPair,
                    localScvfIdx
                );
            }
 
            // build boundary scvfs
            LocalIndexType numBoundaryFaces = 0;
            for (const auto& intersection : intersections(this->gridView(), element))
            {
                if (onDomainBoundary_(intersection))
                {
                    // store information that the face dof is on a boundary
                    const LocalIndexType localFacetIndex = intersection.indexInInside();
                    const auto dofIndex = dofMapper().subIndex(element, localFacetIndex, 1);
                    boundaryDofIndices_[dofIndex] = true;
 
                    // and that the element has a boundary face
                    cache_.hasBoundaryScvf_[eIdx] = true;
 
                    // add boundary face
                    {
                        const auto geo = intersection.geometry();
                        cache_.boundaryFaces_[eIdx].push_back(BoundaryFace{
                            geo.center(),
                            geo.volume(),
                            intersection.centerUnitOuterNormal(),
                            numBoundaryFaces++,
                            static_cast<LocalIndexType>(localFacetIndex),
                            typename BoundaryFace::Traits::BoundaryFlag{intersection}
                        });
 
                        // record the scvf subrange for this boundary face: {offset, count=1}
                        cache_.boundaryFaceScvfRanges_[eIdx].push_back(std::array<LocalIndexType, 2>{{
                            localScvfIdx, 1
                        }});
                    }
 
                    // add boundary scvf
                    const auto geo = intersection.geometry();
                    cache_.scvfs_[eIdx].emplace_back(
                        geo.center(),
                        geo.volume(),
                        intersection.centerUnitOuterNormal(),
                        std::array<LocalIndexType, 2>{{localFacetIndex, localFacetIndex}},
                        localScvfIdx,
                        typename SubControlVolumeFace::Traits::BoundaryFlag{ intersection }
                    );
 
                    // increment local and global counters
                    ++localScvfIdx;
                    ++numBoundaryScvf_;
                    ++numScvf_;
                }
 
                // handle periodic boundaries
                if (onPeriodicBoundary_(intersection))
                {
                    const LocalIndexType localFacetIndex = intersection.indexInInside();
                    const auto dofIndex = dofMapper().subIndex(element, localFacetIndex, 1);
 
                    this->setPeriodic();
 
                    const auto& otherElement = intersection.outside();
 
                    LocalIndexType otherIntersectionLocalIdx = 0;
                    bool periodicFaceFound = false;
 
                    for (const auto& otherIntersection : intersections(this->gridView(), otherElement))
                    {
                        if (periodicFaceFound)
                            continue;
 
                        if (Dune::FloatCmp::eq(intersection.centerUnitOuterNormal()*otherIntersection.centerUnitOuterNormal(), -1.0, 1e-7))
                        {
                            const auto periodicDofIdx = dofMapper().subIndex(otherElement, otherIntersectionLocalIdx, 1);
                            periodicFaceMap_[dofIndex][0] = periodicDofIdx;
                            periodicFaceFound = true;
                        }
 
                        ++otherIntersectionLocalIdx;
                    }
                }
            }
        }
    }
 
    bool onDomainBoundary_(const typename GridView::Intersection& intersection) const
    {
        return !intersection.neighbor() && intersection.boundary();
    }
 
    bool onProcessorBoundary_(const typename GridView::Intersection& intersection) const
    {
        return !intersection.neighbor() && !intersection.boundary();
    }
 
    bool onPeriodicBoundary_(const typename GridView::Intersection& intersection) const
    {
        return periodicGridTraits_.isPeriodic(intersection);
    }
 
    // faces on the boundary
    std::vector<bool> boundaryDofIndices_;
 
    DofMapper dofMapper_;
 
    std::size_t numScv_;
    std::size_t numScvf_;
    std::size_t numBoundaryScvf_;
 
    // a map for periodic boundary vertices
    std::unordered_map<GridIndexType, std::array<GridIndexType, 1>> periodicFaceMap_;
 
    const FeCache feCache_;
 
    Cache cache_;
 
    PeriodicGridTraits<typename GridView::Grid> periodicGridTraits_;
};
 
} // end namespace Dumux
 
#endif