[diamond] Fix and test fvGeometry.geometry(scv,scvf) interface

dumux/discretization/facecentered/diamond/fvelementgeometry.hh

@@ -178,16 +178,39 @@ public:
     { return eIdx_; }
 
     //! Geometry of a sub control volume
-    typename SubControlVolume::Traits::Geometry geometry(const Element& element, const SubControlVolume& scv) const
+    typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const
     {
         assert(isBound());
-        const auto geo = element.geometry();
+        const auto geo = element().geometry();
         return {
-            SubControlVolume::Traits::geometryType(),
+            SubControlVolume::Traits::geometryType(geo.type()),
             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())
+        {
+            // use the information that each boundary scvf corresponds to one scv constructed around the same facet
+            const auto localFacetIndex = scvf.insideScvIdx();
+            return {
+                referenceElement(geo).type(localFacetIndex, 1),
+                GeometryHelper(geo).getBoundaryScvfCorners(localFacetIndex)
+            };
+        }
+        else
+        {
+            return {
+                SubControlVolumeFace::Traits::interiorGeometryType(geo.type()),
+                GeometryHelper(geo).getScvfCorners(scvf.index())
+            };
+        }
+    }
+
 private:
     std::optional<Element> element_;
     GridIndexType eIdx_;

test/discretization/facecentered/diamond/test_diamondgridgeometry.cc

@@ -100,6 +100,9 @@ int main (int argc, char *argv[])
         for (auto&& scv : scvs(fvGeometry))
         {
             std::cout << "-- scv " << scv.localDofIndex() << " center at: " << scv.center() << " , volume: " << scv.volume()  << std::endl;
+            const auto geo = fvGeometry.geometry(scv);
+            if ((geo.center()-scv.center()).two_norm() > 1e-14)
+                DUNE_THROW(Dune::Exception, "Center of scv-geometry and scv do not match! " << geo.center() << ", " << scv.center());
         }
 
         auto range2 = scvfs(fvGeometry);
@@ -111,6 +114,10 @@ int main (int argc, char *argv[])
         for (auto&& scvf : scvfs(fvGeometry))
         {
             std::cout << "-- scvf " << scvf.index() << " ip at: " << scvf.ipGlobal() << " normal: " << scvf.unitOuterNormal();
+            const auto geo = fvGeometry.geometry(scvf);
+            if ((geo.center()-scvf.center()).two_norm() > 1e-14)
+                DUNE_THROW(Dune::Exception, "Center of scvf-geometry and scvf do not match! " << geo.center() << ", " << scvf.center());
+
             if (scvf.boundary())
             {
                 ++boundaryCount;