Skip to content
Snippets Groups Projects
Commit 671c98e7 authored by Kilian Weishaupt's avatar Kilian Weishaupt Committed by Timo Koch
Browse files

[staggeredGrid] Clean up geometryHelper 

*Use separate 2d/3d classes that inherit from a base class
*Add some docu
parent 89a7347b
No related branches found
No related tags found
2 merge requests!617[WIP] Next,!370Feature/staggered grid
Hide whitespace changes
Inline Side-by-side
dumux/discretization/staggered/staggeredgeometryhelper.hh
+ 137
− 69
View file @ 671c98e7
...@@ -43,11 +43,13 @@ struct PairData ...@@ -43,11 +43,13 @@ struct PairData
Scalar normalDistance; Scalar normalDistance;
}; };
//! A class to create sub control volume and sub control volume face geometries per element //! A class to create sub control volume and sub control volume face geometries per element
template <class GridView> template <class GridView, int dim = GridView::dimension >
// class StaggeredGeometryHelper<GridView, dim>
class StaggeredGeometryHelper class StaggeredGeometryHelper
{};
template<class GridView>
class BaseStaggeredGeometryHelper
{ {
using Scalar = typename GridView::ctype; using Scalar = typename GridView::ctype;
static constexpr int dim = GridView::dimension; static constexpr int dim = GridView::dimension;
...@@ -66,20 +68,22 @@ class StaggeredGeometryHelper ...@@ -66,20 +68,22 @@ class StaggeredGeometryHelper
using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>; using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
//TODO include assert that checks for quad geometry //TODO include assert that checks for quad geometry
static constexpr int codimCommonEntity = 2; //TODO: 3d? static constexpr int codimCommonEntity = 2; //TODO: 3d?
static constexpr int numFacetSubEntities = 2; // TODO: 3d?
using Implementation = typename Dumux::StaggeredGeometryHelper<GridView, dim>;
public: public:
BaseStaggeredGeometryHelper(const Intersection& intersection, const GridView& gridView)
StaggeredGeometryHelper(const Intersection& intersection, const GridView& gridView)
: intersection_(intersection), element_(intersection.inside()), elementGeometry_(element_.geometry()), gridView_(gridView), offset_(gridView.size(0)) : intersection_(intersection), element_(intersection.inside()), elementGeometry_(element_.geometry()), gridView_(gridView), offset_(gridView.size(0))
{ {
fillPairData(); fillPairData();
} }
/*!
* \brief Returns the global dofIdx of the intersection itself
*/
int dofIdxSelf() const int dofIdxSelf() const
{ {
//TODO: use proper intersection mapper! //TODO: use proper intersection mapper!
...@@ -87,20 +91,27 @@ public: ...@@ -87,20 +91,27 @@ public:
return gridView_.indexSet().subIndex(intersection_.inside(), inIdx, dim-1) + offset_; return gridView_.indexSet().subIndex(intersection_.inside(), inIdx, dim-1) + offset_;
} }
/*!
* \brief Returns the global dofIdx of the opposing intersection
*/
int dofIdxOpposite() const int dofIdxOpposite() const
{ {
//TODO: use proper intersection mapper! //TODO: use proper intersection mapper!
const auto inIdx = intersection_.indexInInside(); const auto inIdx = intersection_.indexInInside();
return gridView_.indexSet().subIndex(intersection_.inside(), localOppositeIdx_(inIdx), dim-1) + offset_; return gridView_.indexSet().subIndex(this->intersection_.inside(), localOppositeIdx_(inIdx), dim-1) + this->offset_;
} }
/*!
* \brief Returns a copy of the pair data
*/
auto pairData() const auto pairData() const
{ {
return pairData_; return pairData_;
} }
/*!
* \brief Fills all entries of the pair data
*/
void fillPairData() void fillPairData()
{ {
const auto& referenceElement = ReferenceElements::general(element_.geometry().type()); const auto& referenceElement = ReferenceElements::general(element_.geometry().type());
...@@ -114,10 +125,9 @@ public: ...@@ -114,10 +125,9 @@ public:
data.parallelDistance = 0; data.parallelDistance = 0;
} }
// set the inner parts of the normal pairs // set the inner parts of the normal pairs
const auto localInnerNormalDofIndices = getLocalInnerNormalDofIndices_(indexInInside); const auto localInnerNormalDofIndices = asImp_().getLocalInnerNormalDofIndices_(indexInInside);
setInnerNormalPairs_(localInnerNormalDofIndices); asImp_().setInnerNormalPairs_(localInnerNormalDofIndices);
// get the positions of the faces normal to the intersection within the element // get the positions of the faces normal to the intersection within the element
std::vector<GlobalPosition> innerNormalFacePos; std::vector<GlobalPosition> innerNormalFacePos;
...@@ -134,7 +144,6 @@ public: ...@@ -134,7 +144,6 @@ public:
DUNE_THROW(Dune::NotImplemented, "3d not ready yet"); DUNE_THROW(Dune::NotImplemented, "3d not ready yet");
} }
// go into the direct neighbor element // go into the direct neighbor element
if(intersection_.neighbor()) if(intersection_.neighbor())
{ {
...@@ -147,8 +156,8 @@ public: ...@@ -147,8 +156,8 @@ public:
// skip the directly neighboring face itself and its opposing one // skip the directly neighboring face itself and its opposing one
if(neighborIntersectionNormalSide_(neighborIsIdx, intersection_.indexInOutside())) if(neighborIntersectionNormalSide_(neighborIsIdx, intersection_.indexInOutside()))
{ {
// iterate over facets sub-entities // TODO: get number correctly // iterate over facets sub-entities
for(int i = 0; i < 2; ++i) for(int i = 0; i < numFacetSubEntities; ++i)
{ {
int localCommonEntIdx = referenceElement.subEntity(neighborIsIdx, 1, i, dim); int localCommonEntIdx = referenceElement.subEntity(neighborIsIdx, 1, i, dim);
int globalCommonEntIdx = localToGlobalEntityIdx_(localCommonEntIdx, directNeighbor); int globalCommonEntIdx = localToGlobalEntityIdx_(localCommonEntIdx, directNeighbor);
...@@ -166,8 +175,6 @@ public: ...@@ -166,8 +175,6 @@ public:
} }
} }
// go into the adjacent neighbor element // go into the adjacent neighbor element
if(neighborIntersection.neighbor()) if(neighborIntersection.neighbor())
{ {
...@@ -176,7 +183,7 @@ public: ...@@ -176,7 +183,7 @@ public:
{ {
if(neighborIntersectionNormalSide_(dIs.indexInInside(), neighborIntersection.indexInOutside())) if(neighborIntersectionNormalSide_(dIs.indexInInside(), neighborIntersection.indexInOutside()))
{ {
for(int i = 0; i < 2; ++i) for(int i = 0; i < numFacetSubEntities; ++i)
{ {
int localCommonEntIdx = referenceElement.subEntity(dIs.indexInInside(), 1, i, dim); int localCommonEntIdx = referenceElement.subEntity(dIs.indexInInside(), 1, i, dim);
int globalCommonEntIdx = localToGlobalEntityIdx_(localCommonEntIdx, diagonalNeighbor); int globalCommonEntIdx = localToGlobalEntityIdx_(localCommonEntIdx, diagonalNeighbor);
...@@ -202,28 +209,89 @@ public: ...@@ -202,28 +209,89 @@ public:
} }
private: private:
/*!
* \brief Returns the local opposing intersection index
*
* \param idx The local index of the intersection itself
*/
int localOppositeIdx_(const int idx) const int localOppositeIdx_(const int idx) const
{ {
return (idx % 2) ? (idx - 1) : (idx + 1); return (idx % 2) ? (idx - 1) : (idx + 1);
} }
bool neighborIntersectionNormalSide_(const int isIdx, const int neighborIsIdx) const /*!
* \brief Returns true if the intersection lies normal to another given intersection
*
* \param selfIdx The local index of the intersection itself
* \param otherIdx The local index of the other intersection
*/
bool neighborIntersectionNormalSide_(const int selfIdx, const int otherIdx) const
{ {
return !(isIdx == neighborIsIdx || localOppositeIdx_(isIdx) == neighborIsIdx); return !(selfIdx == otherIdx || localOppositeIdx_(selfIdx) == otherIdx);
}; };
int localToGlobalEntityIdx_(const int localIdx, const Element& element) /*!
* \brief Returns the global index of the common entity
*
* \param localIdx The local index of the common entity
* \param element The element
*/
int localToGlobalEntityIdx_(const int localIdx, const Element& element) const
{ {
return gridView_.indexSet().subIndex(element, localIdx, codimCommonEntity); return this->gridView_.indexSet().subIndex(element, localIdx, codimCommonEntity);
}; };
protected:
const Intersection& intersection_; //! The intersection of interest
const Element& element_; //! The respective element
const typename Element::Geometry& elementGeometry_; //! Reference to the element geometry
const GridView gridView_;
const int offset_; //! Offset for intersection dof indexing
std::array<PairData<Scalar>, numPairs> pairData_; //! collection of pair information
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation &asImp_()
{ return *static_cast<Implementation *>(this); }
//! \copydoc asImp_()
const Implementation &asImp_() const
{ return *static_cast<const Implementation *>(this); }
};
// specializations for 2D *************************************************************************************************** template<class GridView>
template<class G = GridView, typename std::enable_if<G::dimension == 2, int>::type = 0> class StaggeredGeometryHelper<GridView, 2> : public BaseStaggeredGeometryHelper<GridView>
auto getLocalInnerNormalDofIndices_(const int directNeighborIsIdx) {
friend class BaseStaggeredGeometryHelper<GridView>;
using Scalar = typename GridView::ctype;
static constexpr int dim = GridView::dimension;
static constexpr 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>;
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>;
//TODO include assert that checks for quad geometry
static constexpr int codimCommonEntity = 2; //TODO: 3d?
using ParentType = BaseStaggeredGeometryHelper<GridView>;
public:
StaggeredGeometryHelper(const Intersection& intersection, const GridView& gridView)
: ParentType(intersection, gridView)
{}
private:
static auto getLocalInnerNormalDofIndices_(const int directNeighborIsIdx)
{ {
struct Indices struct Indices
{ {
...@@ -267,60 +335,60 @@ private: ...@@ -267,60 +335,60 @@ private:
return indices; return indices;
} }
template<class Indices, class G = GridView> template<class Indices>
typename std::enable_if<G::dimension == 2, void>::type void setInnerNormalPairs_(const Indices& indices)
setInnerNormalPairs_(const Indices& indices)
{ {
pairData_[0].normalPair.first = gridView_.indexSet().subIndex(intersection_.inside(), indices.normalLocalDofIdx1, dim-1) + offset_; this->pairData_[0].normalPair.first = this->gridView_.indexSet().subIndex(this->intersection_.inside(), indices.normalLocalDofIdx1, dim-1) + this->offset_;
pairData_[1].normalPair.first = gridView_.indexSet().subIndex(intersection_.inside(), indices.normalLocalDofIdx2, dim-1) + offset_; this->pairData_[1].normalPair.first = this->gridView_.indexSet().subIndex(this->intersection_.inside(), indices.normalLocalDofIdx2, dim-1) + this->offset_;
pairData_[0].globalCommonEntIdx = gridView_.indexSet().subIndex(intersection_.inside(), indices.localCommonEntIdx1, codimCommonEntity); this->pairData_[0].globalCommonEntIdx = this->gridView_.indexSet().subIndex(this->intersection_.inside(), indices.localCommonEntIdx1, codimCommonEntity);
pairData_[1].globalCommonEntIdx = gridView_.indexSet().subIndex(intersection_.inside(), indices.localCommonEntIdx2, codimCommonEntity); this->pairData_[1].globalCommonEntIdx = this->gridView_.indexSet().subIndex(this->intersection_.inside(), indices.localCommonEntIdx2, codimCommonEntity);
} }
};
template<class G = GridView, typename std::enable_if<G::dimension == 3, int>::type = 0> template<class GridView>
auto getLocalInnerNormalDofIndices_(const int directNeighborIsIdx) class StaggeredGeometryHelper<GridView, 3> : public BaseStaggeredGeometryHelper<GridView>
{ {
struct Indices friend class BaseStaggeredGeometryHelper<GridView>;
{ using Scalar = typename GridView::ctype;
int normalLocalDofIdx1; static constexpr int dim = GridView::dimension;
int normalLocalDofIdx2; static constexpr int dimWorld = GridView::dimensionworld;
int normalLocalDofIdx3;
int normalLocalDofIdx4;
int localCommonEntIdx1;
int localCommonEntIdx2;
int localCommonEntIdx3;
int localCommonEntIdx4;
};
Indices indices; using ScvGeometry = Dune::CachedMultiLinearGeometry<Scalar, dim, dimWorld>;
using ScvfGeometry = Dune::CachedMultiLinearGeometry<Scalar, dim-1, dimWorld>;
switch(directNeighborIsIdx) using GlobalPosition = typename ScvGeometry::GlobalCoordinate;
{ using PointVector = std::vector<GlobalPosition>;
default:
DUNE_THROW(Dune::NotImplemented, "3d helper not ready yet");
}
return indices;
}
template<class Indices, class G = GridView> using Element = typename GridView::template Codim<0>::Entity;
typename std::enable_if<G::dimension == 3, void>::type using Intersection = typename GridView::Intersection;
setInnerNormalPairs_(const Indices& indices)
{
// TODO: 3D
DUNE_THROW(Dune::NotImplemented, "3d helper not ready yet");
}
using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
const Intersection& intersection_;
const Element& element_;
const typename Element::Geometry& elementGeometry_; //! Reference to the element geometry
const GridView gridView_;
const int offset_;
std::array<PairData<Scalar>, numPairs> pairData_; //TODO include assert that checks for quad geometry
static constexpr int codimCommonEntity = 2; //TODO: 3d?
using ParentType = BaseStaggeredGeometryHelper<GridView>;
public:
StaggeredGeometryHelper(const Intersection& intersection, const GridView& gridView)
: ParentType(intersection, gridView)
{}
private:
auto getLocalInnerNormalDofIndices_(const int directNeighborIsIdx)
{
// TODO: 3D
DUNE_THROW(Dune::NotImplemented, "3d helper not ready yet");
}
template<class Indices>
void setInnerNormalPairs_(const Indices& indices)
{
// TODO: 3D
DUNE_THROW(Dune::NotImplemented, "3d helper not ready yet");
}
}; };
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment