diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3aa4e623a626a72f7cb6f4ed7388628cbccbf0e0..25da3d04c98c34b547193ddad852ca41016e81ff 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -25,6 +25,7 @@ be evaluated with the function values provided in the same order as the names wh
- __Compositional Staggered__: The staggered discretization method has been updated to fix inconsistencies in handling density when evaluating diffusive boundary fluxes. In addition, a total mass balance is now always considered, rather than a total mole balance.
- __Parameters__: The template argument for `getParam` and `getParamFromGroup` now defaults to `std::string`, such that for string parameters one can simply write `const auto param = getParam("MyParam");`
- __Compositional Freeflow FCSFV__: Compositional transport has been implemented for the overhauled face-centered staggered finite volume scheme (FCSFV) discretization method. Ported tests are available in the testing suite.
+- __GridGeometry__: No longer store the corners/geometry for scv/scvfs also for cc-mpfa and staggered discretizations. This improves the memory footprint a lot. Geometries can be obtained via the local view of the grid geometry.
### Immediate interface changes not allowing/requiring a deprecation period:
@@ -41,6 +42,7 @@ be evaluated with the function values provided in the same order as the names wh
- `HAVE_PVPYTHON` to `DUMUX_HAVE_PVPYTHON`
### Deprecated properties/classes/functions/files, to be removed after 3.8:
+- __MPFA__: `scvf.corners()/facetCorner()/vertexCorner()` have been deprecated. Use functions from the local view of the grid geometry instead.
diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
index b4de9dd31bcb9266a688d47419e2743718b2d948..e1c7531bcdaebc36a86fb1bbdd0bf52071540b1c 100644
--- a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
@@ -19,6 +19,7 @@
#include <dune/common/exceptions.hh>
#include <dune/common/iteratorrange.hh>
+#include <dune/geometry/type.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/indextraits.hh>
@@ -26,6 +27,66 @@
namespace Dumux {
+#ifndef DOXYGEN
+namespace Detail::Mpfa {
+
+template<typename GridGeometry, typename SubControlVolumeFace>
+typename SubControlVolumeFace::Traits::Geometry makeScvfGeometry(const GridGeometry& gridGeometry,
+ const SubControlVolumeFace& scvf)
+{
+ static constexpr int dim = GridGeometry::GridView::dimension;
+
+ const auto& facetInfo = scvf.facetInfo();
+ const auto element = gridGeometry.element(facetInfo.elementIndex);
+ const auto elemGeo = element.geometry();
+ const auto refElement = referenceElement(elemGeo);
+ for (const auto& is : intersections(gridGeometry.gridView(), element))
+ {
+ if (is.indexInInside() == facetInfo.facetIndex)
+ {
+ const auto numCorners = is.geometry().corners();
+ const auto isPositions = GridGeometry::MpfaHelper::computeScvfCornersOnIntersection(
+ elemGeo, refElement, facetInfo.facetIndex, numCorners
+ );
+ return {
+ Dune::GeometryTypes::cube(dim-1),
+ GridGeometry::MpfaHelper::getScvfCorners(
+ isPositions, numCorners, facetInfo.facetCornerIndex
+ )
+ };
+ }
+ }
+ DUNE_THROW(Dune::InvalidStateException, "Could not construct scvf geometry");
+}
+
+template<typename GridGeometry, typename SubControlVolumeFace>
+auto getVertexCorner(const GridGeometry& gridGeometry, const SubControlVolumeFace& scvf)
+{
+ static constexpr int dim = GridGeometry::GridView::dimension;
+
+ const auto& facetInfo = scvf.facetInfo();
+ const auto element = gridGeometry.element(facetInfo.elementIndex);
+ const auto elemGeo = element.geometry();
+ const auto refElement = referenceElement(elemGeo);
+ return elemGeo.global(refElement.position(
+ refElement.subEntity(facetInfo.facetIndex, 1, facetInfo.facetCornerIndex, dim),
+ dim
+ ));
+}
+
+template<typename GridGeometry, typename SubControlVolumeFace>
+auto getFacetCorner(const GridGeometry& gridGeometry, const SubControlVolumeFace& scvf)
+{
+ const auto& facetInfo = scvf.facetInfo();
+ const auto element = gridGeometry.element(facetInfo.elementIndex);
+ const auto elemGeo = element.geometry();
+ const auto refElement = referenceElement(elemGeo);
+ return elemGeo.global(refElement.position(facetInfo.facetIndex, 1));
+}
+
+} // namespace Detail::Mpfa
+#endif // DOXYGEN
+
/*!
* \ingroup CCMpfaDiscretization
* \brief Stencil-local finite volume geometry (scvs and scvfs) for cell-centered mpfa models
@@ -52,6 +113,7 @@ class CCMpfaFVElementGeometry<GG, true>
using GridIndexType = typename IndexTraits<GridView>::GridIndex;
static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
public:
//! export type of the element
@@ -184,16 +246,17 @@ public:
typename Element::Geometry geometry(const SubControlVolume& scv) const
{ return gridGeometryPtr_->element(scv.dofIndex()).geometry(); }
- // suppress warnings due to current implementation
- // these interfaces should be used!
- #pragma GCC diagnostic push
- #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
-
//! Create the geometry of a given sub control volume face
typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const
- { return scvf.geometry(); }
+ { return Detail::Mpfa::makeScvfGeometry(gridGeometry(), scvf); }
+
+ //! Return the position of the scvf corner that coincides with an element vertex
+ typename SubControlVolumeFace::Traits::GlobalPosition vertexCorner(const SubControlVolumeFace& scvf) const
+ { return Detail::Mpfa::getVertexCorner(gridGeometry(), scvf); }
- #pragma GCC diagnostic pop
+ //! Return the corner of the scvf that is inside the facet the scvf is embedded in
+ typename SubControlVolumeFace::Traits::GlobalPosition facetCorner(const SubControlVolumeFace& scvf) const
+ { return Detail::Mpfa::getFacetCorner(gridGeometry(), scvf); }
private:
@@ -215,8 +278,8 @@ class CCMpfaFVElementGeometry<GG, false>
using GridIndexType = typename IndexTraits<GridView>::GridIndex;
using MpfaHelper = typename GG::MpfaHelper;
- static const int dim = GridView::dimension;
- static const int dimWorld = GridView::dimensionworld;
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
using CoordScalar = typename GridView::ctype;
public:
@@ -345,20 +408,21 @@ public:
bool hasBoundaryScvf() const
{ return hasBoundaryScvf_; }
- // suppress warnings due to current implementation
- // these interfaces should be used!
- #pragma GCC diagnostic push
- #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
-
//! Create the geometry of a given sub control volume
typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const
- { return scv.geometry(); }
+ { return gridGeometryPtr_->element(scv.dofIndex()).geometry(); }
//! Create the geometry of a given sub control volume face
typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const
- { return scvf.geometry(); }
+ { return Detail::Mpfa::makeScvfGeometry(gridGeometry(), scvf); }
- #pragma GCC diagnostic pop
+ //! Return the position of the scvf corner that coincides with an element vertex
+ typename SubControlVolumeFace::Traits::GlobalPosition vertexCorner(const SubControlVolumeFace& scvf) const
+ { return Detail::Mpfa::getVertexCorner(gridGeometry(), scvf); }
+
+ //! Return the corner of the scvf that is inside the facet the scvf is embedded in
+ typename SubControlVolumeFace::Traits::GlobalPosition facetCorner(const SubControlVolumeFace& scvf) const
+ { return Detail::Mpfa::getFacetCorner(gridGeometry(), scvf); }
private:
@@ -489,10 +553,15 @@ private:
continue;
hasBoundaryScvf_ = (hasBoundaryScvf_ || is.boundary());
-
+ typename SubControlVolumeFace::FacetInfo facetInfo{
+ gridGeometry().elementMapper().index(e),
+ useNeighbor ? is.indexInOutside() : is.indexInInside(),
+ c
+ };
scvfs_.emplace_back(MpfaHelper(),
MpfaHelper::getScvfCorners(isPositions, numCorners, c),
is,
+ std::move(facetInfo),
vIdxGlobal,
vIdxLocal,
scvFaceIndices[scvfCounter],
@@ -580,9 +649,15 @@ private:
}
// build scvf
+ typename SubControlVolumeFace::FacetInfo facetInfo{
+ gridGeometry().elementMapper().index(e),
+ useNeighbor ? is.indexInOutside() : is.indexInInside(),
+ c
+ };
neighborScvfs_.emplace_back(MpfaHelper(),
MpfaHelper::getScvfCorners(isPositions, numCorners, c),
is,
+ std::move(facetInfo),
vIdxGlobal,
vIdxLocal,
scvFaceIndices[scvfCounter],
diff --git a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
index f89340ab07379df227f2ad2355b35ad62cca05e2..5ddc0345db8ce1887ecc71ae5133b14760942a8e 100644
--- a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
@@ -308,7 +308,7 @@ private:
const bool usesSecondaryIV = secondaryIvIndicator_(element, is, isBranchingPoint);
// make the scv faces belonging to each corner of the intersection
- for (std::size_t c = 0; c < numCorners; ++c)
+ for (int c = 0; c < numCorners; ++c)
{
// get the global vertex index the scv face is connected to
const auto vIdxLocal = refElement.subEntity(indexInElement, 1, c, dim);
@@ -337,9 +337,15 @@ private:
} ();
scvfIndexSet.push_back(scvfIdx);
+ typename SubControlVolumeFace::FacetInfo facetInfo{
+ this->elementMapper().index(e),
+ useNeighbor ? is.indexInOutside() : is.indexInInside(),
+ c
+ };
scvfs_.emplace_back(MpfaHelper(),
MpfaHelper::getScvfCorners(isPositions, numCorners, c),
is,
+ std::move(facetInfo),
vIdxGlobal,
vIdxLocal,
scvfIdx,
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/scvgeometryhelper.hh b/dumux/discretization/cellcentered/mpfa/omethod/scvgeometryhelper.hh
index 4d772a3982d24288adfc4f208fce38675953f4c0..4db457755c16f2df866bade3aec9e8461c734867 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/scvgeometryhelper.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/scvgeometryhelper.hh
@@ -71,9 +71,9 @@ public:
typename MLGTraits::template CornerStorage<dim, dimWorld>::Type corners;
corners[0] = fvGeometry.scv( scv.gridScvIndex() ).center();
- corners[1] = firstGridScvf.facetCorner();
- corners[2] = secondGridScvf.facetCorner();
- corners[3] = secondGridScvf.vertexCorner();
+ corners[1] = fvGeometry.facetCorner(firstGridScvf);
+ corners[2] = fvGeometry.facetCorner(secondGridScvf);
+ corners[3] = fvGeometry.vertexCorner(secondGridScvf);
using std::swap;
typename LocalScvType::LocalBasis basis;
diff --git a/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh
index d677c55da6d041b04515784f3eb6aeb4c1887dc1..7249250cadb874fbc9dc5f62fadbab7d90c5aa51 100644
--- a/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh
+++ b/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh
@@ -12,6 +12,7 @@
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_SUBCONTROLVOLUMEFACE_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_SUBCONTROLVOLUMEFACE_HH
+#include <utility>
#include <vector>
#include <array>
@@ -92,6 +93,14 @@ class CCMpfaSubControlVolumeFace
using BoundaryFlag = typename T::BoundaryFlag;
public:
+ // Information on the intersection from which this scvf was constructed
+ struct FacetInfo
+ {
+ GridIndexType elementIndex;
+ int facetIndex;
+ int facetCornerIndex;
+ };
+
//! export the type used for global coordinates
using GlobalPosition = typename T::GlobalPosition;
//! state the traits public and thus export all types
@@ -103,6 +112,7 @@ public:
* \param helper The helper class for mpfa schemes
* \param corners The corners of the scv face
* \param intersection The intersection
+ * \param facetInfo Information on the facet from which this scvf is constructed
* \param vIdxGlobal The global vertex index the scvf is connected to
* \param vIdxLocal The element-local vertex index the scvf is connected to
* \param scvfIndex The global index of this scv face
@@ -116,6 +126,7 @@ public:
CCMpfaSubControlVolumeFace(const MpfaHelper& helper,
CornerStorage&& corners,
const Intersection& intersection,
+ FacetInfo facetInfo,
GridIndexType vIdxGlobal,
unsigned int vIdxLocal,
GridIndexType scvfIndex,
@@ -133,6 +144,7 @@ public:
, center_(0.0)
, unitOuterNormal_(intersection.centerUnitOuterNormal())
, boundaryFlag_{intersection}
+ , facetInfo_{std::move(facetInfo)}
{
// compute the center of the scvf
for (const auto& corner : corners_)
@@ -182,14 +194,17 @@ public:
{ return outsideScvIndices_; }
//! Returns the number of corners
+ [[deprecated("Will be removed after 3.8. Use fvGeometry.geometry(scvf).corners().")]]
std::size_t corners() const
{ return corners_.size(); }
//! Returns the global position of the vertex the scvf is connected to
+ [[deprecated("Will be removed after 3.8. Use fvGeometry.vertexCorner(scvf)")]]
const GlobalPosition& vertexCorner() const
{ return corners_.back(); }
//! Returns the global position of the center of the element facet this scvf is embedded in
+ [[deprecated("Will be removed after 3.8. Use fvGeometry.facetCorner(scvf)")]]
const GlobalPosition& facetCorner() const
{ return corners_[0]; }
@@ -205,15 +220,14 @@ public:
const GlobalPosition& unitOuterNormal() const
{ return unitOuterNormal_; }
- //! The geometry of the sub control volume face
- [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scvf).")]]
- Geometry geometry() const
- { return Geometry(Dune::GeometryTypes::cube(Geometry::mydimension), corners_); }
-
//! Return the boundary flag
typename BoundaryFlag::value_type boundaryFlag() const
{ return boundaryFlag_.get(); }
+ //! Return information on the facet from which this scvf was constructed
+ const FacetInfo& facetInfo() const
+ { return facetInfo_; }
+
private:
bool boundary_;
GridIndexType vertexIndex_;
@@ -228,6 +242,7 @@ private:
GlobalPosition unitOuterNormal_;
Scalar area_;
BoundaryFlag boundaryFlag_;
+ FacetInfo facetInfo_;
};
} // end namespace Dumux
diff --git a/dumux/discretization/facecentered/staggered/fvelementgeometry.hh b/dumux/discretization/facecentered/staggered/fvelementgeometry.hh
index 5d4bc1024f7b6362b2058a816a614a3372d3d0eb..297df24ade3f394cd87d9b9c04fbb38f92d55bdb 100644
--- a/dumux/discretization/facecentered/staggered/fvelementgeometry.hh
+++ b/dumux/discretization/facecentered/staggered/fvelementgeometry.hh
@@ -13,6 +13,7 @@
#define DUMUX_DISCRETIZATION_FACECENTERED_STAGGERED_FV_ELEMENT_GEOMETRY_HH
#include <utility>
+#include <bitset>
#include <dune/common/rangeutilities.hh>
#include <dune/common/reservedvector.hh>
@@ -26,6 +27,73 @@
namespace Dumux {
+#ifndef DOXYGEN
+namespace Detail::FCStaggered {
+
+template<class FVElementGeometry, class SubControlVolume>
+typename SubControlVolume::Traits::Geometry scvGeometry(const FVElementGeometry& fvGeometry,
+ const SubControlVolume& scv)
+{
+ typename SubControlVolume::Traits::CornerStorage corners{};
+ // select the containing element
+ const auto elementGeometry = (scv.elementIndex() != fvGeometry.elementIndex()) ?
+ fvGeometry.element().geometry() :
+ fvGeometry.gridGeometry().element(scv.elementIndex()).geometry();
+
+ const auto center = elementGeometry.center();
+ const auto dofAxis = scv.dofAxis();
+ for (int i = 0; i < corners.size(); ++i)
+ {
+ auto& corner = corners[i];
+
+ // copy the corner of the corresponding element
+ corner = elementGeometry.corner(i);
+
+ // shift the corner such that the scv covers half of the element
+ // (keep the corner positions at the face with the staggered dof)
+ if ((corner[dofAxis] - center[dofAxis]) * scv.directionSign() < 0.0)
+ corner[dofAxis] = center[dofAxis];
+ }
+
+ return {corners.front(), corners.back()};
+}
+
+template<class FVElementGeometry, class SubControlVolumeFace>
+typename SubControlVolumeFace::Traits::Geometry scvfGeometry(const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf)
+{
+ const auto normalAxis = scvf.normalAxis();
+ const auto center = scvf.center();
+ const auto shift = scvf.ipGlobal() - center;
+ const auto dofAxis = scvf.isLateral() ? Dumux::normalAxis(shift) : normalAxis;
+ const auto insideElementIndex = (fvGeometry.scv(scvf.insideScvIdx())).elementIndex();
+ const auto elementGeometry = (insideElementIndex != fvGeometry.elementIndex()) ?
+ fvGeometry.element().geometry() :
+ fvGeometry.gridGeometry().element(insideElementIndex).geometry();
+
+ auto corners = std::array{
+ elementGeometry.corner(0),
+ elementGeometry.corner(elementGeometry.corners() - 1)
+ };
+
+ // shift corners to scvf plane and halve lateral faces
+ for (int i = 0; i < corners.size(); ++i)
+ {
+ auto& corner = corners[i];
+ corner[normalAxis] = center[normalAxis];
+ if (scvf.isLateral() && (corner - center)*shift < 0.0)
+ corner[dofAxis] = elementGeometry.center()[dofAxis];
+ }
+
+ auto inPlaneAxes = std::move(std::bitset<SubControlVolumeFace::Traits::dimWorld>{}.set());
+ inPlaneAxes.set(normalAxis, false);
+
+ return {corners[0], corners[1], inPlaneAxes};
+}
+
+} // end namespace Detail::FCStaggered
+#endif // DOXYGEN
+
template<class GG, bool cachingEnabled>
class FaceCenteredStaggeredFVElementGeometry;
@@ -209,20 +277,17 @@ public:
DUNE_THROW(Dune::InvalidStateException, "No outside scvf found");
}
- // suppress warnings due to current implementation
- // these interfaces should be used!
- #pragma GCC diagnostic push
- #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
-
//! Create the geometry of a given sub control volume
typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const
- { return scv.geometry(); }
+ {
+ return Detail::FCStaggered::scvGeometry(*this, scv);
+ }
//! Create the geometry of a given sub control volume face
typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const
- { return scvf.geometry(); }
-
- #pragma GCC diagnostic pop
+ {
+ return Detail::FCStaggered::scvfGeometry(*this, scvf);
+ }
private:
@@ -436,20 +501,17 @@ public:
DUNE_THROW(Dune::InvalidStateException, "No outside scvf found");
}
- // suppress warnings due to current implementation
- // these interfaces should be used!
- #pragma GCC diagnostic push
- #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
-
//! Create the geometry of a given sub control volume
typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const
- { return scv.geometry(); }
+ {
+ return Detail::FCStaggered::scvGeometry(*this, scv);
+ }
//! Create the geometry of a given sub control volume face
typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const
- { return scvf.geometry(); }
-
- #pragma GCC diagnostic pop
+ {
+ return Detail::FCStaggered::scvfGeometry(*this, scvf);
+ }
private:
//! Binding of an element preparing the geometries of the whole stencil
diff --git a/dumux/discretization/facecentered/staggered/subcontrolvolume.hh b/dumux/discretization/facecentered/staggered/subcontrolvolume.hh
index a33bc9656bf6183bce4cbb277939f132c086f0f6..9dab0bedfcfff1f663690cd4d78ceeb6d153ef2c 100644
--- a/dumux/discretization/facecentered/staggered/subcontrolvolume.hh
+++ b/dumux/discretization/facecentered/staggered/subcontrolvolume.hh
@@ -86,21 +86,6 @@ public:
, eIdx_(eIdx)
, boundary_(boundary)
{
- for (int i = 0; i < corners_.size(); ++i)
- {
- auto& corner = corners_[i];
-
- // copy the corner of the corresponding element
- corner = elementGeometry.corner(i);
-
- // shift the corner such that the scvf covers half of the lateral facet
- // (keep the outer corner positions)
- using std::abs;
- const auto eps = 1e-8; // TODO
- if (abs(corner[dofAxis] - intersectionGeometry.center()[dofAxis]) > eps)
- corner[dofAxis] = elementGeometry.center()[dofAxis];
- }
-
directionSign_ = (indexInElement % 2) ? 1.0 : -1.0;
}
@@ -139,17 +124,9 @@ public:
bool boundary() const
{ return boundary_; }
- //! The geometry of the sub control volume face
- [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scv).")]]
- Geometry geometry() const
- {
- return Geometry(corners_.front(), corners_.back());
- }
-
private:
GlobalPosition center_;
GlobalPosition dofPosition_;
- CornerStorage corners_;
Scalar volume_;
GridIndexType globalIndex_;
SmallLocalIndexType indexInElement_;
diff --git a/dumux/discretization/facecentered/staggered/subcontrolvolumeface.hh b/dumux/discretization/facecentered/staggered/subcontrolvolumeface.hh
index 9adb53bb86e59e2504a2b38880d81af5e899c237..bc7a5a8ec25dd9bc1e83b802d7d3cd70602bb5f9 100644
--- a/dumux/discretization/facecentered/staggered/subcontrolvolumeface.hh
+++ b/dumux/discretization/facecentered/staggered/subcontrolvolumeface.hh
@@ -13,7 +13,6 @@
#define DUMUX_DISCRETIZATION_FACECENTERED_STAGGERED_SUBCONTROLVOLUMEFACE_HH
#include <array>
-#include <utility>
#include <dune/common/fvector.hh>
#include <dune/geometry/type.hh>
@@ -98,13 +97,6 @@ public:
if (!boundary())
outerNormalSign_ *= -1.0;
-
- // the corners (coincide with intersection corners for boundary scvfs)
- const auto frontalOffSet = boundary() ? GlobalPosition(0.0)
- : intersectionGeometry.center() - elementGeometry.center();
-
- for (int i = 0; i < corners_.size(); ++i)
- corners_[i] = intersectionGeometry.corner(i) + frontalOffSet;
}
//! The constructor for lateral faces
@@ -131,22 +123,6 @@ public:
const auto shift = intersectionGeometry.center() - elementGeometry.center();
ipGlobal_ = lateralFacetGeometry.center() + shift;
center_ = 0.5*(lateralFacetGeometry.center() + ipGlobal_);
-
- const auto dofAxis = Dumux::normalAxis(shift);
- const auto eps = shift.two_norm() * 1e-8;
-
- for (int i = 0; i < corners_.size(); ++i)
- {
- // copy the corner of the corresponding lateral facet
- auto& corner = corners_[i];
- corner = lateralFacetGeometry.corner(i);
-
- // shift the corner such that the scvf covers half of the lateral facet
- // (keep the outer corner positions)
- using std::abs;
- if (abs(corner[dofAxis] - intersectionGeometry.center()[dofAxis]) > eps)
- corner[dofAxis] = elementGeometry.center()[dofAxis];
- }
}
//! The center of the sub control volume face
@@ -203,19 +179,9 @@ public:
std::int_least8_t directionSign() const
{ return outerNormalSign_; }
- //! The geometry of the sub control volume face
- [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scvf).")]]
- Geometry geometry() const
- {
- auto inPlaneAxes = std::move(std::bitset<T::dimWorld>{}.set());
- inPlaneAxes.set(normalAxis_, false);
- return { corners_.front(), corners_.back(), inPlaneAxes };
- }
-
private:
GlobalPosition center_;
GlobalPosition ipGlobal_;
- CornerStorage corners_;
std::array<GridIndexType, 2> globalScvIndices_;
SmallLocalIndexType localScvfIdx_;
GridIndexType globalScvfIdx_;
diff --git a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
index 5caf10a4e67d368336885a8be6f793abcd08d954..219fcc9f7e7c0692598e72df92c01a8720b89862 100644
--- a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
@@ -19,7 +19,6 @@
#include <dune/geometry/multilineargeometry.hh>
#include <dumux/common/indextraits.hh>
-#include <dumux/common/typetraits/isvalid.hh>
#include <dumux/discretization/subcontrolvolumefacebase.hh>
#include <dumux/discretization/staggered/subcontrolvolumeface.hh>
#include <dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh>
@@ -28,22 +27,6 @@
namespace Dumux {
-#ifndef DOXYGEN
-namespace Detail {
-// helper struct detecting if the container class storing the scvf's corners has a resize function
-struct HasResize
-{
- template<class Container> auto operator()(Container&& c)
- -> decltype(c.resize(1)) {}
-};
-
-template<class C>
-constexpr bool hasResize()
-{ return decltype(isValid(HasResize{})(std::declval<C>()))::value; }
-
-} // end namespace Detail
-#endif
-
/*!
* \ingroup StaggeredDiscretization
* \brief Default traits class to be used for the sub-control volume faces
@@ -65,22 +48,9 @@ struct FreeFlowStaggeredDefaultScvfGeometryTraits
static constexpr int dim = Grid::dimension;
static constexpr int dimWorld = Grid::dimensionworld;
- // we use geometry traits that use static corner vectors to and a fixed geometry type
- template <class ct>
- struct ScvfMLGTraits : public Dune::MultiLinearGeometryTraits<ct>
- {
- // we use static vectors to store the corners as we know
- // the number of corners in advance (2^(dim-1) corners (1<<(dim-1))
- template< int mydim, int cdim >
- struct CornerStorage
- {
- using Type = std::array< Dune::FieldVector< ct, cdim >, (1<<(dim-1)) >;
- };
- };
-
- using Geometry = Dune::MultiLinearGeometry<Scalar, dim-1, dimWorld, ScvfMLGTraits<Scalar> >;
- using CornerStorage = typename ScvfMLGTraits<Scalar>::template CornerStorage<dim-1, dimWorld>::Type;
- using GlobalPosition = typename CornerStorage::value_type;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using Geometry = Dune::AxisAlignedCubeGeometry<Scalar, dim-1, dimWorld>;
};
/*!
@@ -115,7 +85,6 @@ class FreeFlowStaggeredSubControlVolumeFace
using Geometry = typename T::Geometry;
using GridIndexType = typename IndexTraits<GV>::GridIndex;
using LocalIndexType = typename IndexTraits<GV>::LocalIndex;
- using CornerStorage = typename T::CornerStorage;
using PairData = typename T::PairData;
using AxisData = typename T::AxisData;
@@ -146,7 +115,6 @@ public:
const std::vector<GridIndexType>& scvIndices,
const typename T::GeometryHelper& geometryHelper)
: ParentType(),
- geomType_(isGeometry.type()),
area_(isGeometry.volume()),
center_(isGeometry.center()),
unitOuterNormal_(is.centerUnitOuterNormal()),
@@ -160,11 +128,9 @@ public:
outerNormalSign_(sign(unitOuterNormal_[directionIndex()])),
isGhostFace_(false)
{
- if constexpr (Detail::hasResize<CornerStorage>())
- corners_.resize(isGeometry.corners());
-
- for (int i = 0; i < isGeometry.corners(); ++i)
- corners_[i] = isGeometry.corner(i);
+ dimensions[0] = (isGeometry.corner(1) - isGeometry.corner(0)).two_norm();
+ if constexpr (dim == 3)
+ dimensions[1] = (isGeometry.corner(2) - isGeometry.corner(0)).two_norm();
}
//! The center of the sub control volume face
@@ -222,13 +188,6 @@ public:
return scvfIndex_;
}
- //! The geometry of the sub control volume face
- [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scvf).")]]
- const Geometry geometry() const
- {
- return Geometry(geomType_, corners_);
- }
-
//! The local index of this sub control volume face
LocalIndexType localFaceIdx() const
{
@@ -280,15 +239,10 @@ public:
//! Returns the length of the face in a certain direction (adaptation of area() for 3d)
Scalar faceLength(const int localSubFaceIdx) const
{
- if (dim == 3)
- {
- if (localSubFaceIdx < 2)
- return (corners_[1] - corners_[0]).two_norm();
- else
- return (corners_[2] - corners_[0]).two_norm();
- }
+ if (dim == 3 && localSubFaceIdx > 1)
+ return dimensions[1];
else
- return (corners_[1] - corners_[0]).two_norm();
+ return dimensions[0];
}
/*!
@@ -442,8 +396,7 @@ public:
{ isGhostFace_ = isGhostFaceFlag; }
private:
- Dune::GeometryType geomType_;
- CornerStorage corners_;
+ std::array<Scalar, dim-1> dimensions;
Scalar area_;
GlobalPosition center_;
GlobalPosition unitOuterNormal_;
diff --git a/dumux/discretization/staggered/fvelementgeometry.hh b/dumux/discretization/staggered/fvelementgeometry.hh
index ee9ec9747ea74296c8a387ec142085998a89acff..8279145c0944034c3c14a5dc7530e914c200c0cc 100644
--- a/dumux/discretization/staggered/fvelementgeometry.hh
+++ b/dumux/discretization/staggered/fvelementgeometry.hh
@@ -15,6 +15,7 @@
#include <optional>
#include <dumux/common/indextraits.hh>
#include <dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh>
+#include <dumux/discretization/facecentered/staggered/normalaxis.hh>
namespace Dumux {
@@ -207,20 +208,32 @@ public:
bool hasBoundaryScvf() const
{ return hasBoundaryScvf_; }
- // suppress warnings due to current implementation
- // these interfaces should be used!
- #pragma GCC diagnostic push
- #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
-
//! Create the geometry of a given sub control volume
typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const
- { return scv.geometry(); }
+ { return gridGeometryPtr_->element(scv.dofIndex()).geometry(); }
//! Create the geometry of a given sub control volume face
typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const
- { return scvf.geometry(); }
+ {
+ const auto insideElementIndex = scvf.insideScvIdx();
+ const auto elementGeometry = (insideElementIndex != scvIndices_[0]) ?
+ element_->geometry() :
+ gridGeometryPtr_->element(insideElementIndex).geometry();
+ const auto center = elementGeometry.center();
+ const auto normalAxis = Dumux::normalAxis(scvf.unitOuterNormal());
+
+ auto lowerLeft = elementGeometry.corner(0);
+ auto upperRight = elementGeometry.corner(elementGeometry.corners()-1);
- #pragma GCC diagnostic pop
+ // shift corners to scvf plane and halve lateral faces
+ lowerLeft[normalAxis] = center[normalAxis];
+ upperRight[normalAxis] = center[normalAxis];
+
+ auto inPlaneAxes = std::move(std::bitset<SubControlVolumeFace::Traits::dimWorld>{}.set());
+ inPlaneAxes.set(normalAxis, false);
+
+ return {lowerLeft, upperRight, inPlaneAxes};
+ }
private:
//! Binding of an element preparing the geometries only inside the element
diff --git a/dumux/discretization/staggered/subcontrolvolumeface.hh b/dumux/discretization/staggered/subcontrolvolumeface.hh
index 6adf5ef48044914edb3a56eea382d3f08ec799e9..1a6b43d1f0a9289acbcd8709f9f3ffc2dadf3cad 100644
--- a/dumux/discretization/staggered/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/subcontrolvolumeface.hh
@@ -85,11 +85,15 @@ private:
template<class GridView>
struct StaggeredDefaultScvfGeometryTraits
{
- using Geometry = typename GridView::template Codim<1>::Geometry;
using GridIndexType = typename IndexTraits<GridView>::GridIndex;
using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
using Scalar = typename GridView::ctype;
- using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimensionworld>;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ static constexpr int dim = GridView::Grid::dimension;
+ static constexpr int dimWorld = GridView::Grid::dimensionworld;
+ using Geometry = Dune::AxisAlignedCubeGeometry<Scalar, dim-1, dimWorld>;
};
/*!
@@ -129,7 +133,6 @@ public:
const std::vector<GridIndexType>& scvIndices,
const GeometryHelper& geometryHelper)
: ParentType()
- , geomType_(isGeometry.type())
, area_(isGeometry.volume())
, center_(isGeometry.center())
, unitOuterNormal_(is.centerUnitOuterNormal())
@@ -137,10 +140,6 @@ public:
, scvIndices_(scvIndices)
, boundary_(is.boundary())
{
- corners_.resize(isGeometry.corners());
- for (int i = 0; i < isGeometry.corners(); ++i)
- corners_[i] = isGeometry.corner(i);
-
dofIdx_ = geometryHelper.dofIndex();
localFaceIdx_ = geometryHelper.localFaceIndex();
}
@@ -200,13 +199,6 @@ public:
return scvfIndex_;
}
- //! The geometry of the sub control volume face
- [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scvf).")]]
- const Geometry geometry() const
- {
- return Geometry(geomType_, corners_);
- }
-
//! The global index of the dof living on this face
GridIndexType dofIndex() const
{
@@ -220,8 +212,6 @@ public:
}
private:
- Dune::GeometryType geomType_;
- std::vector<GlobalPosition> corners_;
Scalar area_;
GlobalPosition center_;
GlobalPosition unitOuterNormal_;
diff --git a/dumux/multidomain/facet/cellcentered/mpfa/couplingmapper.hh b/dumux/multidomain/facet/cellcentered/mpfa/couplingmapper.hh
index 9c9d612d1d8d7c129572e0c0ab5d189429b96be2..54fc7ca8de8c1f034105b3e6784e69302b768e2b 100644
--- a/dumux/multidomain/facet/cellcentered/mpfa/couplingmapper.hh
+++ b/dumux/multidomain/facet/cellcentered/mpfa/couplingmapper.hh
@@ -127,7 +127,7 @@ public:
else
{
const auto eps = lowDimGeometry.volume()*1e-8;
- const auto diffVec = lowDimGeometry.center()-scvf.facetCorner();
+ const auto diffVec = lowDimGeometry.center()-fvGeometry.facetCorner(scvf);
if ( Dune::FloatCmp::eq<GlobalPosition, Dune::FloatCmp::CmpStyle::absolute>(diffVec, GlobalPosition(0.0), eps) )
embeddedScvfIndices.push_back(scvf.index());
diff --git a/test/discretization/cellcentered/CMakeLists.txt b/test/discretization/cellcentered/CMakeLists.txt
index 6e2ca56bcc4d2c3b41ebf2b4c1dd7e3aa68ec146..3e84f26153d4fbe8c2dfd6ad914fe975c85baf76 100644
--- a/test/discretization/cellcentered/CMakeLists.txt
+++ b/test/discretization/cellcentered/CMakeLists.txt
@@ -1,4 +1,5 @@
# SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
# SPDX-License-Identifier: GPL-3.0-or-later
-add_subdirectory("tpfa")
+add_subdirectory(tpfa)
+add_subdirectory(mpfa)
diff --git a/test/discretization/cellcentered/mpfa/CMakeLists.txt b/test/discretization/cellcentered/mpfa/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..606c9160b7d29553f1074f2ae16f0ee262aeb9dc
--- /dev/null
+++ b/test/discretization/cellcentered/mpfa/CMakeLists.txt
@@ -0,0 +1,11 @@
+# SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+dumux_add_test(NAME test_mpfafvgeometry
+ SOURCES test_mpfafvgeometry.cc
+ LABELS unit discretization)
+
+dumux_add_test(NAME test_mpfafvgeometry_caching
+ SOURCES test_mpfafvgeometry.cc
+ LABELS unit discretization
+ COMPILE_DEFINITIONS ENABLECACHING=1)
diff --git a/test/discretization/cellcentered/mpfa/test_mpfafvgeometry.cc b/test/discretization/cellcentered/mpfa/test_mpfafvgeometry.cc
new file mode 100644
index 0000000000000000000000000000000000000000..447e260792a1e9ccbf0113b33e83d91b93a270e0
--- /dev/null
+++ b/test/discretization/cellcentered/mpfa/test_mpfafvgeometry.cc
@@ -0,0 +1,159 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ * \brief Test for the mpfa finite volume element geometry.
+ */
+#include <config.h>
+
+#include <iostream>
+#include <vector>
+#include <utility>
+#include <algorithm>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/float_cmp.hh>
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/common/initialize.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/discretization/ccmpfa.hh>
+
+#ifndef ENABLECACHING
+#define ENABLECACHING 0
+#endif
+
+namespace Dumux::Properties {
+namespace TTag { struct TestTag { using InheritsFrom = std::tuple<CCMpfaModel>; }; }
+
+template<typename TypeTag>
+struct Grid<TypeTag, TTag::TestTag> { using type = Dune::YaspGrid<2>; };
+
+template<typename TypeTag>
+struct Scalar<TypeTag, TTag::TestTag> { using type = double; };
+
+template<typename TypeTag>
+struct EnableGridGeometryCache<TypeTag, TTag::TestTag> {
+ static constexpr bool value = ENABLECACHING;
+};
+
+} // namespace Dumux::Properties
+
+template<class Range>
+auto size(const Range& r)
+{ return std::distance(std::begin(r), std::end(r)); }
+
+template<class FVElementGeometry>
+auto numBoundaryScvfs(const FVElementGeometry& fvGeometry)
+{
+ unsigned int count = 0;
+ for (const auto& scvf : scvfs(fvGeometry))
+ if (scvf.boundary())
+ count++;
+ return count;
+}
+
+template<typename Geometry>
+auto popped(std::vector<typename Geometry::GlobalCoordinate> expected, const Geometry& geo)
+{
+ using T = std::vector<typename Geometry::GlobalCoordinate>;
+ for (int c = 0; c < geo.corners(); ++c)
+ {
+ auto it = std::find_if(
+ expected.begin(),
+ expected.end(),
+ [corner=geo.corner(c)] (const auto& exp) {
+ return Dune::FloatCmp::eq(corner, exp);
+ });
+ if (it == expected.end())
+ {
+ std::cout << "Could not find corner " << geo.corner(c) << std::endl;
+ return std::optional<T>{};
+ }
+ expected.erase(it);
+ }
+ return std::optional<T>{std::move(expected)};
+}
+
+int main (int argc, char *argv[])
+{
+ using namespace Dumux;
+
+ initialize(argc, argv);
+ Parameters::init(argc, argv);
+
+ std::cout << "Checking the FVGeometries, SCVs and SCV faces" << std::endl;
+
+ using TypeTag = Properties::TTag::TestTag;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+
+ GetPropType<TypeTag, Properties::Grid> grid{{1.0, 1.0}, {1, 1}};
+ auto leafGridView = grid.leafGridView();
+ GridGeometry gridGeometry(leafGridView);
+
+ for (const auto& element : elements(leafGridView))
+ {
+ auto fvGeometry = localView(gridGeometry);
+
+ if (fvGeometry.isBound())
+ DUNE_THROW(Dune::Exception, "Local view should not be bound at this point");
+
+ fvGeometry.bind(element);
+ if (!fvGeometry.isBound())
+ DUNE_THROW(Dune::Exception, "Local view should be bound at this point");
+
+ const auto eIdx = gridGeometry.elementMapper().index(element);
+ const auto eIdxBound = gridGeometry.elementMapper().index(fvGeometry.element());
+ if (eIdx != eIdxBound)
+ DUNE_THROW(Dune::Exception, "Bound element index does not match");
+
+ if (size(scvs(fvGeometry)) != 1) DUNE_THROW(Dune::InvalidStateException, "Unexpected number of scvs");
+ if (size(scvfs(fvGeometry)) != 8) DUNE_THROW(Dune::InvalidStateException, "Unexpected number of scvfs");
+ if (numBoundaryScvfs(fvGeometry) != 8) DUNE_THROW(Dune::InvalidStateException, "Unexpected number of boundary scvfs");
+
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ const auto geo = fvGeometry.geometry(scv);
+ if (!Dune::FloatCmp::eq(geo.volume(), 1.0))
+ DUNE_THROW(Dune::InvalidStateException, "Unexpected scv volume");
+ const auto result = popped({
+ {0.0, 0.0},
+ {1.0, 0.0},
+ {1.0, 1.0},
+ {0.0, 1.0}
+ }, geo);
+ if (!result.has_value() || !result.value().empty())
+ DUNE_THROW(Dune::InvalidStateException, "Unexpected scv corners");
+ }
+
+ std::vector<typename SubControlVolumeFace::GlobalPosition> expected{
+ {0.0, 0.0}, {0.5, 0.0},
+ {0.5, 0.0}, {1.0, 0.0},
+ {1.0, 0.0}, {1.0, 0.5},
+ {1.0, 0.5}, {1.0, 1.0},
+ {1.0, 1.0}, {0.5, 1.0},
+ {0.5, 1.0}, {0.0, 1.0},
+ {0.0, 1.0}, {0.0, 0.5},
+ {0.0, 0.5}, {0.0, 0.0}
+ };
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+ const auto geo = fvGeometry.geometry(scvf);
+ if (!Dune::FloatCmp::eq(geo.volume(), 0.5))
+ DUNE_THROW(Dune::InvalidStateException, "Unexpected scvf area " << geo.volume());
+ const auto shrunk = popped(expected, geo);
+ if (!shrunk.has_value())
+ DUNE_THROW(Dune::InvalidStateException, "Unexpected scvf corners");
+ expected = shrunk.value();
+ }
+ if (!expected.empty())
+ DUNE_THROW(Dune::InvalidStateException, "Could not find all expected scvf corners");
+ }
+}