From ba299b504f64b638da8416f3a590686cb42e326a Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Sat, 12 Nov 2016 12:59:16 +0100
Subject: [PATCH] [staggeredGrid] Improve helper class
*Now finds adjacent normal faces
*Put some output in test
---
.../staggered/globalfvgeometry.hh | 2 +-
.../staggered/staggeredgeometryhelper.hh | 184 ++++++++++++++++--
.../staggered/subcontrolvolumeface.hh | 16 ++
.../staggered/test_staggered.cc | 4 +-
4 files changed, 189 insertions(+), 17 deletions(-)
diff --git a/dumux/discretization/staggered/globalfvgeometry.hh b/dumux/discretization/staggered/globalfvgeometry.hh
index 333265e638..d7e16024e5 100644
--- a/dumux/discretization/staggered/globalfvgeometry.hh
+++ b/dumux/discretization/staggered/globalfvgeometry.hh
@@ -64,7 +64,7 @@ class StaggeredGlobalFVGeometry<TypeTag, true>
dimWorld = GridView::dimensionworld
};
- using GeometryHelper = StaggeredGeometryHelper<GridView, dim>;
+ using GeometryHelper = StaggeredGeometryHelper<GridView>;
public:
//! Constructor
diff --git a/dumux/discretization/staggered/staggeredgeometryhelper.hh b/dumux/discretization/staggered/staggeredgeometryhelper.hh
index 5a81f1368f..45799401c9 100644
--- a/dumux/discretization/staggered/staggeredgeometryhelper.hh
+++ b/dumux/discretization/staggered/staggeredgeometryhelper.hh
@@ -28,23 +28,33 @@
#include <dune/geometry/referenceelements.hh>
#include <dumux/common/math.hh>
+#include <type_traits>
namespace Dumux
{
-//! Create sub control volumes and sub control volume face geometries
-template<class GridView, int dim>
-class StaggeredGeometryHelper;
+template<class Scalar>
+struct PairData
+{
+ std::pair<int,int> parallelPair;
+ std::pair<int,int> normalPair;
+ int localCommonEntIdx;
+ Scalar parallelDistance;
+ Scalar normalDistance;
+};
//! A class to create sub control volume and sub control volume face geometries per element
template <class GridView>
-class StaggeredGeometryHelper<GridView, 2>
+// class StaggeredGeometryHelper<GridView, dim>
+class StaggeredGeometryHelper
{
using Scalar = typename GridView::ctype;
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
+ static constexpr int numPairs = (dimWorld == 2) ? 2 : 4;
+
using ScvGeometry = Dune::CachedMultiLinearGeometry<Scalar, dim, dimWorld>;
using ScvfGeometry = Dune::CachedMultiLinearGeometry<Scalar, dim-1, dimWorld>;
@@ -56,40 +66,184 @@ class StaggeredGeometryHelper<GridView, 2>
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;
+
+ //TODO include assert that checks for quad geometry
+
public:
+
StaggeredGeometryHelper(const Intersection& intersection, const GridView& gridView)
- : intersection_(intersection), elementGeometry_(intersection.inside().geometry()), gridView_(gridView)//, corners_(geometry.corners())
+ : intersection_(intersection), elementGeometry_(intersection.inside().geometry()), gridView_(gridView), offset_(gridView.size(0))
{
+ fillPairData();
}
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;
+ return gridView_.indexSet().subIndex(intersection_.inside(), inIdx, dim-1) + offset_;
}
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;
+ return gridView_.indexSet().subIndex(intersection_.inside(), localOppositeIdx_(inIdx), dim-1) + offset_;
+ }
+
+
+ auto pairData() const
+ {
+ return pairData_;
+ }
+
+
+ void fillPairData()
+ {
+ const auto element = intersection_.inside();
+ const auto& referenceElement = ReferenceElements::general(element.geometry().type());
+
+ const int indexInInside = intersection_.indexInInside();
+
+
+ // store the element's own face info first
+ pairData_[0].parallelPair.first = dofIdxSelf();
+ pairData_[1].parallelPair.first = dofIdxSelf();
+
+ setInnerNormalPairs_(indexInInside);
+
+
+
+ // go into the direct neigbor element
+ if(intersection_.neighbor())
+ {
+ const int neighborIsIdx = intersection_.indexInOutside();
+
+ const auto& directNeighbor = intersection_.outside();
+
+
+ for(const auto& is : intersections(gridView_, directNeighbor))
+ {
+ constexpr int codimCommonEntity = 2;// TODO: 3d?
+
+ auto isNormalSide = [&](const int idx)
+ {
+ return !(idx == neighborIsIdx || localOppositeIdx_(idx) == neighborIsIdx);
+ };
+
+ auto localToGlobalEntityIdx = [&](const int localIdx)
+ {
+ return gridView_.indexSet().subIndex(directNeighbor, localIdx, codimCommonEntity);
+ };
+ // skip the directly neighboring face itself
+ if(isNormalSide(is.indexInInside()))
+ {
+ // get facet
+// const auto& facet = directNeighbor.template subEntity < 1 > (is.indexInInside());
+// std::cout <<
+
+ // iterate over facets sub-entities // TODO: get number correctly
+ for(int i = 0; i < 2; ++i)
+ {
+ int localCommonEntIdx = referenceElement.subEntity(is.indexInInside(), 1, i, dim);
+ const auto& commonEnt = directNeighbor.template subEntity < codimCommonEntity > (localCommonEntIdx);
+
+// int globalCommonEntIdx = gridView_.indexSet().subIndex(directNeighbor, localCommonEntIdx, codimCommonEntity);
+ int globalCommonEntIdx = localToGlobalEntityIdx(localCommonEntIdx);
+
+// int dofIdx = gridView_.indexSet().index(commonEnt);
+
+// if(globalCommonEntIdx != dofIdx)
+// std::cout << "error!";
+
+ std::cout << "pos: " << commonEnt.geometry().center() << std::endl;
+
+ if(globalCommonEntIdx == localToGlobalEntityIdx(pairData_[0].localCommonEntIdx))
+ {
+ pairData_[0].normalPair.second = gridView_.indexSet().subIndex(directNeighbor, is.indexInInside(), dim-1) + offset_;
+ }
+ if(globalCommonEntIdx == localToGlobalEntityIdx(pairData_[1].localCommonEntIdx))
+ {
+ pairData_[1].normalPair.second = gridView_.indexSet().subIndex(directNeighbor, is.indexInInside(), dim-1) + offset_;
+ }
+ }
+
+ std::cout << "inters: " << intersection_.geometry().center() <<
+ " || side : " << is.geometry().center() << std::endl;
+
+
+ }
+
+ }
+ }
+
+//TODO: fill parallel faces
}
private:
+
+ int localOppositeIdx_(const int idx) const
+ {
+ return (idx % 2) ? (idx - 1) : (idx + 1);
+ }
+
+
+ template<class G = GridView>
+ typename std::enable_if<G::dimension == 2, void>::type
+ setInnerNormalPairs_(const int isIdx)
+ {
+ auto& data1 = pairData_[0];
+ auto& data2 = pairData_[1];
+
+ switch(isIdx)
+ {
+ case 0:
+ data1.normalPair.first = gridView_.indexSet().subIndex(intersection_.inside(), 3, dim-1) + offset_;
+ data2.normalPair.first = gridView_.indexSet().subIndex(intersection_.inside(), 2, dim-1) + offset_;
+ data1.localCommonEntIdx = 2;
+ data2.localCommonEntIdx = 0;
+ break;
+ case 1:
+ data1.normalPair.first = gridView_.indexSet().subIndex(intersection_.inside(), 2, dim-1) + offset_;
+ data2.normalPair.first = gridView_.indexSet().subIndex(intersection_.inside(), 3, dim-1) + offset_;
+ data1.localCommonEntIdx = 1;
+ data2.localCommonEntIdx = 3;
+ break;
+ case 2:
+ data1.normalPair.first = gridView_.indexSet().subIndex(intersection_.inside(), 0, dim-1) + offset_;
+ data2.normalPair.first = gridView_.indexSet().subIndex(intersection_.inside(), 1, dim-1) + offset_;
+ data1.localCommonEntIdx = 0;
+ data2.localCommonEntIdx = 1;
+ break;
+ case 3:
+ data1.normalPair.first = gridView_.indexSet().subIndex(intersection_.inside(), 1, dim-1) + offset_;
+ data2.normalPair.first = gridView_.indexSet().subIndex(intersection_.inside(), 0, dim-1) + offset_;
+ data1.localCommonEntIdx = 3;
+ data2.localCommonEntIdx = 2;
+ break;
+ default:
+ DUNE_THROW(Dune::InvalidStateException, "foo");
+ }
+ }
+
+ template<class G = GridView>
+ typename std::enable_if<G::dimension == 3, void>::type
+ setInnerNormalPairs_(const int isIdx)
+ {
+ // TODO: 3D
+ DUNE_THROW(Dune::NotImplemented, "3d helper not ready yet");
+ }
+
+
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
+ const int offset_;
+
+ std::array<PairData<Scalar>, numPairs> pairData_;
+
};
diff --git a/dumux/discretization/staggered/subcontrolvolumeface.hh b/dumux/discretization/staggered/subcontrolvolumeface.hh
index b00a83e396..2a65374624 100644
--- a/dumux/discretization/staggered/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/subcontrolvolumeface.hh
@@ -26,8 +26,11 @@
#include <utility>
#include <dune/common/fvector.hh>
#include <dumux/discretization/subcontrolvolumefacebase.hh>
+#include <dumux/discretization/staggered/staggeredgeometryhelper.hh>
#include <dumux/common/optional.hh>
+#include <typeinfo>
+
namespace Dumux
{
@@ -67,6 +70,9 @@ class StaggeredSubControlVolumeFace : public SubControlVolumeFaceBase<StaggeredS
using StaggeredSubFace = Dumux::StaggeredSubFace<G,I>;
+ static constexpr int numPairs = (dimworld == 2) ? 2 : 4;
+
+
public:
// the default constructor
StaggeredSubControlVolumeFace() = default;
@@ -94,6 +100,8 @@ public:
selfIdx_ = geometryHelper.dofIdxSelf();
oppositeIdx_ = geometryHelper.dofIdxOpposite();
+
+ pairData_ = geometryHelper.pairData();
}
/*//! The copy constrcutor
@@ -193,6 +201,12 @@ public:
return oppositeIdx_;
}
+
+ auto pairData(const int idx) const
+ {
+ return pairData_[idx];
+ }
+
private:
Dune::GeometryType geomType_;
std::vector<GlobalPosition> corners_;
@@ -206,6 +220,8 @@ private:
int selfIdx_;
int oppositeIdx_;
std::vector<StaggeredSubFace> subfaces_;
+ std::array<PairData<Scalar>, numPairs> pairData_;
+
};
diff --git a/test/discretization/staggered/test_staggered.cc b/test/discretization/staggered/test_staggered.cc
index 2a798b6412..9998cc4697 100644
--- a/test/discretization/staggered/test_staggered.cc
+++ b/test/discretization/staggered/test_staggered.cc
@@ -131,9 +131,11 @@ int main (int argc, char *argv[]) try
if(0 != testForwardIterator(range2.begin(), range2.end(), op2))
DUNE_THROW(Dune::Exception, "Iterator does not fulfill the forward iterator concept");
+//TODO: beautify output
for (auto&& scvf : scvfs(fvGeometry))
{
- std::cout << "-- scvf " << scvf.index() << " ip at: " << scvf.ipGlobal() << ", doIdx (self): " << scvf.dofIndexSelf() << ", doIdx (opposite): " << scvf.dofIndexOpposite();
+ std::cout << "-- scvf " << scvf.index() << " ip at: " << scvf.ipGlobal() << ", self: " << scvf.dofIndexSelf() << ", opposite: " << scvf.dofIndexOpposite() << ", side1 in: " << scvf.pairData(0).normalPair.first << ", side2 in: " << scvf.pairData(1).normalPair.first << " , side1 out: " << scvf.pairData(0).normalPair.second
+ << ", side2out: " << scvf.pairData(1).normalPair.second;
if (scvf.boundary()) std::cout << " (on boundary).";
std::cout << std::endl;
}
--
GitLab