diff --git a/dumux/common/pointsource.hh b/dumux/common/pointsource.hh
index eaf3b3f6bd901c2b07a52964bdb3fcf9b63ed353..2ca28690edaf53502abd10fbb1fbd0abb7f84ffe 100644
--- a/dumux/common/pointsource.hh
+++ b/dumux/common/pointsource.hh
@@ -316,7 +316,7 @@ public:
constexpr int dim = GridGeometry::GridView::dimension;
Dune::ReservedVector<std::size_t, 1<<dim> scvIndices;
for (const auto& scv : scvs(fvGeometry))
- if (intersectsPointGeometry(globalPos, scv.geometry()))
+ if (intersectsPointGeometry(globalPos, fvGeometry.geometry(scv)))
scvIndices.push_back(scv.indexInElement());
// for all scvs that tested positive add the point sources
diff --git a/dumux/discretization/box/fvelementgeometry.hh b/dumux/discretization/box/fvelementgeometry.hh
index 4424e4c64127ffc9c74c60327cd572a8ff50bacc..aeaf3c984e01a2db731d3d4370334e4b355ee3a9 100644
--- a/dumux/discretization/box/fvelementgeometry.hh
+++ b/dumux/discretization/box/fvelementgeometry.hh
@@ -28,6 +28,10 @@
#include <optional>
#include <utility>
+#include <unordered_map>
+#include <array>
+#include <vector>
+
#include <dune/geometry/type.hh>
#include <dune/localfunctions/lagrange/pqkfactory.hh>
@@ -59,6 +63,10 @@ class BoxFVElementGeometry<GG, true>
using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
using CoordScalar = typename GridView::ctype;
using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType;
+
+ using GeometryHelper = BoxGeometryHelper<GridView, dim,
+ typename GG::SubControlVolume,
+ typename GG::SubControlVolumeFace>;
public:
//! export the element type
using Element = typename GridView::template Codim<0>::Entity;
@@ -185,6 +193,20 @@ public:
bool hasBoundaryScvf() const
{ return gridGeometry().hasBoundaryScvf(eIdx_); }
+ //! Geometry of a sub control volume
+ typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const
+ {
+ assert(isBound());
+ return gridGeometryPtr_->geometry(element(), scv);
+ }
+
+ //! Geometry of a sub control volume face
+ typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const
+ {
+ assert(isBound());
+ return gridGeometryPtr_->geometry(element(), scvf);
+ }
+
private:
const GridGeometry* gridGeometryPtr_;
GridIndexType eIdx_;
@@ -331,7 +353,32 @@ public:
bool hasBoundaryScvf() const
{ return hasBoundaryScvf_; }
+ //! Geometry of a sub control volume
+ typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const
+ {
+ assert(isBound());
+ const auto geo = element().geometry();
+ return { Dune::GeometryTypes::cube(dim), 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())
+ {
+ const auto localBoundaryIndex = scvf.index() - numInnerScvf_(element());
+ const auto& key = scvfBoundaryGeometryKeys_[localBoundaryIndex];
+ return { Dune::GeometryTypes::cube(dim-1), GeometryHelper(geo).getBoundaryScvfCorners(key[0], key[1]) };
+ }
+ else
+ return { Dune::GeometryTypes::cube(dim-1), GeometryHelper(geo).getScvfCorners(scvf.index()) };
+ }
+
private:
+ unsigned int numInnerScvf_(const Element& element) const
+ { return (dim==1) ? 1 : element.subEntities(dim-1); }
void makeElementGeometries_()
{
@@ -349,7 +396,7 @@ private:
scvs_.resize(elementGeometry.corners());
for (LocalIndexType scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx)
{
- // get asssociated dof index
+ // get associated dof index
const auto dofIdxGlobal = gridGeometry().vertexMapper().subIndex(element, scvLocalIdx, dim);
// add scv to the local container
@@ -360,8 +407,9 @@ private:
}
// construct the sub control volume faces
- const auto numInnerScvf = (dim==1) ? 1 : element.subEntities(dim-1);
+ const auto numInnerScvf = numInnerScvf_(element);
scvfs_.resize(numInnerScvf);
+ scvfBoundaryGeometryKeys_.clear();
LocalIndexType scvfLocalIdx = 0;
for (; scvfLocalIdx < numInnerScvf; ++scvfLocalIdx)
@@ -400,6 +448,11 @@ private:
std::move(localScvIndices),
true);
+ scvfBoundaryGeometryKeys_.emplace_back(std::array<LocalIndexType, 2>{{
+ static_cast<LocalIndexType>(intersection.indexInInside()),
+ static_cast<LocalIndexType>(isScvfLocalIdx)
+ }});
+
// increment local counter
scvfLocalIdx++;
}
@@ -417,6 +470,7 @@ private:
//! vectors to store the geometries locally after binding an element
std::vector<SubControlVolume> scvs_;
std::vector<SubControlVolumeFace> scvfs_;
+ std::vector<std::array<LocalIndexType, 2>> scvfBoundaryGeometryKeys_;
bool hasBoundaryScvf_ = false;
};
diff --git a/dumux/discretization/box/fvgridgeometry.hh b/dumux/discretization/box/fvgridgeometry.hh
index b26b54d3c28b628e491361870960b93b90776968..3e74e7bf25084382c1fac22c92866584cdd8cf26 100644
--- a/dumux/discretization/box/fvgridgeometry.hh
+++ b/dumux/discretization/box/fvgridgeometry.hh
@@ -28,6 +28,8 @@
#include <utility>
#include <unordered_map>
+#include <array>
+#include <vector>
#include <dune/localfunctions/lagrange/lagrangelfecache.hh>
@@ -199,11 +201,36 @@ public:
bool hasBoundaryScvf(GridIndexType eIdx) const
{ return hasBoundaryScvf_[eIdx]; }
+ //! Geometry of a sub control volume
+ typename SubControlVolume::Traits::Geometry geometry(const Element& element, const SubControlVolume& scv) const
+ {
+ assert(isBound());
+ const auto geo = element.geometry();
+ return { Dune::GeometryTypes::cube(dim), GeometryHelper(geo).getScvCorners(scv.indexInElement()) };
+ }
+
+ //! Geometry of a sub control volume face
+ typename SubControlVolumeFace::Traits::Geometry geometry(const Element& element, const SubControlVolumeFace& scvf) const
+ {
+ assert(isBound());
+ const auto eIdx = this->elementMapper().index(element);
+ const auto geo = element.geometry();
+ if (scvf.boundary())
+ {
+ const auto localBoundaryIndex = scvf.index() - numInnerScvf_(element);
+ const auto& key = scvfBoundaryGeometryKeys_.at(eIdx)[localBoundaryIndex];
+ return { Dune::GeometryTypes::cube(dim-1), GeometryHelper(geo).getBoundaryScvfCorners(key[0], key[1]) };
+ }
+ else
+ return { Dune::GeometryTypes::cube(dim-1), GeometryHelper(geo).getScvfCorners(scvf.index()) };
+ }
+
private:
void update_()
{
scvs_.clear();
scvfs_.clear();
+ scvfBoundaryGeometryKeys_.clear();
auto numElements = this->gridView().size(0);
scvs_.resize(numElements);
@@ -287,6 +314,11 @@ private:
std::move(localScvIndices),
true);
+ scvfBoundaryGeometryKeys_[eIdx].emplace_back(std::array<LocalIndexType, 2>{{
+ static_cast<LocalIndexType>(intersection.indexInInside()),
+ static_cast<LocalIndexType>(isScvfLocalIdx)
+ }});
+
// increment local counter
scvfLocalIdx++;
}
@@ -347,10 +379,15 @@ private:
DUNE_THROW(Dune::NotImplemented, "Periodic boundaries for box method for parallel simulations!");
}
+ unsigned int numInnerScvf_(const Element& element) const
+ { return (dim==1) ? 1 : element.subEntities(dim-1); }
+
const FeCache feCache_;
std::vector<std::vector<SubControlVolume>> scvs_;
std::vector<std::vector<SubControlVolumeFace>> scvfs_;
+ std::unordered_map<GridIndexType, std::vector<std::array<LocalIndexType, 2>>> scvfBoundaryGeometryKeys_;
+
// TODO do we need those?
std::size_t numScv_;
std::size_t numScvf_;
diff --git a/test/discretization/box/test_boxfvgeometry.cc b/test/discretization/box/test_boxfvgeometry.cc
index 98ae68f315d698a7415b334424143bbe38a233fa..7550d624eda1eb8de9837e3776c915b3f43d47d6 100644
--- a/test/discretization/box/test_boxfvgeometry.cc
+++ b/test/discretization/box/test_boxfvgeometry.cc
@@ -100,6 +100,20 @@ 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 scvGeo = fvGeometry.geometry(scv);
+
+ if ((scv.center() - scvGeo.center()).two_norm2() > 1e-15)
+ DUNE_THROW(Dune::Exception,
+ "scv.center() and fvGeometry.geometry(scv).center() do not match! "
+ << scv.center() << " :: " << scvGeo.center()
+ );
+
+ if (std::abs(scv.volume() - scvGeo.volume()) > 1e-7)
+ DUNE_THROW(Dune::Exception,
+ "scv.volume() and fvGeometry.geometry(scv).volume() do not match! "
+ << scv.volume() << " :: " << scvGeo.volume()
+ );
}
auto range2 = scvfs(fvGeometry);
@@ -111,6 +125,21 @@ 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 scvfGeo = fvGeometry.geometry(scvf);
+
+ if ((scvf.ipGlobal() - scvfGeo.center()).two_norm2() > 1e-15)
+ DUNE_THROW(Dune::Exception,
+ "scvf.ipGlobal() and fvGeometry.geometry(scvf).center() do not match! "
+ << scvf.ipGlobal() << " :: " << scvfGeo.center()
+ );
+
+ if (std::abs(scvf.area() - scvfGeo.volume()) > 1e-7)
+ DUNE_THROW(Dune::Exception,
+ "scvf.area() and fvGeometry.geometry(scvf).volume() do not match! "
+ << scvf.area() << " :: " << scvfGeo.volume()
+ );
+
if (scvf.boundary())
{
++boundaryCount;