[disc][wmpfa] Add draft of gridgeometry and elementgeometry

dumux/discretization/cellcentered/wmpfa/fvelementgeometry.hh

+// -*- 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/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup CCWMpfaDiscretization
+ * \brief Stencil-local finite volume geometry (scvs and scvfs) for cell-centered TPFA models
+ *        This builds up the sub control volumes and sub control volume faces
+ *        for each element in the local scope we are restricting to, e.g. stencil or element.
+ */
+#ifndef DUMUX_DISCRETIZATION_CCWMPFA_FV_ELEMENT_GEOMETRY_HH
+#define DUMUX_DISCRETIZATION_CCWMPFA_FV_ELEMENT_GEOMETRY_HH
+
+#include <algorithm>
+#include <array>
+#include <vector>
+
+#include <dune/common/exceptions.hh>
+#include <dumux/common/indextraits.hh>
+#include <dune/common/iteratorrange.hh>
+#include <dumux/discretization/scvandscvfiterators.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup CCWMpfaDiscretization
+ * \brief Stencil-local finite volume geometry (scvs and scvfs) for cell-centered TPFA models
+ *        This builds up the sub control volumes and sub control volume faces
+ *        for each element in the local scope we are restricting to, e.g. stencil or element.
+ * \tparam GG the finite volume grid geometry type
+ * \tparam enableGridGeometryCache if the grid geometry is cached or not
+ * \note This class is specialized for versions with and without caching the fv geometries on the grid view
+ */
+template<class GG, bool enableGridGeometryCache>
+class CCWMpfaFVElementGeometry;
+
+/*!
+ * \ingroup CCWMpfaDiscretization
+ * \brief Stencil-local finite volume geometry (scvs and scvfs) for cell-centered TPFA models
+ *        Specialization for grid caching enabled
+ * \note The finite volume geometries are stored in the corresponding GridGeometry
+ */
+template<class GG>
+class CCWMpfaFVElementGeometry<GG, true>
+{
+    using ThisType = CCWMpfaFVElementGeometry<GG, true>;
+    using GridView = typename GG::GridView;
+    using GridIndexType = typename IndexTraits<GridView>::GridIndex;
+
+public:
+    //! export type of the element
+    using Element = typename GridView::template Codim<0>::Entity;
+    //! export type of subcontrol volume
+    using SubControlVolume = typename GG::SubControlVolume;
+    //! export type of subcontrol volume face
+    using SubControlVolumeFace = typename GG::SubControlVolumeFace;
+    //! export type of finite volume grid geometry
+    using GridGeometry = GG;
+
+    //! the maximum number of scvs per element
+    static constexpr std::size_t maxNumElementScvs = 1;
+    //! the maximum number of scvfs per element (use cubes for maximum)
+    static constexpr std::size_t maxNumElementScvfs = 2*GridView::dimension;
+
+    //! Constructor
+    CCWMpfaFVElementGeometry(const GridGeometry& gridGeometry)
+    : gridGeometryPtr_(&gridGeometry) {}
+
+    //! Get an elment sub control volume with a global scv index
+    //! We separate element and neighbor scvs to speed up mapping
+    const SubControlVolume& scv(GridIndexType scvIdx) const
+    {
+        return gridGeometry().scv(scvIdx);
+    }
+
+    //! Get an element sub control volume face with a global scvf index
+    //! We separate element and neighbor scvfs to speed up mapping
+    const SubControlVolumeFace& scvf(GridIndexType scvfIdx) const
+    {
+        return gridGeometry().scvf(scvfIdx);
+    }
+
+    //! Get the scvf on the same face but from the other side
+    //! Note that e.g. the normals might be different in the case of surface grids
+    const SubControlVolumeFace& flipScvf(GridIndexType scvfIdx, unsigned int outsideScvIdx = 0) const
+    {
+        return gridGeometry().flipScvf(scvfIdx, outsideScvIdx);
+    }
+
+    //! iterator range for sub control volumes. Iterates over
+    //! all scvs of the bound element (not including neighbor scvs)
+    //! This is a free function found by means of ADL
+    //! To iterate over all sub control volumes of this FVElementGeometry use
+    //! for (auto&& scv : scvs(fvGeometry))
+    friend inline Dune::IteratorRange< ScvIterator<SubControlVolume, std::array<GridIndexType, 1>, ThisType> >
+    scvs(const CCWMpfaFVElementGeometry& 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));
+    }
+
+    //! iterator range for sub control volumes faces. Iterates over
+    //! all scvfs of the bound element (not including neighbor scvfs)
+    //! This is a free function found by means of ADL
+    //! To iterate over all sub control volume faces of this FVElementGeometry use
+    //! for (auto&& scvf : scvfs(fvGeometry))
+    friend inline Dune::IteratorRange< ScvfIterator<SubControlVolumeFace, std::vector<GridIndexType>, ThisType> >
+    scvfs(const CCWMpfaFVElementGeometry& 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));
+    }
+
+    //! number of sub control volumes in this fv element geometry
+    std::size_t numScv() const
+    {
+        return scvIndices_.size();
+    }
+
+    //! number of sub control volumes in this fv element geometry
+    std::size_t numScvf() const
+    {
+        return gridGeometry().scvfIndicesOfScv(scvIndices_[0]).size();
+    }
+
+    //! Binding of an element, called by the local jacobian to prepare element assembly
+    void bind(const Element& element)
+    {
+        this->bindElement(element);
+    }
+
+    //! Bind only element-local
+    void bindElement(const Element& element)
+    {
+        elementPtr_ = &element;
+        scvIndices_[0] = gridGeometry().elementMapper().index(*elementPtr_);
+    }
+
+    //! The global finite volume geometry we are a restriction of
+    const GridGeometry& gridGeometry() const
+    { return *gridGeometryPtr_; }
+
+    //! Returns whether one of the geometry's scvfs lies on a boundary
+    bool hasBoundaryScvf() const
+    { return gridGeometry().hasBoundaryScvf(scvIndices_[0]); }
+
+private:
+
+    const Element* elementPtr_;
+    std::array<GridIndexType, 1> scvIndices_;
+    const GridGeometry* gridGeometryPtr_;
+};
+
+} // end namespace Dumux
+
+#endif

dumux/discretization/cellcentered/wmpfa/fvgridgeometry.hh

+// -*- 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/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup CCWMpfaDiscretization
+ * \brief The finite volume geometry (scvs and scvfs) for cell-centered TPFA models on a grid view
+ *        This builds up the sub control volumes and sub control volume faces
+ *        for each element of the grid partition.
+ */
+#ifndef DUMUX_DISCRETIZATION_CCWMpfa_FV_GRID_GEOMETRY_HH
+#define DUMUX_DISCRETIZATION_CCWMpfa_FV_GRID_GEOMETRY_HH
+
+#include <algorithm>
+
+#include <dumux/common/indextraits.hh>
+#include <dumux/common/defaultmappertraits.hh>
+
+#include <dumux/discretization/method.hh>
+#include <dumux/discretization/basegridgeometry.hh>
+#include <dumux/discretization/checkoverlapsize.hh>
+#include <dumux/discretization/cellcentered/subcontrolvolume.hh>
+#include <dumux/discretization/cellcentered/connectivitymap.hh>
+#include <dumux/discretization/cellcentered/wmpfa/fvelementgeometry.hh>
+#include <dumux/discretization/cellcentered/wmpfa/subcontrolvolumeface.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup CCWMpfaDiscretization
+ * \brief The default traits for the tpfa finite volume grid geometry
+ *        Defines the scv and scvf types and the mapper types
+ * \tparam the grid view type
+ */
+template<class GridView, class MapperTraits = DefaultMapperTraits<GridView>>
+struct CCWMpfaDefaultGridGeometryTraits
+: public MapperTraits
+{
+    using SubControlVolume = CCSubControlVolume<GridView>;
+    using SubControlVolumeFace = CCWMpfaSubControlVolumeFace<GridView>;
+
+    template<class GridGeometry>
+    using ConnectivityMap = CCSimpleConnectivityMap<GridGeometry>;
+
+    template<class GridGeometry, bool enableCache>
+    using LocalView = CCWMpfaFVElementGeometry<GridGeometry, enableCache>;
+
+    //! State the maximum admissible number of neighbors per scvf
+    //! Per default, we allow for 8 branches on network/surface grids, where
+    //! conformity is assumed. For normal grids, we allow a maximum of one
+    //! hanging node per scvf. Use different traits if you need more.
+    static constexpr int maxNumScvfNeighbors = int(GridView::dimension)<int(GridView::dimensionworld) ? 8 : 1<<(GridView::dimension-1);
+};
+
+/*!
+ * \ingroup CCWMpfaDiscretization
+ * \brief The finite volume geometry (scvs and scvfs) for cell-centered weighted MPFA models on a grid view
+ *        This builds up the sub control volumes and sub control volume faces
+ * \note This class is specialized for versions with and without caching the fv geometries on the grid view
+ */
+template<class GridView,
+         bool enableGridGeometryCache = false,
+         class Traits = CCWMpfaDefaultGridGeometryTraits<GridView> >
+class CCWMpfaFVGridGeometry;
+
+/*!
+ * \ingroup CCWMpfaDiscretization
+ * \brief The finite volume geometry (scvs and scvfs) for cell-centered weighted MPFA models on a grid view
+ *        This builds up the sub control volumes and sub control volume faces
+ * \note For caching enabled we store the fv geometries for the whole grid view which is memory intensive but faster
+ */
+template<class GV, class Traits>
+class CCWMpfaFVGridGeometry<GV, true, Traits>
+: public BaseGridGeometry<GV, Traits>
+{
+    using ThisType = CCWMpfaFVGridGeometry<GV, true, Traits>;
+    using ParentType = BaseGridGeometry<GV, Traits>;
+    using ConnectivityMap = typename Traits::template ConnectivityMap<ThisType>;
+    using GridIndexType = typename IndexTraits<GV>::GridIndex;
+    using LocalIndexType = typename IndexTraits<GV>::LocalIndex;
+    using Element = typename GV::template Codim<0>::Entity;
+
+    static const int dim = GV::dimension;
+    static const int dimWorld = GV::dimensionworld;
+
+public:
+    //! export the type of the fv element geometry (the local view type)
+    using LocalView = typename Traits::template LocalView<ThisType, true>;
+    //! export the type of sub control volume
+    using SubControlVolume = typename Traits::SubControlVolume;
+    //! export the type of sub control volume
+    using SubControlVolumeFace = typename Traits::SubControlVolumeFace;
+    //! export dof mapper type
+    using DofMapper = typename Traits::ElementMapper;
+
+    //! export the discretization method this geometry belongs to
+    static constexpr DiscretizationMethod discMethod = DiscretizationMethod::ccwmpfa;
+
+    //! The maximum admissible stencil size (used for static memory allocation during assembly)
+    static constexpr int maxElementStencilSize = LocalView::maxNumElementScvfs*Traits::maxNumScvfNeighbors + 1;
+
+    //! export the grid view type
+    using GridView = GV;
+
+    //! Constructor
+    CCWMpfaFVGridGeometry(const GridView& gridView)
+    : ParentType(gridView)
+    {
+        // Check if the overlap size is what we expect
+        if (!CheckOverlapSize<DiscretizationMethod::ccwmpfa>::isValid(gridView))
+            DUNE_THROW(Dune::InvalidStateException, "The ccwmpfa discretization method needs at least an overlap of 1 for parallel computations. "
+                                                     << " Set the parameter \"Grid.Overlap\" in the input file.");
+    }
+
+    //! the element mapper is the dofMapper
+    //! this is convenience to have better chance to have the same main files for box/tpfa/mpfa...
+    const DofMapper& dofMapper() const
+    { return this->elementMapper(); }
+
+    //! The total number of sub control volumes
+    std::size_t numScv() const
+    {
+        return scvs_.size();
+    }
+
+    //! The total number of sub control volume faces
+    std::size_t numScvf() const
+    {
+        return scvfs_.size();
+    }
+
+    //! The total number of boundary sub control volume faces
+    std::size_t numBoundaryScvf() const
+    {
+        return numBoundaryScvf_;
+    }
+
+    //! The total number of degrees of freedom
+    std::size_t numDofs() const
+    { return this->gridView().size(0); }
+
+    //! update all fvElementGeometries (do this again after grid adaption)
+    void update()
+    {
+        ParentType::update();
+
+        // clear containers (necessary after grid refinement)
+        scvs_.clear();
+        scvfs_.clear();
+        scvfIndicesOfScv_.clear();
+        flipScvfIndices_.clear();
+
+        // determine size of containers
+        std::size_t numScvs = numDofs();
+        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);
+        hasBoundaryScvf_.resize(numScvs, false);
+
+        // Build the scvs and scv faces
+        GridIndexType scvfIdx = 0;
+        numBoundaryScvf_ = 0;
+        for (const auto& element : elements(this->gridView()))
+        {
+            const auto eIdx = this->elementMapper().index(element);
+            scvs_[eIdx] = SubControlVolume(element.geometry(), eIdx);
+
+            // the element-wise index sets for finite volume geometry
+            std::vector<GridIndexType> scvfsIndexSet;
+            scvfsIndexSet.reserve(element.subEntities(1));
+
+            // for network grids there might be multiple intersection with the same geometryInInside
+            // we indentify those by the indexInInside for now (assumes conforming grids at branching facets)
+            using ScvfGridIndexStorage = typename SubControlVolumeFace::Traits::GridIndexStorage;
+
+            LocalIndexType localScvfIdx = 0;
+            for (const auto& intersection : intersections(this->gridView(), element))
+            {
+                // inner sub control volume faces (includes periodic boundaries)
+                if (intersection.neighbor())
+                {
+                    // update the grid geometry if we have periodic boundaries
+                    if (intersection.boundary())
+                        this->setPeriodic();
+
+                    if (dim == dimWorld)
+                    {
+                        const auto nIdx = this->elementMapper().index(intersection.outside());
+                        scvfs_.emplace_back(intersection,
+                                            intersection.geometry(),
+                                            scvfIdx,
+                                            localScvfIdx,
+                                            ScvfGridIndexStorage({eIdx, nIdx}),
+                                            false);
+                        scvfsIndexSet.push_back(scvfIdx++);
+                    }
+                    // this is for network grids
+                    // (will be optimized away of dim == dimWorld)
+                    else
+                        DUNE_THROW(Dune::NotImplemented, "Weighted MPFA schemes not implemented for network grids");
+                }
+                // boundary sub control volume faces
+                else if (intersection.boundary())
+                {
+                    scvfs_.emplace_back(intersection,
+                                        intersection.geometry(),
+                                        scvfIdx,
+                                        localScvfIdx,
+                                        ScvfGridIndexStorage({eIdx, static_cast<GridIndexType>(this->gridView().size(0) + numBoundaryScvf_++)}),
+                                        true);
+                    scvfsIndexSet.push_back(scvfIdx++);
+
+                    hasBoundaryScvf_[eIdx] = true;
+                }
+                ++localScvfIdx;
+            }
+
+            // Save the scvf indices belonging to this scv to build up fv element geometries fast
+            scvfIndicesOfScv_[eIdx] = scvfsIndexSet;
+        }
+
+        // Make the flip index set for network, surface, and periodic grids
+        if (this->isPeriodic())
+        {
+            flipScvfIndices_.resize(scvfs_.size());
+            for (auto&& scvf : scvfs_)
+            {
+                if (scvf.boundary())
+                    continue;
+
+                flipScvfIndices_[scvf.index()].resize(scvf.numOutsideScvs());
+                const auto insideScvIdx = scvf.insideScvIdx();
+                // check which outside scvf has the insideScvIdx index in its outsideScvIndices
+                for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i)
+                    flipScvfIndices_[scvf.index()][i] = findFlippedScvfIndex_(insideScvIdx, scvf.outsideScvIdx(i));
+            }
+        }
+
+        // build the connectivity map for an effecient assembly
+        connectivityMap_.update(*this);
+    }
+
+    //! Get a sub control volume with a global scv index
+    const SubControlVolume& scv(GridIndexType scvIdx) const
+    {
+        return scvs_[scvIdx];
+    }
+
+    //! Get a sub control volume face with a global scvf index
+    const SubControlVolumeFace& scvf(GridIndexType scvfIdx) const
+    {
+        return scvfs_[scvfIdx];
+    }
+
+    //! Get the scvf on the same face but from the other side
+    //! Note that e.g. the normals might be different in the case of surface grids
+    const SubControlVolumeFace& flipScvf(GridIndexType scvfIdx, unsigned int outsideScvfIdx = 0) const
+    {
+        return scvfs_[flipScvfIndices_[scvfIdx][outsideScvfIdx]];
+    }
+
+    //! Get the sub control volume face indices of an scv by global index
+    const std::vector<GridIndexType>& scvfIndicesOfScv(GridIndexType scvIdx) const
+    {
+        return scvfIndicesOfScv_[scvIdx];
+    }
+
+    /*!
+     * \brief Returns the connectivity map of which dofs have derivatives with respect
+     *        to a given dof.
+     */
+    const ConnectivityMap &connectivityMap() const
+    { return connectivityMap_; }
+
+    //! Returns whether one of the geometry's scvfs lies on a boundary
+    bool hasBoundaryScvf(GridIndexType eIdx) const
+    { return hasBoundaryScvf_[eIdx]; }
+
+private:
+    // find the scvf that has insideScvIdx in its outsideScvIdx list and outsideScvIdx as its insideScvIdx
+    GridIndexType findFlippedScvfIndex_(GridIndexType insideScvIdx, GridIndexType outsideScvIdx)
+    {
+        // go over all potential scvfs of the outside scv
+        for (auto outsideScvfIndex : scvfIndicesOfScv_[outsideScvIdx])
+        {
+            const auto& outsideScvf = this->scvf(outsideScvfIndex);
+            for (unsigned int j = 0; j < outsideScvf.numOutsideScvs(); ++j)
+                if (outsideScvf.outsideScvIdx(j) == insideScvIdx)
+                    return outsideScvf.index();
+        }
+
+        DUNE_THROW(Dune::InvalidStateException, "No flipped version of this scvf found!");
+    }
+
+    //! connectivity map for efficient assembly
+    ConnectivityMap connectivityMap_;
+
+    //! containers storing the global data
+    std::vector<SubControlVolume> scvs_;
+    std::vector<SubControlVolumeFace> scvfs_;
+    std::vector<std::vector<GridIndexType>> scvfIndicesOfScv_;
+    std::size_t numBoundaryScvf_;
+    std::vector<bool> hasBoundaryScvf_;
+
+    //! needed for embedded surface and network grids (dim < dimWorld)
+    std::vector<std::vector<GridIndexType>> flipScvfIndices_;
+};
+
+} // end namespace Dumux
+
+#endif