From 2899ede610362b50cd6e28d74c03cf6eed54e445 Mon Sep 17 00:00:00 2001
From: Timo Koch <timokoch@math.uio.no>
Date: Fri, 12 Aug 2022 09:39:17 +0000
Subject: [PATCH] [diamond] Fix and test fvGeometry.geometry(scv,scvf)
interface
---
.../facecentered/diamond/fvelementgeometry.hh | 29 +++++++++++++++++--
.../diamond/test_diamondgridgeometry.cc | 7 +++++
2 files changed, 33 insertions(+), 3 deletions(-)
diff --git a/dumux/discretization/facecentered/diamond/fvelementgeometry.hh b/dumux/discretization/facecentered/diamond/fvelementgeometry.hh
index 507c0b1e27..03eed2cb39 100644
--- a/dumux/discretization/facecentered/diamond/fvelementgeometry.hh
+++ b/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_;
diff --git a/test/discretization/facecentered/diamond/test_diamondgridgeometry.cc b/test/discretization/facecentered/diamond/test_diamondgridgeometry.cc
index 43815b8853..f992d83998 100644
--- a/test/discretization/facecentered/diamond/test_diamondgridgeometry.cc
+++ b/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;
--
GitLab