diff --git a/dumux/discretization/facecentered/diamond/fvelementgeometry.hh b/dumux/discretization/facecentered/diamond/fvelementgeometry.hh index 507c0b1e275905da110f8cb0d614666be9479559..03eed2cb39da551dc6e538bceaacf6ed91d0c477 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 43815b8853133220d9e9805efcd7d824609753ae..f992d839986cf7ad804ded8efc15c9efce612b07 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;