diff --git a/dumux/discretization/staggered/globalfvgeometry.hh b/dumux/discretization/staggered/globalfvgeometry.hh
index bcfb8d03c1629d4b6b543f8cc2702d60372d8e17..333265e638312f163b071905f21784e0d36b5daf 100644
--- a/dumux/discretization/staggered/globalfvgeometry.hh
+++ b/dumux/discretization/staggered/globalfvgeometry.hh
@@ -28,6 +28,7 @@
#include <dumux/common/elementmap.hh>
#include <dumux/implicit/staggered/properties.hh>
#include <dumux/discretization/staggered/fvelementgeometry.hh>
+#include <dumux/discretization/staggered/staggeredgeometryhelper.hh>
namespace Dumux
{
@@ -63,6 +64,8 @@ class StaggeredGlobalFVGeometry<TypeTag, true>
dimWorld = GridView::dimensionworld
};
+ using GeometryHelper = StaggeredGeometryHelper<GridView, dim>;
+
public:
//! Constructor
StaggeredGlobalFVGeometry(const GridView& gridView)
@@ -110,7 +113,6 @@ public:
elementMap_.clear();
// determine size of containers
- const int numElements = gridView_.size(0);
IndexType numScvs = gridView_.size(0);
IndexType numScvf = 0;
for (const auto& element : elements(gridView_))
@@ -128,6 +130,12 @@ public:
for (const auto& element : elements(gridView_))
{
auto eIdx = problem.elementMapper().index(element);
+
+ // get the element geometry
+// auto elementGeometry = element.geometry();
+// GeometryHelper geometryHelper(elementGeometry);
+
+
scvs_[eIdx] = SubControlVolume(element.geometry(), eIdx);
// fill the element map with seeds
@@ -138,22 +146,8 @@ public:
scvfsIndexSet.reserve(element.subEntities(1));
for (const auto& intersection : intersections(gridView_, element))
{
- //TODO: use proper intersection mapper!
- const auto inIdx = intersection.indexInInside();
- const auto globalIsIdx = gridView_.indexSet().subIndex(element, inIdx, dim-1) + numElements;
-
- int oppoSiteIdx;
-
- if(inIdx % 2) // face index is odd
- {
- oppoSiteIdx = inIdx -1;
- }
- else
- {
- oppoSiteIdx = inIdx +1;
- }
+ GeometryHelper geometryHelper(intersection, gridView_);
-// std::cout << "faceSelf: " << inIdx << ", oppo: " << oppoSiteIdx << std::endl;
// inner sub control volume faces
if (intersection.neighbor())
{
@@ -162,7 +156,7 @@ public:
intersection.geometry(),
scvfIdx,
std::vector<IndexType>({eIdx, nIdx}),
- globalIsIdx
+ geometryHelper
);
scvfsIndexSet.push_back(scvfIdx++);
}
@@ -173,7 +167,7 @@ public:
intersection.geometry(),
scvfIdx,
std::vector<IndexType>({eIdx, gridView_.size(0) + numBoundaryScvf_++}),
- globalIsIdx
+ geometryHelper
);
scvfsIndexSet.push_back(scvfIdx++);
}
diff --git a/dumux/discretization/staggered/staggeredgeometryhelper.hh b/dumux/discretization/staggered/staggeredgeometryhelper.hh
index 6d7339bbde20ce589f657a5e78c8abad3d27a054..5a81f1368f8bc4cfbaf7236c18e223a282ad35bf 100644
--- a/dumux/discretization/staggered/staggeredgeometryhelper.hh
+++ b/dumux/discretization/staggered/staggeredgeometryhelper.hh
@@ -32,8 +32,67 @@
namespace Dumux
{
+//! Create sub control volumes and sub control volume face geometries
+template<class GridView, int dim>
+class StaggeredGeometryHelper;
+//! A class to create sub control volume and sub control volume face geometries per element
+template <class GridView>
+class StaggeredGeometryHelper<GridView, 2>
+{
+ using Scalar = typename GridView::ctype;
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+
+ using ScvGeometry = Dune::CachedMultiLinearGeometry<Scalar, dim, dimWorld>;
+ using ScvfGeometry = Dune::CachedMultiLinearGeometry<Scalar, dim-1, dimWorld>;
+
+ using GlobalPosition = typename ScvGeometry::GlobalCoordinate;
+ using PointVector = std::vector<GlobalPosition>;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Intersection = typename GridView::Intersection;
+
+ 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
+ static constexpr int maxPoints = 9;
+
+public:
+
+ StaggeredGeometryHelper(const Intersection& intersection, const GridView& gridView)
+ : intersection_(intersection), elementGeometry_(intersection.inside().geometry()), gridView_(gridView)//, corners_(geometry.corners())
+ {
+ }
+
+ int dofIdxSelf() const
+ {
+ //TODO: use proper intersection mapper!
+ const auto inIdx = intersection_.indexInInside();
+ const int numElements = gridView_.size(0);
+ return gridView_.indexSet().subIndex(intersection_.inside(), inIdx, dim-1) + numElements;
+ }
+
+ int dofIdxOpposite() const
+ {
+ //TODO: use proper intersection mapper!
+ const auto inIdx = intersection_.indexInInside();
+ const int numElements = gridView_.size(0);
+ const int localOppositeIdx = (inIdx % 2) ? (inIdx - 1) : (inIdx + 1);
+ return gridView_.indexSet().subIndex(intersection_.inside(), localOppositeIdx, dim-1) + numElements;
+ }
+
+private:
+ const Intersection& intersection_;
+ const typename Element::Geometry& elementGeometry_; //! Reference to the element geometry
+ const GridView gridView_;
+// GlobalPosition p[maxPoints]; // the points needed for construction of the geometries
+// std::size_t corners_; // number of element corners
+
+};
+
} // end namespace Dumux
#endif
diff --git a/dumux/discretization/staggered/subcontrolvolumeface.hh b/dumux/discretization/staggered/subcontrolvolumeface.hh
index a137e17722c9ebb04dd5ce1bf593570ce23b0922..b00a83e396dca4d84a58990bba58db13e2d8ad74 100644
--- a/dumux/discretization/staggered/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/subcontrolvolumeface.hh
@@ -72,12 +72,12 @@ public:
StaggeredSubControlVolumeFace() = default;
//! Constructor with intersection
- template <class Intersection>
+ template <class Intersection, class GeometryHelper>
StaggeredSubControlVolumeFace(const Intersection& is,
const typename Intersection::Geometry& isGeometry,
IndexType scvfIndex,
const std::vector<IndexType>& scvIndices,
- const int selfIdx
+ const GeometryHelper& geometryHelper
)
: ParentType(),
geomType_(isGeometry.type()),
@@ -86,19 +86,14 @@ public:
unitOuterNormal_(is.centerUnitOuterNormal()),
scvfIndex_(scvfIndex),
scvIndices_(scvIndices),
- boundary_(is.boundary()),
- selfIdx_(selfIdx)
+ boundary_(is.boundary())
{
corners_.resize(isGeometry.corners());
for (int i = 0; i < isGeometry.corners(); ++i)
corners_[i] = isGeometry.corner(i);
-// subfaces_.resize(2);
-
-// const auto& element = is.inside();
-// const int inIdx = is.indexInInside();
-// gridView.indexSet().subIndex(element, inIdx, dimWorld-1);
-
+ selfIdx_ = geometryHelper.dofIdxSelf();
+ oppositeIdx_ = geometryHelper.dofIdxOpposite();
}
/*//! The copy constrcutor
diff --git a/test/discretization/staggered/test_staggered.cc b/test/discretization/staggered/test_staggered.cc
index 1c3a0f9217926162121325e84f2d01c54d5d7b51..2a798b6412ae151daaa1179fbc8942a555721a2c 100644
--- a/test/discretization/staggered/test_staggered.cc
+++ b/test/discretization/staggered/test_staggered.cc
@@ -133,7 +133,7 @@ int main (int argc, char *argv[]) try
for (auto&& scvf : scvfs(fvGeometry))
{
- std::cout << "-- scvf " << scvf.index() << " ip at: " << scvf.ipGlobal();
+ std::cout << "-- scvf " << scvf.index() << " ip at: " << scvf.ipGlobal() << ", doIdx (self): " << scvf.dofIndexSelf() << ", doIdx (opposite): " << scvf.dofIndexOpposite();
if (scvf.boundary()) std::cout << " (on boundary).";
std::cout << std::endl;
}