diff --git a/dumux/discretization/staggered/staggeredgeometryhelper.hh b/dumux/discretization/staggered/staggeredgeometryhelper.hh
index 614507c36c30fa0e13b104888af74a2842c0a6e4..a7fb1d7787a61f26bc0d2cf9432572c06b81d259 100644
--- a/dumux/discretization/staggered/staggeredgeometryhelper.hh
+++ b/dumux/discretization/staggered/staggeredgeometryhelper.hh
@@ -79,7 +79,7 @@ public:
BaseStaggeredGeometryHelper(const Intersection& intersection, const GridView& gridView)
: intersection_(intersection), element_(intersection.inside()), elementGeometry_(element_.geometry()), gridView_(gridView), offset_(gridView.size(0))
{
- fillPairData();
+ fillPairData_();
}
/*!
@@ -121,10 +121,11 @@ public:
return pairData_;
}
+private:
/*!
* \brief Fills all entries of the pair data
*/
- void fillPairData()
+ void fillPairData_()
{
const auto& referenceElement = ReferenceElements::general(element_.geometry().type());
const int indexInInside = intersection_.indexInInside();
@@ -132,10 +133,10 @@ public:
// initialize values that could remain unitialized if the intersection lies on a boundary
for(auto& data : pairData_)
{
- data.outerParallelFaceDofIdx = 0;
- data.outerParallelElementDofIdx = 0;
- data.normalDistance = 0;
- data.parallelDistance = 0;
+ data.outerParallelFaceDofIdx = -1;
+ data.outerParallelElementDofIdx = -1;
+ data.normalDistance = -1;
+ data.parallelDistance = -1;
}
// set the inner parts of the normal pairs
@@ -167,13 +168,13 @@ public:
{
const int neighborIsIdx = neighborIntersection.indexInInside();
// skip the directly neighboring face itself and its opposing one
- if(neighborIntersectionNormalSide_(neighborIsIdx, intersection_.indexInOutside()))
+ if(facetIsNormal_(neighborIsIdx, intersection_.indexInOutside()))
{
// iterate over facets sub-entities
for(int i = 0; i < numFacetSubEntities; ++i)
{
int localCommonEntIdx = referenceElement.subEntity(neighborIsIdx, 1, i, dim);
- int globalCommonEntIdx = localToGlobalEntityIdx_(localCommonEntIdx, directNeighbor);
+ int globalCommonEntIdx = localToGlobalCommonEntityIdx_(localCommonEntIdx, directNeighbor);
// fill the normal pair entries
for(int pairIdx = 0; pairIdx < numPairs; ++pairIdx)
@@ -194,12 +195,12 @@ public:
const auto& diagonalNeighbor = neighborIntersection.outside();
for(const auto& dIs : intersections(gridView_, diagonalNeighbor))
{
- if(neighborIntersectionNormalSide_(dIs.indexInInside(), neighborIntersection.indexInOutside()))
+ if(facetIsNormal_(dIs.indexInInside(), neighborIntersection.indexInOutside()))
{
for(int i = 0; i < numFacetSubEntities; ++i)
{
int localCommonEntIdx = referenceElement.subEntity(dIs.indexInInside(), 1, i, dim);
- int globalCommonEntIdx = localToGlobalEntityIdx_(localCommonEntIdx, diagonalNeighbor);
+ int globalCommonEntIdx = localToGlobalCommonEntityIdx_(localCommonEntIdx, diagonalNeighbor);
for(int pairIdx = 0; pairIdx < numPairs; ++pairIdx)
@@ -207,7 +208,7 @@ public:
if(globalCommonEntIdx == pairData_[pairIdx].globalCommonEntIdx)
{
pairData_[pairIdx].outerParallelFaceDofIdx = gridView_.indexSet().subIndex(diagonalNeighbor, dIs.indexInInside(), dim-1) + offset_;
- pairData_[pairIdx].outerParallelElementDofIdx = dIs.indexInOutside();
+ pairData_[pairIdx].outerParallelElementDofIdx = gridView_.indexSet().index(dIs.outside());
const auto& selfFacet = element_.template subEntity <1> (indexInInside);
const auto& parallelFacet = diagonalNeighbor.template subEntity <1> (dIs.indexInInside());
pairData_[pairIdx].parallelDistance = (selfFacet.geometry().center() - parallelFacet.geometry().center()).two_norm();
@@ -220,9 +221,42 @@ public:
}
}
}
+ else // intersection is on boundary
+ {
+ // find an intersection normal to the face
+ for(const auto& normalIntersection : intersections(gridView_, element_))
+ {
+ if(facetIsNormal_(normalIntersection.indexInInside(), intersection_.indexInInside()) && normalIntersection.neighbor())
+ {
+ const auto& neighbor = normalIntersection.outside();
+
+ for(const auto& neighborIs : intersections(gridView_, neighbor))
+ {
+ // iterate over facets sub-entities
+ if(neighborIs.indexInInside() != normalIntersection.indexInOutside())
+ {
+ for(int i = 0; i < numFacetSubEntities; ++i)
+ {
+ int localCommonEntIdx = referenceElement.subEntity(neighborIs.indexInInside(), 1, i, dim);
+ int globalCommonEntIdx = localToGlobalCommonEntityIdx_(localCommonEntIdx, neighbor);
+
+ // fill the parallel pair entries
+ for(int pairIdx = 0; pairIdx < numPairs; ++pairIdx)
+ {
+ if(globalCommonEntIdx == pairData_[pairIdx].globalCommonEntIdx)
+ {
+ pairData_[pairIdx].outerParallelFaceDofIdx = gridView_.indexSet().subIndex(neighbor, neighborIs.indexInInside(), dim-1) + offset_;
+ pairData_[pairIdx].outerParallelElementDofIdx = gridView_.indexSet().index(neighbor);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
}
-private:
/*!
* \brief Returns the local opposing intersection index
*
@@ -239,7 +273,7 @@ private:
* \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
+ bool facetIsNormal_(const int selfIdx, const int otherIdx) const
{
return !(selfIdx == otherIdx || localOppositeIdx_(selfIdx) == otherIdx);
};
@@ -250,7 +284,7 @@ private:
* \param localIdx The local index of the common entity
* \param element The element
*/
- int localToGlobalEntityIdx_(const int localIdx, const Element& element) const
+ int localToGlobalCommonEntityIdx_(const int localIdx, const Element& element) const
{
return this->gridView_.indexSet().subIndex(element, localIdx, codimCommonEntity);
};
diff --git a/dumux/freeflow/staggered/fluxvariables.hh b/dumux/freeflow/staggered/fluxvariables.hh
index 5107f731df43734153610e68cd89a3f7f6e59e42..f832c6c321fe2a5ec1042f25a739303f85cd57ac 100644
--- a/dumux/freeflow/staggered/fluxvariables.hh
+++ b/dumux/freeflow/staggered/fluxvariables.hh
@@ -112,18 +112,15 @@ public:
const SubControlVolumeFace& scvFace)
{
Stencil cellCenterStencil;
+
+ // the cell center dof indices
+ cellCenterStencil.push_back(scvFace.insideScvIdx());
+
+ // the face dof indices
+ cellCenterStencil.push_back(scvFace.dofIndexSelf());
if (!scvFace.boundary())
- {
- // the cell center dof indices
- cellCenterStencil.push_back(scvFace.insideScvIdx());
cellCenterStencil.push_back(scvFace.outsideScvIdx());
- // the face dof indices
- cellCenterStencil.push_back(scvFace.dofIndexSelf());
- }
- else
- cellCenterStencil.push_back(scvFace.insideScvIdx());
-
return cellCenterStencil;
}
@@ -131,27 +128,29 @@ public:
const SubControlVolumeFace& scvFace)
{
Stencil faceStencil;
+
+ // normal element dof indices
+ faceStencil.push_back(scvFace.insideScvIdx());
if (!scvFace.boundary())
- {
- // the cell center dof indices normal to the face
- faceStencil.push_back(scvFace.insideScvIdx());
faceStencil.push_back(scvFace.outsideScvIdx());
- // the cell center dof indices parallel to the face
- for(const auto& data : scvFace.pairData())
- {
- faceStencil.push_back(data.outerParallelElementDofIdx);
- faceStencil.push_back(data.outerParallelFaceDofIdx);
- faceStencil.push_back(data.normalPair.first);
- faceStencil.push_back(data.normalPair.second);
- }
+ // the normal face dof indices
+ faceStencil.push_back(scvFace.dofIndexSelf());
+ faceStencil.push_back(scvFace.dofIndexOpposite());
- // the face dof indices
- faceStencil.push_back(scvFace.dofIndexSelf());
- faceStencil.push_back(scvFace.dofIndexOpposite());
+ for(const auto& data : scvFace.pairData())
+ {
+ const auto& outerParallelElementDofIdx = data.outerParallelElementDofIdx;
+ const auto& outerParallelFaceDofIdx = data.outerParallelFaceDofIdx;
+ if(outerParallelElementDofIdx >= 0)
+ faceStencil.push_back(outerParallelElementDofIdx);
+ if(outerParallelFaceDofIdx >= 0)
+ faceStencil.push_back(outerParallelFaceDofIdx);
+
+ faceStencil.push_back(data.normalPair.first);
+ if(!scvFace.boundary())
+ faceStencil.push_back(data.normalPair.second);
}
- else
- faceStencil.push_back(scvFace.dofIndexSelf());
return faceStencil;
}