discretization/staggered/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_STAGGERED_FV_GRID_GEOMETRY
#define DUMUX_DISCRETIZATION_STAGGERED_FV_GRID_GEOMETRY
 
#include <utility>
 
#include <dumux/common/indextraits.hh>
#include <dumux/discretization/basegridgeometry.hh>
#include <dumux/discretization/checkoverlapsize.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/extrusion.hh>
 
namespace Dumux {

template<class ActualGridGeometry>
class GridGeometryView
{
public:
 
    explicit GridGeometryView(const ActualGridGeometry* actualGridGeometry)
    : gridGeometry_(actualGridGeometry) {}

    using GridView = typename ActualGridGeometry::GridView;

    using DiscretizationMethod = DiscretizationMethods::Staggered;
    static constexpr DiscretizationMethod discMethod{};
 
    using LocalView = typename ActualGridGeometry::LocalView;

    static constexpr bool isCellCenter() { return false; }

    static constexpr bool isFace() {return false; }

    static constexpr auto cellCenterIdx()
    { return typename ActualGridGeometry::DofTypeIndices::CellCenterIdx{}; }

    static constexpr auto faceIdx()
    { return typename ActualGridGeometry::DofTypeIndices::FaceIdx{}; }

    const auto& gridView() const
    { return gridGeometry_->gridView(); }

    const auto& connectivityMap() const // TODO return correct map
    { return gridGeometry_->connectivityMap(); }

    const auto& vertexMapper() const
    { return gridGeometry_->vertexMapper(); }

    const auto& elementMapper() const
    { return gridGeometry_->elementMapper(); }

    const ActualGridGeometry& actualGridGeometry() const
    { return *gridGeometry_; }
 
protected:
    const ActualGridGeometry* gridGeometry_;
 
};

template <class ActualGridGeometry>
class CellCenterFVGridGeometry : public GridGeometryView<ActualGridGeometry>
{
    using ParentType = GridGeometryView<ActualGridGeometry>;
public:
 
    using ParentType::ParentType;

    static constexpr bool isCellCenter() { return true; }

    std::size_t numDofs() const
    { return this->gridGeometry_->numCellCenterDofs(); }
};

template <class ActualGridGeometry>
class FaceFVGridGeometry : public GridGeometryView<ActualGridGeometry>
{
    using ParentType = GridGeometryView<ActualGridGeometry>;
public:
 
    using ParentType::ParentType;

    static constexpr bool isFace() {return true; }

    std::size_t numDofs() const
    { return this->gridGeometry_->numFaceDofs(); }
};

template<class GridView,
         bool cachingEnabled,
         class Traits>
class StaggeredFVGridGeometry;

template<class GV, class T>
class StaggeredFVGridGeometry<GV, true, T>
: public BaseGridGeometry<GV, T>
{
    using ThisType = StaggeredFVGridGeometry<GV, true, T>;
    using ParentType = BaseGridGeometry<GV, T>;
    using GridIndexType = typename IndexTraits<GV>::GridIndex;
    using LocalIndexType = typename IndexTraits<GV>::LocalIndex;
    using Element = typename GV::template Codim<0>::Entity;
 
    using IntersectionMapper = typename T::IntersectionMapper;
    using GeometryHelper = typename T::GeometryHelper;
    using ConnectivityMap = typename T::template ConnectivityMap<ThisType>;
 
public:
    using Traits = typename T::PublicTraits;

    using DiscretizationMethod = DiscretizationMethods::Staggered;
    static constexpr DiscretizationMethod discMethod{};
 
    static constexpr int upwindSchemeOrder = T::upwindSchemeOrder;
    static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
    static constexpr bool cachingEnabled = true;

    using LocalView = typename T::template LocalView<ThisType, true>;
    using SubControlVolume = typename T::SubControlVolume;
    using SubControlVolumeFace = typename T::SubControlVolumeFace;
    using Extrusion = Extrusion_t<T>;
    using GridView = GV;
    using DofTypeIndices = typename T::DofTypeIndices;

    static constexpr auto cellCenterIdx()
    { return typename DofTypeIndices::CellCenterIdx{}; }

    static constexpr auto faceIdx()
    { return typename DofTypeIndices::FaceIdx{}; }

    static constexpr int upwindStencilOrder()
    {   return upwindSchemeOrder; }
 
    using CellCenterFVGridGeometryType = CellCenterFVGridGeometry<ThisType>;
    using FaceFVGridGeometryType = FaceFVGridGeometry<ThisType>;
 
    using FVGridGeometryTuple = std::tuple< CellCenterFVGridGeometry<ThisType>, FaceFVGridGeometry<ThisType> >;

    StaggeredFVGridGeometry(const GridView& gridView, const std::string& paramGroup = "")
    : ParentType(gridView)
    , intersectionMapper_(gridView)
    {
        // Check if the overlap size is what we expect
        if (!CheckOverlapSize<DiscretizationMethod>::isValid(gridView))
            DUNE_THROW(Dune::InvalidStateException, "The staggered discretization method needs at least an overlap of 1 for parallel computations. "
                                                     << " Set the parameter \"Grid.Overlap\" in the input file.");
 
        update_();
    }

    std::size_t numScv() const
    {
        return scvs_.size();
    }

    std::size_t numScvf() const
    {
        return scvfs_.size();
    }

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

    std::size_t numIntersections() const
    {
        return intersectionMapper_.numIntersections();
    }

    std::size_t numDofs() const
    { return numCellCenterDofs() + numFaceDofs(); }
 
    std::size_t numCellCenterDofs() const
    { return this->gridView().size(0); }
 
    std::size_t numFaceDofs() const
    { return this->gridView().size(1); }

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

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

    const SubControlVolume& scv(GridIndexType scvIdx) const
    {
        return scvs_[scvIdx];
    }

    const SubControlVolumeFace& scvf(GridIndexType scvfIdx) const
    {
        return scvfs_[scvfIdx];
    }

    const std::vector<GridIndexType>& scvfIndicesOfScv(GridIndexType scvIdx) const
    {
        return scvfIndicesOfScv_[scvIdx];
    }
 
    GridIndexType localToGlobalScvfIndex(GridIndexType eIdx, LocalIndexType localScvfIdx) const
    {
        return localToGlobalScvfIndices_[eIdx][localScvfIdx];
    }
 
    const SubControlVolumeFace& scvf(GridIndexType eIdx, LocalIndexType localScvfIdx) const
    {
        return scvf(localToGlobalScvfIndex(eIdx, localScvfIdx));
    }

    const ConnectivityMap &connectivityMap() const
    { return connectivityMap_; }

    std::unique_ptr<CellCenterFVGridGeometry<ThisType>> cellCenterFVGridGeometryPtr() const
    {
        return std::make_unique<CellCenterFVGridGeometry<ThisType>>(this);
    }

    std::unique_ptr<FaceFVGridGeometry<ThisType>> faceFVGridGeometryPtr() const
    {
        return std::make_unique<FaceFVGridGeometry<ThisType>>(this);
    }

    CellCenterFVGridGeometry<ThisType> cellCenterFVGridGeometry() const
    {
        return CellCenterFVGridGeometry<ThisType>(this);
    }

    FaceFVGridGeometry<ThisType> faceFVGridGeometry() const
    {
        return FaceFVGridGeometry<ThisType>(this);
    }

    bool hasBoundaryScvf(GridIndexType eIdx) const
    { return hasBoundaryScvf_[eIdx]; }
 
private:
 
    void updateIntersectionMapper_()
    {
        intersectionMapper_.update(this->gridView());
    }
 
    void update_()
    {
        // clear containers (necessary after grid refinement)
        scvs_.clear();
        scvfs_.clear();
        scvfIndicesOfScv_.clear();
        localToGlobalScvfIndices_.clear();
 
        // determine size of containers
        std::size_t numScvs = this->gridView().size(0);
        std::size_t numScvf = 0;
        for (const auto& element : elements(this->gridView()))
            numScvf += element.subEntities(1);
 
        // reserve memory
        scvs_.resize(numScvs);
        scvfs_.reserve(numScvf);
        scvfIndicesOfScv_.resize(numScvs);
        localToGlobalScvfIndices_.resize(numScvs);
        hasBoundaryScvf_.assign(numScvs, false);
 
        // Build the scvs and scv faces
        GridIndexType scvfIdx = 0;
        numBoundaryScvf_ = 0;
        for (const auto& element : elements(this->gridView()))
        {
            auto eIdx = this->elementMapper().index(element);
 
            // reserve memory for the localToGlobalScvfIdx map
            auto numLocalFaces = intersectionMapper_.numFaces(element);
            localToGlobalScvfIndices_[eIdx].resize(numLocalFaces);
 
            scvs_[eIdx] = SubControlVolume(element.geometry(), eIdx);
 
            // the element-wise index sets for finite volume geometry
            std::vector<GridIndexType> scvfsIndexSet;
            scvfsIndexSet.reserve(numLocalFaces);
 
            GeometryHelper geometryHelper(element, this->gridView());
 
            for (const auto& intersection : intersections(this->gridView(), element))
            {
                geometryHelper.updateLocalFace(intersectionMapper_, intersection);
                const int localFaceIndex = geometryHelper.localFaceIndex();
 
                // inner sub control volume faces
                if (intersection.neighbor())
                {
                    auto nIdx = this->elementMapper().index(intersection.outside());
                    scvfs_.emplace_back(intersection,
                                        intersection.geometry(),
                                        scvfIdx,
                                        std::vector<GridIndexType>({eIdx, nIdx}),
                                        geometryHelper);
                    localToGlobalScvfIndices_[eIdx][localFaceIndex] = scvfIdx;
                    scvfsIndexSet.push_back(scvfIdx++);
                }
                // boundary sub control volume faces
                else if (intersection.boundary())
                {
                    scvfs_.emplace_back(intersection,
                                        intersection.geometry(),
                                        scvfIdx,
                                        std::vector<GridIndexType>({eIdx, this->gridView().size(0) + numBoundaryScvf_++}),
                                        geometryHelper);
                    localToGlobalScvfIndices_[eIdx][localFaceIndex] = scvfIdx;
                    scvfsIndexSet.push_back(scvfIdx++);
 
                    hasBoundaryScvf_[eIdx] = true;
                }
            }
 
            // Save the scvf indices belonging to this scv to build up fv element geometries fast
            scvfIndicesOfScv_[eIdx] = scvfsIndexSet;
        }
 
        // build the connectivity map for an efficient assembly
        connectivityMap_.update(*this);
    }
 
    // mappers
    ConnectivityMap connectivityMap_;
    IntersectionMapper intersectionMapper_;
 
    std::vector<SubControlVolume> scvs_;
    std::vector<SubControlVolumeFace> scvfs_;
    std::vector<std::vector<GridIndexType>> scvfIndicesOfScv_;
    std::vector<std::vector<GridIndexType>> localToGlobalScvfIndices_;
    GridIndexType numBoundaryScvf_;
    std::vector<bool> hasBoundaryScvf_;
};

template<class GV, class T>
class StaggeredFVGridGeometry<GV, false, T>
: public BaseGridGeometry<GV, T>
{
    using ThisType = StaggeredFVGridGeometry<GV, false, T>;
    using ParentType = BaseGridGeometry<GV, T>;
    using GridIndexType = typename IndexTraits<GV>::GridIndex;
    using LocalIndexType = typename IndexTraits<GV>::LocalIndex;
    using Element = typename GV::template Codim<0>::Entity;
 
    using IntersectionMapper = typename T::IntersectionMapper;
    using ConnectivityMap = typename T::template ConnectivityMap<ThisType>;
 
public:
    using Traits = typename T::PublicTraits;

    using DiscretizationMethod = DiscretizationMethods::Staggered;
    static constexpr DiscretizationMethod discMethod{};
 
    static constexpr int upwindSchemeOrder = T::upwindSchemeOrder;
    static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
    static constexpr bool cachingEnabled = false;
 
    using GeometryHelper = typename T::GeometryHelper;

    using LocalView = typename T::template LocalView<ThisType, false>;
    using SubControlVolume = typename T::SubControlVolume;
    using SubControlVolumeFace = typename T::SubControlVolumeFace;
    using Extrusion = Extrusion_t<T>;
    using GridView = GV;
    using DofTypeIndices = typename T::DofTypeIndices;

    static constexpr auto cellCenterIdx()
    { return typename DofTypeIndices::CellCenterIdx{}; }

    static constexpr auto faceIdx()
    { return typename DofTypeIndices::FaceIdx{}; }

    static constexpr int upwindStencilOrder()
    {   return upwindSchemeOrder; }
 
    using CellCenterFVGridGeometryType = CellCenterFVGridGeometry<ThisType>;
    using FaceFVGridGeometryType = FaceFVGridGeometry<ThisType>;
 
    using FVGridGeometryTuple = std::tuple< CellCenterFVGridGeometry<ThisType>, FaceFVGridGeometry<ThisType> >;

    StaggeredFVGridGeometry(const GridView& gridView, const std::string& paramGroup = "")
    : ParentType(gridView)
    , intersectionMapper_(gridView)
    {
        // Check if the overlap size is what we expect
        if (!CheckOverlapSize<DiscretizationMethod>::isValid(gridView))
            DUNE_THROW(Dune::InvalidStateException, "The staggered discretization method needs at least an overlap of 1 for parallel computations. "
                                                     << " Set the parameter \"Grid.Overlap\" in the input file.");
 
        update_();
    }
 

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

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

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

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

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

    std::size_t numIntersections() const
    {
        return intersectionMapper_.numIntersections();
    }

    std::size_t numDofs() const
    { return numCellCenterDofs() + numFaceDofs(); }
 
    std::size_t numCellCenterDofs() const
    { return this->gridView().size(0); }
 
    std::size_t numFaceDofs() const
    { return this->gridView().size(1); }
 
    const std::vector<GridIndexType>& scvfIndicesOfScv(GridIndexType scvIdx) const
    { return scvfIndicesOfScv_[scvIdx]; }
 
    GridIndexType localToGlobalScvfIndex(GridIndexType eIdx, LocalIndexType localScvfIdx) const
    {
        return localToGlobalScvfIndices_[eIdx][localScvfIdx];
    }

    const ConnectivityMap &connectivityMap() const
    { return connectivityMap_; }

    std::unique_ptr<CellCenterFVGridGeometry<ThisType>> cellCenterFVGridGeometryPtr() const
    {
        return std::make_unique<CellCenterFVGridGeometry<ThisType>>(this);
    }

    std::unique_ptr<FaceFVGridGeometry<ThisType>> faceFVGridGeometryPtr() const
    {
        return std::make_unique<FaceFVGridGeometry<ThisType>>(this);
    }

    CellCenterFVGridGeometry<ThisType> cellCenterFVGridGeometry() const
    {
        return CellCenterFVGridGeometry<ThisType>(this);
    }

    FaceFVGridGeometry<ThisType> faceFVGridGeometry() const
    {
        return FaceFVGridGeometry<ThisType>(this);
    }

    const IntersectionMapper& intersectionMapper() const
    {
        return intersectionMapper_;
    }

    const std::vector<GridIndexType>& neighborVolVarIndices(GridIndexType scvIdx) const
    { return neighborVolVarIndices_[scvIdx]; }
 
private:
 
    void updateIntersectionMapper_()
    {
        intersectionMapper_.update(this->gridView());
    }
 
    void update_()
    {
        // clear containers (necessary after grid refinement)
        scvfIndicesOfScv_.clear();
        neighborVolVarIndices_.clear();
        localToGlobalScvfIndices_.clear();
 
        numScvs_ = numCellCenterDofs();
        numScvf_ = 0;
        numBoundaryScvf_ = 0;
        scvfIndicesOfScv_.resize(numScvs_);
        localToGlobalScvfIndices_.resize(numScvs_);
        neighborVolVarIndices_.resize(numScvs_);
 
        // Build the scvs and scv faces
        for (const auto& element : elements(this->gridView()))
        {
            auto eIdx = this->elementMapper().index(element);
 
            // the element-wise index sets for finite volume geometry
            auto numLocalFaces = intersectionMapper_.numFaces(element);
            std::vector<GridIndexType> scvfsIndexSet;
            scvfsIndexSet.reserve(numLocalFaces);
            localToGlobalScvfIndices_[eIdx].resize(numLocalFaces);
 
            std::vector<GridIndexType> neighborVolVarIndexSet;
            neighborVolVarIndexSet.reserve(numLocalFaces);
 
            for (const auto& intersection : intersections(this->gridView(), element))
            {
                const auto localFaceIndex = intersection.indexInInside();
                localToGlobalScvfIndices_[eIdx][localFaceIndex] = numScvf_;
                scvfsIndexSet.push_back(numScvf_++);
 
                if (intersection.neighbor())
                {
                    const auto nIdx = this->elementMapper().index(intersection.outside());
                    neighborVolVarIndexSet.emplace_back(nIdx);
                }
                else
                    neighborVolVarIndexSet.emplace_back(numScvs_ + numBoundaryScvf_++);
            }
 
            // Save the scvf indices belonging to this scv to build up fv element geometries fast
            scvfIndicesOfScv_[eIdx] = scvfsIndexSet;
            neighborVolVarIndices_[eIdx] = neighborVolVarIndexSet;
        }
 
        // build the connectivity map for an efficient assembly
        connectivityMap_.update(*this);
    }

    std::size_t numScvs_;
    std::size_t numScvf_;
    std::size_t numBoundaryScvf_;
    std::vector<std::vector<GridIndexType>> localToGlobalScvfIndices_;
    std::vector<std::vector<GridIndexType>> neighborVolVarIndices_;
 
    // mappers
    ConnectivityMap connectivityMap_;
    IntersectionMapper intersectionMapper_;

    std::vector<std::vector<GridIndexType>> scvfIndicesOfScv_;
};
 
} // end namespace
 
#endif