discretization/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:
/*****************************************************************************
 *   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 3 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/>.   *
 *****************************************************************************/
#ifndef DUMUX_DISCRETIZATION_BOX_GRID_FVGEOMETRY_HH
#define DUMUX_DISCRETIZATION_BOX_GRID_FVGEOMETRY_HH
 
#include <utility>
#include <unordered_map>
#include <array>
#include <vector>
 
#include <dune/localfunctions/lagrange/lagrangelfecache.hh>
 
#include <dumux/discretization/method.hh>
#include <dumux/common/indextraits.hh>
#include <dumux/common/defaultmappertraits.hh>
#include <dumux/discretization/basegridgeometry.hh>
#include <dumux/discretization/box/boxgeometryhelper.hh>
#include <dumux/discretization/box/fvelementgeometry.hh>
#include <dumux/discretization/box/subcontrolvolume.hh>
#include <dumux/discretization/box/subcontrolvolumeface.hh>
#include <dumux/discretization/extrusion.hh>
 
namespace Dumux {

template<class GridView, class MapperTraits = DefaultMapperTraits<GridView>>
struct BoxDefaultGridGeometryTraits
: public MapperTraits
{
    using SubControlVolume = BoxSubControlVolume<GridView>;
    using SubControlVolumeFace = BoxSubControlVolumeFace<GridView>;
 
    template<class GridGeometry, bool enableCache>
    using LocalView = BoxFVElementGeometry<GridGeometry, enableCache>;
};

template<class Scalar,
         class GridView,
         bool enableGridGeometryCache = false,
         class Traits = BoxDefaultGridGeometryTraits<GridView> >
class BoxFVGridGeometry;

template<class Scalar, class GV, class Traits>
class BoxFVGridGeometry<Scalar, GV, true, Traits>
: public BaseGridGeometry<GV, Traits>
{
    using ThisType = BoxFVGridGeometry<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;
 
    using GeometryHelper = BoxGeometryHelper<GV, dim,
                                             typename Traits::SubControlVolume,
                                             typename Traits::SubControlVolumeFace>;
 
    // Remove this after release 3.6 when the deprecated local view constructor is removed
    using LV_ = typename Traits::template LocalView<ThisType, true>;
    friend LV_;
 
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;

    BoxFVGridGeometry(const GridView gridView)
    : ParentType(gridView)
    , cache_(*this)
    {
        update_();
    }

    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(); }
 

    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_; }

    [[deprecated("Will be removed after release 3.6")]]
    const std::vector<SubControlVolume>& scvs(GridIndexType eIdx) const
    { return cache_.scvs_[eIdx]; }

    [[deprecated("Will be removed after release 3.6")]]
    const std::vector<SubControlVolumeFace>& scvfs(GridIndexType eIdx) const
    { return cache_.scvfs_[eIdx]; }

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

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

    GridIndexType periodicallyMappedDof(GridIndexType dofIdx) const
    { return periodicVertexMap_.at(dofIdx); }

    const std::unordered_map<GridIndexType, GridIndexType>& periodicVertexMap() const
    { return periodicVertexMap_; }

    [[deprecated("Will be removed after release 3.6")]]
    bool hasBoundaryScvf(GridIndexType eIdx) const
    { return cache_.hasBoundaryScvf_[eIdx]; }

    friend inline LocalView localView(const BoxFVGridGeometry& gg)
    { return { gg.cache_ }; }
 
private:
 
    class BoxGridGeometryCache
    {
        friend class BoxFVGridGeometry;
    public:
        explicit BoxGridGeometryCache(const BoxFVGridGeometry& gg)
        : gridGeometry_(&gg)
        {}
 
        const BoxFVGridGeometry& 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 std::vector<std::array<LocalIndexType, 2>>& scvfBoundaryGeometryKeys(GridIndexType eIdx) const
        { return scvfBoundaryGeometryKeys_.at(eIdx); }
 
    private:
        void clear_()
        {
            scvs_.clear();
            scvfs_.clear();
            hasBoundaryScvf_.clear();
            scvfBoundaryGeometryKeys_.clear();
        }
 
        std::vector<std::vector<SubControlVolume>> scvs_;
        std::vector<std::vector<SubControlVolumeFace>> scvfs_;
        std::vector<bool> hasBoundaryScvf_;
        std::unordered_map<GridIndexType, std::vector<std::array<LocalIndexType, 2>>> scvfBoundaryGeometryKeys_;
 
        const BoxFVGridGeometry* gridGeometry_;
    };
 
public:
    using Cache = BoxGridGeometryCache;
 
private:
    void update_()
    {
        cache_.clear_();
 
        const auto numElements = this->gridView().size(0);
        cache_.scvs_.resize(numElements);
        cache_.scvfs_.resize(numElements);
        cache_.hasBoundaryScvf_.resize(numElements, false);
 
        boundaryDofIndices_.assign(numDofs(), false);
 
        numScv_ = 0;
        numScvf_ = 0;
        numBoundaryScvf_ = 0;
        // Build the SCV and SCV faces
        for (const auto& element : elements(this->gridView()))
        {
            // fill the element map with seeds
            const auto eIdx = this->elementMapper().index(element);
 
            // count
            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)
            {
                const auto dofIdxGlobal = this->vertexMapper().subIndex(element, scvLocalIdx, dim);
 
                cache_.scvs_[eIdx][scvLocalIdx] = SubControlVolume(
                    geometryHelper,
                    scvLocalIdx,
                    eIdx,
                    dofIdxGlobal
                );
            }
 
            // 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::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
                                                             static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
 
                cache_.scvfs_[eIdx][scvfLocalIdx] = SubControlVolumeFace(
                    geometryHelper,
                    element,
                    elementGeometry,
                    scvfLocalIdx,
                    std::move(localScvIndices),
                    false
                );
            }
 
            // construct the sub control volume faces on the domain boundary
            for (const auto& intersection : intersections(this->gridView(), element))
            {
                if (intersection.boundary() && !intersection.neighbor())
                {
                    const auto isGeometry = intersection.geometry();
                    cache_.hasBoundaryScvf_[eIdx] = true;
 
                    // count
                    numScvf_ += isGeometry.corners();
                    numBoundaryScvf_ += isGeometry.corners();
 
                    for (unsigned int isScvfLocalIdx = 0; isScvfLocalIdx < isGeometry.corners(); ++isScvfLocalIdx)
                    {
                        // find the scvs this scvf is belonging to
                        const LocalIndexType insideScvIdx = static_cast<LocalIndexType>(refElement.subEntity(intersection.indexInInside(), 1, isScvfLocalIdx, dim));
                        std::vector<LocalIndexType> localScvIndices = {insideScvIdx, insideScvIdx};
 
                        cache_.scvfs_[eIdx].emplace_back(
                            geometryHelper,
                            intersection,
                            isGeometry,
                            isScvfLocalIdx,
                            scvfLocalIdx,
                            std::move(localScvIndices),
                            true
                        );
 
                        cache_.scvfBoundaryGeometryKeys_[eIdx].emplace_back(std::array<LocalIndexType, 2>{{
                            static_cast<LocalIndexType>(intersection.indexInInside()),
                            static_cast<LocalIndexType>(isScvfLocalIdx)
                        }});
 
                        // increment local counter
                        scvfLocalIdx++;
                    }
 
                    // add all vertices on the intersection to the set of
                    // boundary vertices
                    const auto fIdx = intersection.indexInInside();
                    const auto numFaceVerts = refElement.size(fIdx, 1, dim);
                    for (int localVIdx = 0; localVIdx < numFaceVerts; ++localVIdx)
                    {
                        const auto vIdx = refElement.subEntity(fIdx, 1, localVIdx, dim);
                        const auto vIdxGlobal = this->vertexMapper().subIndex(element, vIdx, dim);
                        boundaryDofIndices_[vIdxGlobal] = true;
                    }
                }
 
                // inform the grid geometry if we have periodic boundaries
                else if (intersection.boundary() && intersection.neighbor())
                {
                    this->setPeriodic();
 
                    // find the mapped periodic vertex of all vertices on periodic boundaries
                    const auto fIdx = intersection.indexInInside();
                    const auto numFaceVerts = refElement.size(fIdx, 1, dim);
                    const auto eps = 1e-7*(elementGeometry.corner(1) - elementGeometry.corner(0)).two_norm();
                    for (int localVIdx = 0; localVIdx < numFaceVerts; ++localVIdx)
                    {
                        const auto vIdx = refElement.subEntity(fIdx, 1, localVIdx, dim);
                        const auto vIdxGlobal = this->vertexMapper().subIndex(element, vIdx, dim);
                        const auto vPos = elementGeometry.corner(vIdx);
 
                        const auto& outside = intersection.outside();
                        const auto outsideGeometry = outside.geometry();
                        for (const auto& isOutside : intersections(this->gridView(), outside))
                        {
                            // only check periodic vertices of the periodic neighbor
                            if (isOutside.boundary() && isOutside.neighbor())
                            {
                                const auto fIdxOutside = isOutside.indexInInside();
                                const auto numFaceVertsOutside = refElement.size(fIdxOutside, 1, dim);
                                for (int localVIdxOutside = 0; localVIdxOutside < numFaceVertsOutside; ++localVIdxOutside)
                                {
                                    const auto vIdxOutside = refElement.subEntity(fIdxOutside, 1, localVIdxOutside, dim);
                                    const auto vPosOutside = outsideGeometry.corner(vIdxOutside);
                                    const auto shift = std::abs((this->bBoxMax()-this->bBoxMin())*intersection.centerUnitOuterNormal());
                                    if (std::abs((vPosOutside-vPos).two_norm() - shift) < eps)
                                        periodicVertexMap_[vIdxGlobal] = this->vertexMapper().subIndex(outside, vIdxOutside, dim);
                                }
                            }
                        }
                    }
                }
            }
        }
 
        // error check: periodic boundaries currently don't work for box in parallel
        if (this->isPeriodic() && this->gridView().comm().size() > 1)
            DUNE_THROW(Dune::NotImplemented, "Periodic boundaries for box method for parallel simulations!");
    }
 
    const FeCache feCache_;
 
    std::size_t numScv_;
    std::size_t numScvf_;
    std::size_t numBoundaryScvf_;
 
    // vertices on the boundary
    std::vector<bool> boundaryDofIndices_;
 
    // a map for periodic boundary vertices
    std::unordered_map<GridIndexType, GridIndexType> periodicVertexMap_;
 
    Cache cache_;
};

template<class Scalar, class GV, class Traits>
class BoxFVGridGeometry<Scalar, GV, false, Traits>
: public BaseGridGeometry<GV, Traits>
{
    using ThisType = BoxFVGridGeometry<Scalar, GV, false, Traits>;
    using ParentType = BaseGridGeometry<GV, Traits>;
    using GridIndexType = typename IndexTraits<GV>::GridIndex;
 
    static const int dim = GV::dimension;
    static const int dimWorld = GV::dimensionworld;
 
    using Element = typename GV::template Codim<0>::Entity;
    using CoordScalar = typename GV::ctype;
 
    // Remove this after release 3.6 when the deprecated local view constructor is removed
    using LV_ = typename Traits::template LocalView<ThisType, false>;
    friend LV_;
 
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;

    BoxFVGridGeometry(const GridView gridView)
    : ParentType(gridView)
    , cache_(*this)
    {
        update_();
    }

    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(); }
 

    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]; }

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

    GridIndexType periodicallyMappedDof(GridIndexType dofIdx) const
    { return periodicVertexMap_.at(dofIdx); }

    const std::unordered_map<GridIndexType, GridIndexType>& periodicVertexMap() const
    { return periodicVertexMap_; }

    friend inline LocalView localView(const BoxFVGridGeometry& gg)
    { return { gg.cache_ }; }
 
private:
 
    class BoxGridGeometryCache
    {
        friend class BoxFVGridGeometry;
    public:
        explicit BoxGridGeometryCache(const BoxFVGridGeometry& gg)
        : gridGeometry_(&gg)
        {}
 
        const BoxFVGridGeometry& gridGeometry() const
        { return *gridGeometry_; }
 
    private:
        const BoxFVGridGeometry* gridGeometry_;
    };
 
public:
    using Cache = BoxGridGeometryCache;
 
private:
 
    void update_()
    {
        boundaryDofIndices_.assign(numDofs(), false);
 
        // save global data on the grid's scvs and scvfs
        // TODO do we need those information?
        numScv_ = 0;
        numScvf_ = 0;
        numBoundaryScvf_ = 0;
        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 boundary
            for (const auto& intersection : intersections(this->gridView(), element))
            {
                if (intersection.boundary() && !intersection.neighbor())
                {
                    const auto isGeometry = intersection.geometry();
                    numScvf_ += isGeometry.corners();
                    numBoundaryScvf_ += isGeometry.corners();
 
                    // add all vertices on the intersection to the set of
                    // boundary vertices
                    const auto fIdx = intersection.indexInInside();
                    const auto numFaceVerts = refElement.size(fIdx, 1, dim);
                    for (int localVIdx = 0; localVIdx < numFaceVerts; ++localVIdx)
                    {
                        const auto vIdx = refElement.subEntity(fIdx, 1, localVIdx, dim);
                        const auto vIdxGlobal = this->vertexMapper().subIndex(element, vIdx, dim);
                        boundaryDofIndices_[vIdxGlobal] = true;
                    }
                }
 
                // inform the grid geometry if we have periodic boundaries
                else if (intersection.boundary() && intersection.neighbor())
                {
                    this->setPeriodic();
 
                    // find the mapped periodic vertex of all vertices on periodic boundaries
                    const auto fIdx = intersection.indexInInside();
                    const auto numFaceVerts = refElement.size(fIdx, 1, dim);
                    const auto eps = 1e-7*(elementGeometry.corner(1) - elementGeometry.corner(0)).two_norm();
                    for (int localVIdx = 0; localVIdx < numFaceVerts; ++localVIdx)
                    {
                        const auto vIdx = refElement.subEntity(fIdx, 1, localVIdx, dim);
                        const auto vIdxGlobal = this->vertexMapper().subIndex(element, vIdx, dim);
                        const auto vPos = elementGeometry.corner(vIdx);
 
                        const auto& outside = intersection.outside();
                        const auto outsideGeometry = outside.geometry();
                        for (const auto& isOutside : intersections(this->gridView(), outside))
                        {
                            // only check periodic vertices of the periodic neighbor
                            if (isOutside.boundary() && isOutside.neighbor())
                            {
                                const auto fIdxOutside = isOutside.indexInInside();
                                const auto numFaceVertsOutside = refElement.size(fIdxOutside, 1, dim);
                                for (int localVIdxOutside = 0; localVIdxOutside < numFaceVertsOutside; ++localVIdxOutside)
                                {
                                    const auto vIdxOutside = refElement.subEntity(fIdxOutside, 1, localVIdxOutside, dim);
                                    const auto vPosOutside = outsideGeometry.corner(vIdxOutside);
                                    const auto shift = std::abs((this->bBoxMax()-this->bBoxMin())*intersection.centerUnitOuterNormal());
                                    if (std::abs((vPosOutside-vPos).two_norm() - shift) < eps)
                                        periodicVertexMap_[vIdxGlobal] = this->vertexMapper().subIndex(outside, vIdxOutside, dim);
                                }
                            }
                        }
                    }
                }
            }
        }
 
        // error check: periodic boundaries currently don't work for box in parallel
        if (this->isPeriodic() && this->gridView().comm().size() > 1)
            DUNE_THROW(Dune::NotImplemented, "Periodic boundaries for box method for parallel simulations!");
    }
 
    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 the boundary
    std::vector<bool> boundaryDofIndices_;
 
    // a map for periodic boundary vertices
    std::unordered_map<GridIndexType, GridIndexType> periodicVertexMap_;
 
    Cache cache_;
};
 
} // end namespace Dumux
 
#endif