Merge branch 'feature/localgeometry-element-interface' into 'master'

[discretization] Add interfaces element/isBound to grid geometry local views

See merge request !2509

CHANGELOG.md

@@ -32,6 +32,11 @@ auto fluidMatrixInteractionAtPos(const GlobalPosition &globalPos) const
     return makeFluidMatrixInteraction(FluidMatrix::MPAdapter(pcKrSwCurve_));
 }
 ```
+- __Grid-geometry__: The local views of grid geometries are now required to implement the interfaces
+`element()` (returning the bound element) and `isBound()` returning a `bool` which is `true` if the
+functions `bind` or `bindElement` have been called (i.e. the local geometry is in a bound state). These
+interfaces are currently not used (except in the unit tests) but will be required
+by the assembler in future Dumux versions.
 
 Differences Between DuMu<sup>x</sup> 3.3 and DuMu<sup>x</sup> 3.2
 =============================================

dumux/discretization/box/fvelementgeometry.hh

@@ -26,6 +26,7 @@
 #ifndef DUMUX_DISCRETIZATION_BOX_FV_ELEMENT_GEOMETRY_HH
 #define DUMUX_DISCRETIZATION_BOX_FV_ELEMENT_GEOMETRY_HH
 
+#include <optional>
 #include <dune/geometry/type.hh>
 #include <dune/localfunctions/lagrange/pqkfactory.hh>
 
@@ -114,7 +115,7 @@ public:
     //! Get a local finite element basis
     const FeLocalBasis& feLocalBasis() const
     {
-        return gridGeometry().feCache().get(elemGeometryType_).localBasis();
+        return gridGeometry().feCache().get(element_->type()).localBasis();
     }
 
     //! The total number of sub control volumes
@@ -142,10 +143,19 @@ public:
     //! For compatibility reasons with the FVGeometry cache being disabled
     void bindElement(const Element& element)
     {
-        elemGeometryType_ = element.type();
+        element_ = element;
+        // cache element index
         eIdx_ = gridGeometry().elementMapper().index(element);
     }
 
+    //! Returns true if bind/bindElement has already been called
+    bool isBound() const
+    { return static_cast<bool>(element_); }
+
+    //! The bound element
+    const Element& element() const
+    { return *element_; }
+
     //! The grid geometry we are a restriction of
     const GridGeometry& gridGeometry() const
     { return *gridGeometryPtr_; }
@@ -155,10 +165,10 @@ public:
     { return gridGeometry().hasBoundaryScvf(eIdx_); }
 
 private:
-    Dune::GeometryType elemGeometryType_;
     const GridGeometry* gridGeometryPtr_;
-
     GridIndexType eIdx_;
+
+    std::optional<Element> element_;
 };
 
 //! specialization in case the FVElementGeometries are not stored
@@ -231,7 +241,7 @@ public:
     //! Get a local finite element basis
     const FeLocalBasis& feLocalBasis() const
     {
-        return gridGeometry().feCache().get(elemGeometryType_).localBasis();
+        return gridGeometry().feCache().get(element_->type()).localBasis();
     }
 
     //! The total number of sub control volumes
@@ -259,10 +269,19 @@ public:
     //! For compatibility reasons with the FVGeometry cache being disabled
     void bindElement(const Element& element)
     {
+        element_ = element;
         eIdx_ = gridGeometry().elementMapper().index(element);
-        makeElementGeometries(element);
+        makeElementGeometries_();
     }
 
+    //! Returns true if bind/bindElement has already been called
+    bool isBound() const
+    { return static_cast<bool>(element_); }
+
+    //! The bound element
+    const Element& element() const
+    { return *element_; }
+
     //! The grid geometry we are a restriction of
     const GridGeometry& gridGeometry() const
     { return *gridGeometryPtr_; }
@@ -273,13 +292,13 @@ public:
 
 private:
 
-    void makeElementGeometries(const Element& element)
+    void makeElementGeometries_()
     {
         hasBoundaryScvf_ = false;
 
         // get the element geometry
-        auto elementGeometry = element.geometry();
-        elemGeometryType_ = elementGeometry.type();
+        const auto& element = *element_;
+        const auto elementGeometry = element.geometry();
         const auto refElement = referenceElement(elementGeometry);
 
         // get the sub control volume geometries of this element
@@ -348,8 +367,8 @@ private:
     }
 
     //! The bound element
-    Dune::GeometryType elemGeometryType_;
     GridIndexType eIdx_;
+    std::optional<Element> element_;
 
     //! The global geometry this is a restriction of
     const GridGeometry* gridGeometryPtr_;

dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh

@@ -26,6 +26,7 @@
 #ifndef DUMUX_DISCRETIZATION_CCMPFA_FV_ELEMENT_GEOMETRY_HH
 #define DUMUX_DISCRETIZATION_CCMPFA_FV_ELEMENT_GEOMETRY_HH
 
+#include <optional>
 #include <dune/common/exceptions.hh>
 #include <dune/common/iteratorrange.hh>
 
@@ -148,9 +149,18 @@ public:
     //! Bind only element-local
     void bindElement(const Element& element)
     {
+        element_ = element;
         scvIndices_[0] = gridGeometry().elementMapper().index(element);
     }
 
+    //! Returns true if bind/bindElement has already been called
+    bool isBound() const
+    { return static_cast<bool>(element_); }
+
+    //! The bound element
+    const Element& element() const
+    { return *element_; }
+
     //! The global finite volume geometry we are a restriction of
     const GridGeometry& gridGeometry() const
     { return *gridGeometryPtr_; }
@@ -161,6 +171,7 @@ public:
 
 private:
 
+    std::optional<Element> element_;
     std::array<GridIndexType, 1> scvIndices_;
     const GridGeometry* gridGeometryPtr_;
 };
@@ -307,9 +318,18 @@ public:
     void bindElement(const Element& element)
     {
         clear();
+        element_ = element;
         makeElementGeometries(element);
     }
 
+    //! Returns true if bind/bindElement has already been called
+    bool isBound() const
+    { return static_cast<bool>(element_); }
+
+    //! The bound element
+    const Element& element() const
+    { return *element_; }
+
     //! The global finite volume geometry we are a restriction of
     const GridGeometry& gridGeometry() const
     { return *gridGeometryPtr_; }
@@ -530,6 +550,7 @@ private:
     }
 
     const GridGeometry* gridGeometryPtr_;
+    std::optional<Element> element_;
 
     // local storage after binding an element
     std::array<GridIndexType, 1> scvIndices_;

dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh

@@ -26,6 +26,7 @@
 #ifndef DUMUX_DISCRETIZATION_CCTPFA_FV_ELEMENT_GEOMETRY_HH
 #define DUMUX_DISCRETIZATION_CCTPFA_FV_ELEMENT_GEOMETRY_HH
 
+#include <optional>
 #include <algorithm>
 #include <array>
 #include <vector>
@@ -151,10 +152,18 @@ public:
     //! Bind only element-local
     void bindElement(const Element& element)
     {
-        elementPtr_ = &element;
-        scvIndices_[0] = gridGeometry().elementMapper().index(*elementPtr_);
+        element_ = element;
+        scvIndices_[0] = gridGeometry().elementMapper().index(*element_);
     }
 
+    //! Returns true if bind/bindElement has already been called
+    bool isBound() const
+    { return static_cast<bool>(element_); }
+
+    //! The bound element
+    const Element& element() const
+    { return *element_; }
+
     //! The global finite volume geometry we are a restriction of
     const GridGeometry& gridGeometry() const
     { return *gridGeometryPtr_; }
@@ -165,7 +174,7 @@ public:
 
 private:
 
-    const Element* elementPtr_;
+    std::optional<Element> element_;
     std::array<GridIndexType, 1> scvIndices_;
     const GridGeometry* gridGeometryPtr_;
 };
@@ -350,12 +359,20 @@ public:
     void bindElement(const Element& element)
     {
         clear();
-        elementPtr_ = &element;
+        element_ = element;
         scvfs_.reserve(element.subEntities(1));
         scvfIndices_.reserve(element.subEntities(1));
         makeElementGeometries(element);
     }
 
+    //! Returns true if bind/bindElement has already been called
+    bool isBound() const
+    { return static_cast<bool>(element_); }
+
+    //! The bound element
+    const Element& element() const
+    { return *element_; }
+
     //! The global finite volume geometry we are a restriction of
     const GridGeometry& gridGeometry() const
     { return *gridGeometryPtr_; }
@@ -479,7 +496,7 @@ private:
                     // only create subcontrol faces where the outside element is the bound element
                     if (dim == dimWorld)
                     {
-                        if (scvfNeighborVolVarIndices[0] == gridGeometry().elementMapper().index(*elementPtr_))
+                        if (scvfNeighborVolVarIndices[0] == gridGeometry().elementMapper().index(*element_))
                         {
                             ScvfGridIndexStorage scvIndices({eIdx, scvfNeighborVolVarIndices[0]});
                             neighborScvfs_.emplace_back(intersection,
@@ -499,7 +516,7 @@ private:
                     {
                         for (unsigned outsideScvIdx = 0; outsideScvIdx < scvfNeighborVolVarIndices.size(); ++outsideScvIdx)
                         {
-                            if (scvfNeighborVolVarIndices[outsideScvIdx] == gridGeometry().elementMapper().index(*elementPtr_))
+                            if (scvfNeighborVolVarIndices[outsideScvIdx] == gridGeometry().elementMapper().index(*element_))
                             {
                                 ScvfGridIndexStorage scvIndices;
                                 scvIndices.resize(scvfNeighborVolVarIndices.size() + 1);
@@ -560,7 +577,7 @@ private:
         hasBoundaryScvf_ = false;
     }
 
-    const Element* elementPtr_; //!< the element to which this fvgeometry is bound
+    std::optional<Element> element_; //!< the element to which this fvgeometry is bound
     const GridGeometry* gridGeometryPtr_;  //!< the grid fvgeometry
 
     // local storage after binding an element

dumux/discretization/fem/feelementgeometry.hh

@@ -62,6 +62,14 @@ public:
     void bindElement(const Element& element)
     { bind(element); }
 
+    //! Returns true if bind/bindElement has already been called
+    bool isBound() const
+    { return feBasisLocalView_.isBound(); }
+
+    //! The bound element
+    const Element& element() const
+    { return feBasisLocalView_.element(); }
+
     //! Return the finite element basis local view
     const FEBasisLocalView& feBasisLocalView() const
     { return feBasisLocalView_; }

dumux/discretization/porenetwork/fvelementgeometry.hh

@@ -24,6 +24,7 @@
 #ifndef DUMUX_DISCRETIZATION_PNM_FV_ELEMENT_GEOMETRY_HH
 #define DUMUX_DISCRETIZATION_PNM_FV_ELEMENT_GEOMETRY_HH
 
+#include <optional>
 #include <dumux/common/indextraits.hh>
 #include <dumux/discretization/scvandscvfiterators.hh>
 
@@ -103,7 +104,7 @@ public:
     //! Get a local finite element basis
     const FeLocalBasis& feLocalBasis() const
     {
-        return gridGeometry().feCache().get(elementPtr_->geometry().type()).localBasis();
+        return gridGeometry().feCache().get(element_->type()).localBasis();
     }
 
     //! The total number of sub control volumes
@@ -131,10 +132,18 @@ public:
     //! For compatibility reasons with the FVGeometry cache being disabled
     void bindElement(const Element& element)
     {
-        elementPtr_ = &element;
+        element_ = element;
         eIdx_ = gridGeometry().elementMapper().index(element);
     }
 
+    //! Returns true if bind/bindElement has already been called
+    bool isBound() const
+    { return static_cast<bool>(element_); }
+
+    //! The bound element
+    const Element& element() const
+    { return *element_; }
+
     //! The global finite volume geometry we are a restriction of
     const GridGeometry& gridGeometry() const
     { return *gridGeometryPtr_; }
@@ -144,7 +153,7 @@ public:
     { return gridGeometry().hasBoundaryScvf(eIdx_); }
 
 private:
-    const Element* elementPtr_;
+    std::optional<Element> element_;
     const GridGeometry* gridGeometryPtr_;
 
     GridIndexType eIdx_;
@@ -214,7 +223,7 @@ public:
     //! Get a local finite element basis
     const FeLocalBasis& feLocalBasis() const
     {
-        return gridGeometry().feCache().get(elementPtr_->geometry().type()).localBasis();
+        return gridGeometry().feCache().get(element_->type()).localBasis();
     }
 
     //! The total number of sub control volumes
@@ -242,11 +251,19 @@ public:
     //! For compatibility reasons with the FVGeometry cache being disabled
     void bindElement(const Element& element)
     {
-        elementPtr_ = &element;
+        element_ = element;
         eIdx_ = gridGeometry().elementMapper().index(element);
         makeElementGeometries(element);
     }
 
+    //! Returns true if bind/bindElement has already been called
+    bool isBound() const
+    { return static_cast<bool>(element_); }
+
+    //! The bound element
+    const Element& element() const
+    { return *element_; }
+
     //! The global finite volume geometry we are a restriction of
     const GridGeometry& gridGeometry() const
     { return *gridGeometryPtr_; }
@@ -299,7 +316,7 @@ private:
     }
 
     //! The bound element
-    const Element* elementPtr_;
+    std::optional<Element> element_;
     GridIndexType eIdx_;
 
     //! The global geometry this is a restriction of

dumux/discretization/staggered/fvelementgeometry.hh

@@ -24,6 +24,7 @@
 #ifndef DUMUX_DISCRETIZATION_STAGGERED_FV_ELEMENT_GEOMETRY_HH
 #define DUMUX_DISCRETIZATION_STAGGERED_FV_ELEMENT_GEOMETRY_HH
 
+#include <optional>
 #include <dumux/common/indextraits.hh>
 #include <dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh>
 
@@ -205,12 +206,20 @@ public:
     void bindElement(const Element& element)
     {
         clear_();
-        elementPtr_ = &element;
+        element_ = element;
         scvfs_.reserve(element.subEntities(1));
         scvfIndices_.reserve(element.subEntities(1));
-        makeElementGeometries_(element);
+        makeElementGeometries_();
     }
 
+    //! Returns true if bind/bindElement has already been called
+    bool isBound() const
+    { return static_cast<bool>(element_); }
+
+    //! The bound element
+    const Element& element() const
+    { return *element_; }
+
     //! The grid finite volume geometry we are a restriction of
     const GridGeometry& gridGeometry() const
     { return *gridGeometryPtr_; }
@@ -222,8 +231,9 @@ public:
 private:
 
     //! create scvs and scvfs of the bound element
-    void makeElementGeometries_(const Element& element)
+    void makeElementGeometries_()
     {
+        const auto& element = *element_;
         const auto eIdx = gridGeometry().elementMapper().index(element);
         scvs_[0] = SubControlVolume(element.geometry(), eIdx);
         scvIndices_[0] = eIdx;
@@ -279,7 +289,7 @@ private:
             if (intersection.neighbor())
             {
                 // only create subcontrol faces where the outside element is the bound element
-                if (intersection.outside() == *elementPtr_)
+                if (intersection.outside() == *element_)
                 {
                     std::vector<GridIndexType> scvIndices{eIdx, scvfNeighborVolVarIndex};
                     neighborScvfs_.emplace_back(intersection,
@@ -319,7 +329,7 @@ private:
         hasBoundaryScvf_ = false;
     }
 
-    const Element* elementPtr_; //!< the element to which this fvgeometry is bound
+    std::optional<Element> element_; //!< the element to which this fvgeometry is bound
     const GridGeometry* gridGeometryPtr_;  //!< the grid fvgeometry
 
     // local storage after binding an element

dumux/multidomain/facet/box/fvelementgeometry.hh

@@ -25,6 +25,7 @@
 #define DUMUX_FACETCOUPLING_BOX_FV_ELEMENT_GEOMETRY_HH
 
 #include <algorithm>
+#include <optional>
 
 #include <dune/geometry/type.hh>
 
@@ -110,7 +111,7 @@ public:
 
     //! Get a local finite element basis
     const FeLocalBasis& feLocalBasis() const
-    { return gridGeometry().feCache().get(elemGeometryType_).localBasis(); }
+    { return gridGeometry().feCache().get(element_->type()).localBasis(); }
 
     //! The total number of sub control volumes
     std::size_t numScv() const
@@ -133,7 +134,7 @@ public:
     //! For compatibility reasons with the FVGeometry cache being disabled
     void bindElement(const Element& element)
     {
-        elemGeometryType_ = element.type();
+        element_ = element;
         eIdx_ = gridGeometry().elementMapper().index(element);
     }
 
@@ -141,11 +142,19 @@ public:
     const GridGeometry& gridGeometry() const
     { return *gridGeometryPtr_; }
 
+    //! Returns true if bind/bindElement has already been called
+    bool isBound() const
+    { return static_cast<bool>(element_); }
+
+    //! The bound element
+    const Element& element() const
+    { return *element_; }
+
 private:
-    Dune::GeometryType elemGeometryType_;
     const GridGeometry* gridGeometryPtr_;
 
     GridIndexType eIdx_;
+    std::optional<Element> element_;
 };
 
 //! specialization in case the geometries are not stored grid-wide
@@ -216,7 +225,7 @@ public:
 
     //! Get a local finite element basis
     const FeLocalBasis& feLocalBasis() const
-    { return gridGeometry().feCache().get(elemGeometryType_).localBasis(); }
+    { return gridGeometry().feCache().get(element_->type()).localBasis(); }
 
     //! The total number of sub control volumes
     std::size_t numScv() const
@@ -239,23 +248,30 @@ public:
     //! For compatibility reasons with the FVGeometry cache being disabled
     void bindElement(const Element& element)
     {
+        element_ = element;
         eIdx_ = gridGeometry().elementMapper().index(element);
-        makeElementGeometries(element);
+        makeElementGeometries_();
     }
 
     //! The global finite volume geometry we are a restriction of
     const GridGeometry& gridGeometry() const
     { return *gridGeometryPtr_; }
 
+    //! Returns true if bind/bindElement has already been called
+    bool isBound() const
+    { return static_cast<bool>(element_); }
+
+    //! The bound element
+    const Element& element() const
+    { return *element_; }
+
 private:
 
-    void makeElementGeometries(const Element& element)
+    void makeElementGeometries_()
     {
-        auto eIdx = gridGeometry().elementMapper().index(element);
-
         // get the element geometry
-        auto elementGeometry = element.geometry();
-        elemGeometryType_ = elementGeometry.type();
+        const auto& element = *element_;
+        const auto elementGeometry = element.geometry();
         const auto refElement = referenceElement(elementGeometry);
 
         // get the sub control volume geometries of this element
@@ -268,7 +284,7 @@ private:
         for (LocalIndexType scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx)
             scvs_.emplace_back(geometryHelper,
                                scvLocalIdx,
-                               eIdx,
+                               eIdx_,
                                gridGeometry().vertexMapper().subIndex(element, scvLocalIdx, dim));
 
         // construct the sub control volume faces
@@ -338,13 +354,13 @@ private:
         }
     }
 
-    //! The bound element
-    Dune::GeometryType elemGeometryType_;
-    GridIndexType eIdx_;
-
     //! The global geometry this is a restriction of
     const GridGeometry* gridGeometryPtr_;
 
+    //! The bound element
+    GridIndexType eIdx_;
+    std::optional<Element> element_;
+
     //! vectors to store the geometries locally after binding an element
     std::vector<SubControlVolume> scvs_;
     std::vector<SubControlVolumeFace> scvfs_;

dumux/porousmediumflow/boxdfm/fvelementgeometry.hh

@@ -29,6 +29,8 @@
 #ifndef DUMUX_POROUSMEDIUMFLOW_BOXDFM_FV_ELEMENT_GEOMETRY_HH
 #define DUMUX_POROUSMEDIUMFLOW_BOXDFM_FV_ELEMENT_GEOMETRY_HH
 
+#include <optional>
+
 #include <dune/common/version.hh>
 #include <dune/geometry/type.hh>
 #include <dune/localfunctions/lagrange/pqkfactory.hh>
@@ -119,7 +121,7 @@ public:
 
     //! Get a local finite element basis
     const FeLocalBasis& feLocalBasis() const
-    { return gridGeometry().feCache().get(elemGeometryType_).localBasis(); }
+    { return gridGeometry().feCache().get(element_->type()).localBasis(); }
 
     //! The total number of sub control volumes
     std::size_t numScv() const
@@ -144,7 +146,7 @@ public:
     */
     void bindElement(const Element& element)
     {
-        elemGeometryType_ = element.type();
+        element_ = element;
         eIdx_ = gridGeometry().elementMapper().index(element);
     }
 
@@ -152,9 +154,18 @@ public:
     const GridGeometry& gridGeometry() const
     { return *gridGeometryPtr_; }
 
+    //! Returns true if bind/bindElement has already been called
+    bool isBound() const
+    { return static_cast<bool>(element_); }
+
+    //! The bound element
+    const Element& element() const
+    { return *element_; }
+
 private:
-    Dune::GeometryType elemGeometryType_;
     const GridGeometry* gridGeometryPtr_;
+
+    std::optional<Element> element_;
     GridIndexType eIdx_;
 };
 
@@ -229,7 +240,7 @@ public:
 
     //! Get a local finite element basis
     const FeLocalBasis& feLocalBasis() const
-    { return gridGeometry().feCache().get(elemGeometryType_).localBasis(); }
+    { return gridGeometry().feCache().get(element_->type()).localBasis(); }
 
     //! The total number of sub control volumes
     std::size_t numScv() const
@@ -239,8 +250,12 @@ public:
     std::size_t numScvf() const
     { return scvfs_.size(); }
 
-    //! This function is for compatibility reasons with cc methods
-    //! The box stencil is always element-local so bind and bindElement are identical.
+    /*!
+     * \brief Binding of an element, has to be called before using the fvgeometries
+     *        Prepares all the volume variables within the element.
+     * \note For the box scheme, bind() and bindElement() are identical, but the
+     *       distinction is here for the sake of compatibility with cc schemes.
+     */
     void bind(const Element& element)
     {
         this->bindElement(element);
@@ -248,29 +263,34 @@ public: