From b1878faff351e19c1c8efdca26954f78fc2333e5 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Thu, 14 Jul 2022 14:49:18 +0200
Subject: [PATCH] [box] Add geometry(scv)/geometry(scvf) interfaces to the
 localview

---
 dumux/discretization/box/fvelementgeometry.hh | 70 ++++++++++++++++++-
 dumux/discretization/box/fvgridgeometry.hh    | 15 ++++
 2 files changed, 83 insertions(+), 2 deletions(-)

diff --git a/dumux/discretization/box/fvelementgeometry.hh b/dumux/discretization/box/fvelementgeometry.hh
index 4424e4c641..6b8d2c3f96 100644
--- a/dumux/discretization/box/fvelementgeometry.hh
+++ b/dumux/discretization/box/fvelementgeometry.hh
@@ -28,6 +28,10 @@
 
 #include <optional>
 #include <utility>
+#include <unordered_map>
+#include <array>
+#include <vector>
+
 #include <dune/geometry/type.hh>
 #include <dune/localfunctions/lagrange/pqkfactory.hh>
 
@@ -59,6 +63,10 @@ class BoxFVElementGeometry<GG, true>
     using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
     using CoordScalar = typename GridView::ctype;
     using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType;
+
+    using GeometryHelper = BoxGeometryHelper<GridView, dim,
+                                             typename GG::SubControlVolume,
+                                             typename GG::SubControlVolumeFace>;
 public:
     //! export the element type
     using Element = typename GridView::template Codim<0>::Entity;
@@ -185,7 +193,33 @@ public:
     bool hasBoundaryScvf() const
     { return gridGeometry().hasBoundaryScvf(eIdx_); }
 
+    //! Geometry of a sub control volume
+    typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const
+    {
+        assert(isBound());
+        const auto geo = element().geometry();
+        return { Dune::GeometryTypes::cube(dim), GeometryHelper(geo).getScvCorners(scv.indexInElement()) };
+    }
+
+    //! Geometry of a sub control volume face
+    typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const
+    {
+        assert(isBound());
+        const auto geo = element().geometry();
+        if (scvf.boundary())
+        {
+            const auto localBoundaryIndex = scvf.index() - numInnerScvf_(element());
+            const auto& key = gridGeometryPtr_->scvfBoundaryGeometryKeys(eIdx_)[localBoundaryIndex];
+            return { Dune::GeometryTypes::cube(dim-1), GeometryHelper(geo).getBoundaryScvfCorners(key[0], key[1]) };
+        }
+        else
+            return { Dune::GeometryTypes::cube(dim-1), GeometryHelper(geo).getScvfCorners(scvf.index()) };
+    }
+
 private:
+    unsigned int numInnerScvf_(const Element& element) const
+    { return (dim==1) ? 1 : element.subEntities(dim-1); }
+
     const GridGeometry* gridGeometryPtr_;
     GridIndexType eIdx_;
 
@@ -331,7 +365,32 @@ public:
     bool hasBoundaryScvf() const
     { return hasBoundaryScvf_; }
 
+    //! Geometry of a sub control volume
+    typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const
+    {
+        assert(isBound());
+        const auto geo = element().geometry();
+        return { Dune::GeometryTypes::cube(dim), GeometryHelper(geo).getScvCorners(scv.indexInElement()) };
+    }
+
+    //! Geometry of a sub control volume face
+    typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const
+    {
+        assert(isBound());
+        const auto geo = element().geometry();
+        if (scvf.boundary())
+        {
+            const auto localBoundaryIndex = scvf.index() - numInnerScvf_(element());
+            const auto& key = scvfBoundaryGeometryKeys_[localBoundaryIndex];
+            return { Dune::GeometryTypes::cube(dim-1), GeometryHelper(geo).getBoundaryScvfCorners(key[0], key[1]) };
+        }
+        else
+            return { Dune::GeometryTypes::cube(dim-1), GeometryHelper(geo).getScvfCorners(scvf.index()) };
+    }
+
 private:
+    unsigned int numInnerScvf_(const Element& element) const
+    { return (dim==1) ? 1 : element.subEntities(dim-1); }
 
     void makeElementGeometries_()
     {
@@ -349,7 +408,7 @@ private:
         scvs_.resize(elementGeometry.corners());
         for (LocalIndexType scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx)
         {
-            // get asssociated dof index
+            // get associated dof index
             const auto dofIdxGlobal = gridGeometry().vertexMapper().subIndex(element, scvLocalIdx, dim);
 
             // add scv to the local container
@@ -360,8 +419,9 @@ private:
         }
 
         // construct the sub control volume faces
-        const auto numInnerScvf = (dim==1) ? 1 : element.subEntities(dim-1);
+        const auto numInnerScvf = numInnerScvf_(element);
         scvfs_.resize(numInnerScvf);
+        scvfBoundaryGeometryKeys_.clear();
 
         LocalIndexType scvfLocalIdx = 0;
         for (; scvfLocalIdx < numInnerScvf; ++scvfLocalIdx)
@@ -400,6 +460,11 @@ private:
                                         std::move(localScvIndices),
                                         true);
 
+                    scvfBoundaryGeometryKeys_.emplace_back(std::array<LocalIndexType, 2>{{
+                        static_cast<LocalIndexType>(intersection.indexInInside()),
+                        static_cast<LocalIndexType>(isScvfLocalIdx)
+                    }});
+
                     // increment local counter
                     scvfLocalIdx++;
                 }
@@ -417,6 +482,7 @@ private:
     //! vectors to store the geometries locally after binding an element
     std::vector<SubControlVolume> scvs_;
     std::vector<SubControlVolumeFace> scvfs_;
+    std::vector<std::array<LocalIndexType, 2>> scvfBoundaryGeometryKeys_;
 
     bool hasBoundaryScvf_ = false;
 };
diff --git a/dumux/discretization/box/fvgridgeometry.hh b/dumux/discretization/box/fvgridgeometry.hh
index b26b54d3c2..0ce3c0b1dd 100644
--- a/dumux/discretization/box/fvgridgeometry.hh
+++ b/dumux/discretization/box/fvgridgeometry.hh
@@ -28,6 +28,8 @@
 
 #include <utility>
 #include <unordered_map>
+#include <array>
+#include <vector>
 
 #include <dune/localfunctions/lagrange/lagrangelfecache.hh>
 
@@ -199,11 +201,17 @@ public:
     bool hasBoundaryScvf(GridIndexType eIdx) const
     { return hasBoundaryScvf_[eIdx]; }
 
+    //! Returns a key of local indices to obtain boundary scvf geometries
+    //! Used by the local view
+    const std::vector<std::array<LocalIndexType, 2>>& scvfBoundaryGeometryKeys(GridIndexType eIdx) const
+    { return scvfBoundaryGeometryKeys_.at(eIdx); }
+
 private:
     void update_()
     {
         scvs_.clear();
         scvfs_.clear();
+        scvfBoundaryGeometryKeys_.clear();
 
         auto numElements = this->gridView().size(0);
         scvs_.resize(numElements);
@@ -287,6 +295,11 @@ private:
                                                   std::move(localScvIndices),
                                                   true);
 
+                        scvfBoundaryGeometryKeys_[eIdx].emplace_back(std::array<LocalIndexType, 2>{{
+                            static_cast<LocalIndexType>(intersection.indexInInside()),
+                            static_cast<LocalIndexType>(isScvfLocalIdx)
+                        }});
+
                         // increment local counter
                         scvfLocalIdx++;
                     }
@@ -351,6 +364,8 @@ private:
 
     std::vector<std::vector<SubControlVolume>> scvs_;
     std::vector<std::vector<SubControlVolumeFace>> scvfs_;
+    std::unordered_map<GridIndexType, std::vector<std::array<LocalIndexType, 2>>> scvfBoundaryGeometryKeys_;
+
     // TODO do we need those?
     std::size_t numScv_;
     std::size_t numScvf_;
-- 
GitLab