diff --git a/dumux/discretization/box/boxgeometryhelper.hh b/dumux/discretization/box/boxgeometryhelper.hh
index 9b525b5226ce7c9b873ad4237518a0356b38d1b8..f70785c255d14532ad98f27dfdf9915741a18e86 100644
--- a/dumux/discretization/box/boxgeometryhelper.hh
+++ b/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]],
diff --git a/dumux/discretization/box/subcontrolvolumeface.hh b/dumux/discretization/box/subcontrolvolumeface.hh
index 843ce184e0d3bee776b4bc3a14c25f68e683f714..f5709a5d395591ee124c0a776f5cfdb68d850fd8 100644
--- a/dumux/discretization/box/subcontrolvolumeface.hh
+++ b/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_)),
diff --git a/test/porousmediumflow/2p/implicit/incompressible/spatialparams.hh b/test/porousmediumflow/2p/implicit/incompressible/spatialparams.hh
index 043cdbc147eeda956af6d7729929ec59c6c449dd..0e5a73b95dfcbb77732b34f7f401fdcd1015479e 100644
--- a/test/porousmediumflow/2p/implicit/incompressible/spatialparams.hh
+++ b/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>