Merge branch 'fix/box-scvf-construction' into 'master'

Fix/box scvf construction

See merge request !924

dumux/discretization/box/boxgeometryhelper.hh

@@ -81,7 +81,8 @@ public:
     }
 
     //! Create the sub control volume face geometries on the boundary
-    ScvfCornerStorage getBoundaryScvfCorners(const typename Intersection::Geometry& geometry,
+    ScvfCornerStorage getBoundaryScvfCorners(const Intersection& is,
+                                             const typename Intersection::Geometry& geometry,
                                              unsigned int indexInIntersection) const
     {
         return ScvfCornerStorage{{geometry.corner(0)}};
@@ -126,8 +127,9 @@ class BoxGeometryHelper<GridView, 2, ScvType, ScvfType>
     using Element = typename GridView::template Codim<0>::Entity;
     using Intersection = typename GridView::Intersection;
 
-    using ReferenceElements = typename Dune::ReferenceElements<Scalar, GridView::dimension>;
-    static constexpr int dimWorld = GridView::dimensionworld;
+    static constexpr auto dim = GridView::dimension;
+    static constexpr auto dimWorld = GridView::dimensionworld;
+    using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
 
     //! the maximum number of helper points used to construct the geometries
     //! Using a statically sized point array is much faster than dynamic allocation
@@ -136,9 +138,9 @@ class BoxGeometryHelper<GridView, 2, ScvType, ScvfType>
 public:
 
     BoxGeometryHelper(const typename Element::Geometry& geometry)
-    : elementGeometry_(geometry), corners_(geometry.corners())
+    : elementGeometry_(geometry)
+    , corners_(geometry.corners())
     {
-        // extract the corners of the sub control volumes
 #if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
         const auto referenceElement = ReferenceElements::general(geometry.type());
 #else
@@ -199,8 +201,8 @@ public:
                                       p[map[localScvIdx][3]]} };
         }
         default:
-            DUNE_THROW(Dune::NotImplemented, "Box scv geometries for dim=" << GridView::dimension
-                                                            << " dimWorld=" << GridView::dimensionworld
+            DUNE_THROW(Dune::NotImplemented, "Box scv geometries for dim=" << dim
+                                                            << " dimWorld=" << dimWorld
                                                             << " corners=" << corners_);
         }
     }
@@ -242,20 +244,28 @@ public:
                                        p[map[localScvfIdx][1]]} };
         }
         default:
-            DUNE_THROW(Dune::NotImplemented, "Box scvf geometries for dim=" << GridView::dimension
-                                                            << " dimWorld=" << GridView::dimensionworld
+            DUNE_THROW(Dune::NotImplemented, "Box scvf geometries for dim=" << dim
+                                                            << " dimWorld=" << dimWorld
                                                             << " corners=" << corners_);
         }
     }
 
     //! Create the sub control volume face geometries on the boundary
-    ScvfCornerStorage getBoundaryScvfCorners(const typename Intersection::Geometry& geometry,
+    ScvfCornerStorage getBoundaryScvfCorners(const Intersection& is,
+                                             const typename Intersection::Geometry& isGeom,
                                              unsigned int indexInIntersection) const
     {
+#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
+        const auto referenceElement = ReferenceElements::general(elementGeometry_.type());
+#else
+        const auto& referenceElement = ReferenceElements::general(elementGeometry_.type());
+#endif
+
+        const auto vIdxLocal = referenceElement.subEntity(is.indexInInside(), 1, indexInIntersection, dim);
         if (indexInIntersection == 0)
-            return ScvfCornerStorage({geometry.corner(0), geometry.center()});
+            return ScvfCornerStorage({p[vIdxLocal+1], isGeom.center()});
         else if (indexInIntersection == 1)
-            return ScvfCornerStorage({geometry.center(), geometry.corner(1)});
+            return ScvfCornerStorage({isGeom.center(), p[vIdxLocal+1]});
         else
             DUNE_THROW(Dune::InvalidStateException, "local index exceeds the number of corners of 2d intersections");
     }
@@ -355,9 +365,9 @@ class BoxGeometryHelper<GridView, 3, ScvType, ScvfType>
     static constexpr int maxPoints = 27;
 public:
     BoxGeometryHelper(const typename Element::Geometry& geometry)
-    : elementGeometry_(geometry), corners_(geometry.corners())
+    : elementGeometry_(geometry)
+    , corners_(geometry.corners())
     {
-        // extract the corners of the sub control volumes
 #if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
         const auto referenceElement = ReferenceElements::general(geometry.type());
 #else
@@ -503,29 +513,39 @@ public:
     }
 
     //! Create the sub control volume face geometries on the boundary
-    ScvfCornerStorage getBoundaryScvfCorners(const typename Intersection::Geometry& geometry,
+    ScvfCornerStorage getBoundaryScvfCorners(const Intersection& is,
+                                             const typename Intersection::Geometry& geometry,
                                              unsigned int indexInIntersection) const
     {
-        // extract the corners of the sub control volumes
+
 #if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
-        const auto referenceElement = FaceReferenceElements::general(geometry.type());
+        const auto referenceElement = ReferenceElements::general(elementGeometry_.type());
+        const auto faceRefElem = FaceReferenceElements::general(geometry.type());
 #else
-        const auto& referenceElement = FaceReferenceElements::general(geometry.type());
+        const auto& referenceElement = ReferenceElements::general(elementGeometry_.type());
+        const auto& faceRefElem = FaceReferenceElements::general(geometry.type());
 #endif
 
         GlobalPosition pi[9];
         auto corners = geometry.corners();
 
-        // the element center
+        // the facet center
         pi[0] = geometry.center();
 
-        // vertices
+        // corners
+        const auto idxInInside = is.indexInInside();
         for (int i = 0; i < corners; ++i)
-            pi[i+1] = geometry.corner(i);
+        {
+            const auto vIdxLocal = referenceElement.subEntity(idxInInside, 1, i, dim);
+            pi[i+1] = elementGeometry_.corner(vIdxLocal);
+        }
 
-        // face midpoints
-        for (int i = 0; i < referenceElement.size(1); ++i)
-            pi[i+corners+1] = geometry.global(referenceElement.position(i, 1));
+        // edge midpoints
+        for (int i = 0; i < faceRefElem.size(1); ++i)
+        {
+            const auto edgeIdxLocal = referenceElement.subEntity(idxInInside, 1, i, dim-1);
+            pi[i+corners+1] = p[edgeIdxLocal+corners_+1];
+        }
 
         // procees according to number of corners
         switch (corners)
@@ -533,13 +553,13 @@ public:
         case 3: // triangle
         {
             //! Only build the maps the first time we encounter a triangle
-            static const std::uint8_t vo = 1; //!< vertex offset in point vector p
-            static const std::uint8_t fo = 4; //!< face offset in point vector p
+            static const std::uint8_t vo = 1; //!< vertex offset in point vector pi
+            static const std::uint8_t eo = 4; //!< edge offset in point vector pi
             static const std::uint8_t map[3][4] =
             {
-                {vo+0, fo+0, fo+1, 0},
-                {vo+1, fo+2, fo+0, 0},
-                {vo+2, fo+1, fo+2, 0}
+                {vo+0, eo+0, eo+1, 0},
+                {vo+1, eo+2, eo+0, 0},
+                {vo+2, eo+1, eo+2, 0}
             };
 
             return ScvfCornerStorage{ {pi[map[indexInIntersection][0]],
@@ -550,14 +570,14 @@ public:
         case 4: // quadrilateral
         {
             //! Only build the maps the first time we encounter a quadrilateral
-            static const std::uint8_t vo = 1; //!< vertex offset in point vector p
-            static const std::uint8_t fo = 5; //!< face offset in point vector p
+            static const std::uint8_t vo = 1; //!< vertex offset in point vector pi
+            static const std::uint8_t eo = 5; //!< edge offset in point vector pi
             static const std::uint8_t map[4][4] =
             {
-                {vo+0, fo+2, fo+0, 0},
-                {vo+1, fo+1, fo+2, 0},
-                {vo+2, fo+0, fo+3, 0},
-                {vo+3, fo+3, fo+1, 0}
+                {vo+0, eo+2, eo+0, 0},
+                {vo+1, eo+1, eo+2, 0},
+                {vo+2, eo+0, eo+3, 0},
+                {vo+3, eo+3, eo+1, 0}
             };
 
             return ScvfCornerStorage{ {pi[map[indexInIntersection][0]],

dumux/discretization/box/subcontrolvolumeface.hh

@@ -138,7 +138,7 @@ public:
                             GridIndexType scvfIndex,
                             std::vector<LocalIndexType>&& scvIndices,
                             bool boundary = false)
-    : corners_(geometryHelper.getBoundaryScvfCorners(isGeometry, indexInIntersection)),
+    : corners_(geometryHelper.getBoundaryScvfCorners(intersection, isGeometry, indexInIntersection)),
       center_(0.0),
       unitOuterNormal_(intersection.centerUnitOuterNormal()),
       area_(geometryHelper.scvfArea(corners_)),

test/porousmediumflow/2p/implicit/incompressible/spatialparams.hh

@@ -20,8 +20,8 @@
  * \ingroup TwoPTests
  * \brief The spatial params for the incompressible 2p test
  */
-#ifndef DUMUX_COMPRESSIBLE_ONEP_TEST_SPATIAL_PARAMS_HH
-#define DUMUX_COMPRESSIBLE_ONEP_TEST_SPATIAL_PARAMS_HH
+#ifndef DUMUX_INCOMPRESSIBLE_TWOP_TEST_SPATIAL_PARAMS_HH
+#define DUMUX_INCOMPRESSIBLE_TWOP_TEST_SPATIAL_PARAMS_HH
 
 #include <dumux/discretization/methods.hh>