Commit 134d602c authored by Melanie Lipp's avatar Melanie Lipp Committed by Ned Coltman
Browse files

[staggered][cleanup] normal to lateral

parent 57427c3d
......@@ -136,10 +136,10 @@ private:
for (const auto& data : scvf.pairData())
{
auto& normalFace = fvGeometry.scvf(eIdx, data.localNormalFaceIdx);
if (!normalFace.boundary())
auto& lateralFace = fvGeometry.scvf(eIdx, data.localLateralFaceIdx);
if (!lateralFace.boundary())
{
const auto firstParallelElementDofIdx = normalFace.outsideScvIdx();
const auto firstParallelElementDofIdx = lateralFace.outsideScvIdx();
stencil.push_back(firstParallelElementDofIdx);
}
}
......@@ -159,9 +159,9 @@ private:
for (const auto& data : scvf.pairData())
{
// add normal dofs
stencil.push_back(data.normalPair.first);
stencil.push_back(data.lateralPair.first);
if (!scvf.boundary())
stencil.push_back(data.normalPair.second);
stencil.push_back(data.lateralPair.second);
// add parallel dofs
for (SmallLocalIndex i = 0; i < upwindSchemeOrder; i++)
......
......@@ -110,10 +110,10 @@ public:
const auto& subFaceData = scvf.pairData(i);
// treat the velocities normal to the face
velocityNormalInside_[i] = faceSol[subFaceData.normalPair.first];
velocityLateralInside_[i] = faceSol[subFaceData.lateralPair.first];
if (scvf.hasOuterNormal(i))
velocityNormalOutside_[i] = faceSol[subFaceData.normalPair.second];
if (scvf.hasOuterLateral(i))
velocityLateralOutside_[i] = faceSol[subFaceData.lateralPair.second];
// treat the velocities parallel to the self face
for (int j = 0; j < upwindSchemeOrder; j++)
......@@ -178,9 +178,9 @@ public:
*
* \param localSubFaceIdx The local index of the subface
*/
Scalar velocityNormalInside(const int localSubFaceIdx) const
Scalar velocityLateralInside(const int localSubFaceIdx) const
{
return velocityNormalInside_[localSubFaceIdx];
return velocityLateralInside_[localSubFaceIdx];
}
/*!
......@@ -188,9 +188,9 @@ public:
*
* \param localSubFaceIdx The local index of the subface
*/
Scalar velocityNormalOutside(const int localSubFaceIdx) const
Scalar velocityLateralOutside(const int localSubFaceIdx) const
{
return velocityNormalOutside_[localSubFaceIdx];
return velocityLateralOutside_[localSubFaceIdx];
}
private:
......@@ -221,8 +221,8 @@ private:
InAxisVelocities inAxisVelocities_;
std::array<std::array<Scalar, upwindSchemeOrder>, numPairs> velocityParallel_;
std::array<Scalar, numPairs> velocityNormalInside_;
std::array<Scalar, numPairs> velocityNormalOutside_;
std::array<Scalar, numPairs> velocityLateralInside_;
std::array<Scalar, numPairs> velocityLateralOutside_;
};
......
......@@ -53,11 +53,11 @@ struct PairData
std::bitset<upwindSchemeOrder> hasParallelNeighbor;
std::array<GridIndexType, upwindSchemeOrder> parallelDofs;
std::array<Scalar, upwindSchemeOrder> parallelCellWidths;
bool hasNormalNeighbor = false;
std::pair<GridIndexType, GridIndexType> normalPair;
SmallLocalIndexType localNormalFaceIdx;
Scalar normalDistance;
GlobalPosition virtualFirstParallelFaceDofPos;
bool hasOuterLateral = false;
std::pair<GridIndexType, GridIndexType> lateralPair;
SmallLocalIndexType localLateralFaceIdx;
Scalar lateralDistance;
GlobalPosition virtualBoundaryFaceDofPos;
};
......@@ -350,7 +350,7 @@ private:
}
/*!
* \brief Fills the pair data with the normal dofs and distances
* \brief Fills the pair data with the lateral dofs and distances
* and calls a further function to collect the parallel dofs and distances
*/
void fillPairData_()
......@@ -359,23 +359,23 @@ private:
pairData_ = {};
// set basic global positions
const auto& elementCenter = element_.geometry().center();
const auto& FacetCenter = intersection_.geometry().center();
const auto& selfElementCenter = element_.geometry().center();
const auto& selfFacetCenter = intersection_.geometry().center();
// get the inner normal Dof Index
SmallLocalIndexType numPairInnerNormalIdx = 0;
// get the inner lateral Dof Index
SmallLocalIndexType numPairInnerLateralIdx = 0;
for (const auto& innerElementIntersection : intersections(gridView_, element_))
{
if (facetIsNormal_(innerElementIntersection.indexInInside(), intersection_.indexInInside()))
{
const auto innerElementIntersectionIdx = innerElementIntersection.indexInInside();
setNormalPairFirstInfo_(innerElementIntersectionIdx, element_, numPairInnerNormalIdx);
numPairInnerNormalIdx++;
setLateralPairFirstInfo_(innerElementIntersectionIdx, element_, numPairInnerLateralIdx);
numPairInnerLateralIdx++;
}
}
// get the outer normal Dof Index
SmallLocalIndexType numPairOuterNormalIdx = 0;
// get the outer lateral Dof Index
SmallLocalIndexType numPairOuterLateralIdx = 0;
if (intersection_.neighbor())
{
// the direct neighbor element and the respective intersection index
......@@ -387,20 +387,25 @@ private:
if (facetIsNormal_(directNeighborElementIntersection.indexInInside(), intersection_.indexInOutside()))
{
const auto directNeighborElemIsIdx = directNeighborElementIntersection.indexInInside();
setNormalPairSecondInfo_(directNeighborElemIsIdx, directNeighborElement, numPairOuterNormalIdx);
numPairOuterNormalIdx++;
setLateralPairSecondInfo_(directNeighborElemIsIdx, directNeighborElement, numPairOuterLateralIdx);
numPairOuterLateralIdx++;
}
}
}
else // intersection is on boundary
{
// fill the normal pair entries
for(SmallLocalIndexType pairIdx = 0; pairIdx < numPairs; ++pairIdx)
for (const auto& intersection : intersections(gridView_, element_))
{
assert(!pairData_[pairIdx].hasNormalNeighbor);
const auto distance = FacetCenter - elementCenter;
pairData_[pairIdx].normalDistance = std::move(distance.two_norm());
numPairOuterNormalIdx++;
if (facetIsNormal_(intersection.indexInInside(), intersection_.indexInInside()))
{
assert(!pairData_[numPairOuterLateralIdx].hasOuterLateral);
const auto boundaryDistanceOffset = intersection.geometry().center() - selfElementCenter;
pairData_[numPairOuterLateralIdx].virtualBoundaryFaceDofPos = std::move(boundaryDistanceOffset + selfElementCenter);
const auto normalDistanceoffset = selfFacetCenter - selfElementCenter;
pairData_[numPairOuterLateralIdx].lateralDistance = std::move(normalDistanceoffset.two_norm());
numPairOuterLateralIdx++;
}
}
}
......@@ -427,7 +432,7 @@ private:
auto parallelAxisIdx = directionIndex(intersection);
if (intersection.neighbor())
{
// If the normal intersection has a neighboring cell, go in and store the parallel information.
// If the lateral intersection has a neighboring cell, go in and store the parallel information.
const auto& outerElement = intersection.outside();
pairData_[numPairParallelIdx].hasParallelNeighbor.set(0, true);
pairData_[numPairParallelIdx].parallelDofs[0] = gridView_.indexSet().subIndex(outerElement, parallelLocalIdx, codimIntersection);
......@@ -438,9 +443,9 @@ private:
// If the intersection has no neighbor we have to deal with the virtual outer parallel dof
const auto& boundaryFacetCenter = intersection.geometry().center();
const auto distance = boundaryFacetCenter - elementCenter;
const auto virtualFirstParallelFaceDofPos = intersection_.geometry().center() + distance;
const auto virtualBoundaryFaceDofPos = intersection_.geometry().center() + distance;
pairData_[numPairParallelIdx].virtualFirstParallelFaceDofPos = std::move(virtualFirstParallelFaceDofPos);
pairData_[numPairParallelIdx].virtualBoundaryFaceDofPos = std::move(virtualBoundaryFaceDofPos);
}
numPairParallelIdx++;
}
......@@ -468,19 +473,19 @@ private:
if( intersection.neighbor() )
{
auto parallelAxisIdx = directionIndex(intersection);
auto localNormalIntersectionIndex = intersection.indexInInside();
auto localLateralIntersectionIndex = intersection.indexInInside();
auto e = element_;
bool keepStacking = (parallelElementStack.size() < numParallelFaces);
while(keepStacking)
{
for(const auto& normalIntersection : intersections(gridView_, e))
for(const auto& lateralIntersection : intersections(gridView_, e))
{
if( normalIntersection.indexInInside() == localNormalIntersectionIndex )
if( lateralIntersection.indexInInside() == localLateralIntersectionIndex )
{
if( normalIntersection.neighbor() )
if( lateralIntersection.neighbor() )
{
parallelElementStack.push(normalIntersection.outside());
parallelElementStack.push(lateralIntersection.outside());
keepStacking = (parallelElementStack.size() < numParallelFaces);
}
else
......@@ -507,9 +512,9 @@ private:
const auto& boundaryFacetCenter = intersection.geometry().center();
const auto distance = boundaryFacetCenter - elementCenter;
const auto virtualFirstParallelFaceDofPos = intersection_.geometry().center() + distance;
const auto virtualBoundaryFaceDofPos = intersection_.geometry().center() + distance;
pairData_[numPairParallelIdx].virtualFirstParallelFaceDofPos = std::move(virtualFirstParallelFaceDofPos);
pairData_[numPairParallelIdx].virtualBoundaryFaceDofPos = std::move(virtualBoundaryFaceDofPos);
}
numPairParallelIdx++;
}
......@@ -542,29 +547,29 @@ private:
return element.template subEntity <1> (localFacetIdx);
};
//! Sets the information about the normal faces (within the element)
void setNormalPairFirstInfo_(const int isIdx, const Element& element, const int numPairsIdx)
//! Sets the information about the lateral faces (within the element)
void setLateralPairFirstInfo_(const int isIdx, const Element& element, const int numPairsIdx)
{
// store the inner normal dofIdx
pairData_[numPairsIdx].normalPair.first = gridView_.indexSet().subIndex(element, isIdx, codimIntersection);
// store the inner lateral dofIdx
pairData_[numPairsIdx].lateralPair.first = gridView_.indexSet().subIndex(element, isIdx, codimIntersection);
// store the local normal facet index
pairData_[numPairsIdx].localNormalFaceIdx = isIdx;
// store the local lateral facet index
pairData_[numPairsIdx].localLateralFaceIdx = isIdx;
}
//! Sets the information about the normal faces (in the neighbor element)
void setNormalPairSecondInfo_(const int isIdx, const Element& element, const int numPairsIdx)
//! Sets the information about the lateral faces (in the neighbor element)
void setLateralPairSecondInfo_(const int isIdx, const Element& element, const int numPairsIdx)
{
// store the dofIdx
pairData_[numPairsIdx].hasNormalNeighbor = true;
pairData_[numPairsIdx].normalPair.second = gridView_.indexSet().subIndex(element, isIdx, codimIntersection);
pairData_[numPairsIdx].hasOuterLateral = true;
pairData_[numPairsIdx].lateralPair.second = gridView_.indexSet().subIndex(element, isIdx, codimIntersection);
// store the element distance
const auto& outerNormalFacet = getFacet_(isIdx, element);
const auto outerNormalFacetPos = outerNormalFacet.geometry().center();
const auto& innerNormalFacet = getFacet_(isIdx, element_);
const auto innerNormalFacetPos = innerNormalFacet.geometry().center();
pairData_[numPairsIdx].normalDistance = (innerNormalFacetPos - outerNormalFacetPos).two_norm();
const auto& outerLateralFacet = getFacet_(isIdx, element);
const auto outerLateralFacetPos = outerLateralFacet.geometry().center();
const auto& innerLateralFacet = getFacet_(isIdx, element_);
const auto innerLateralFacetPos = innerLateralFacet.geometry().center();
pairData_[numPairsIdx].lateralDistance = (innerLateralFacetPos - outerLateralFacetPos).two_norm();
}
//! Sets the information about the parallel distances
......@@ -597,7 +602,7 @@ private:
const typename Element::Geometry elementGeometry_; //!< Reference to the element geometry
const GridView gridView_; //!< The grid view
AxisData axisData_; //!< Data related to forward and backward faces
std::array<PairData, numPairs> pairData_; //!< Collection of pair information related to normal and parallel faces
std::array<PairData, numPairs> pairData_; //!< Collection of pair information related to lateral and parallel faces
};
} // end namespace Dumux
......
......@@ -267,9 +267,9 @@ public:
*
* \param localSubFaceIdx The local index of the subface
*/
bool hasOuterNormal(const int localSubFaceIdx) const
bool hasOuterLateral(const int localSubFaceIdx) const
{
return pairData(localSubFaceIdx).hasNormalNeighbor;
return pairData(localSubFaceIdx).hasOuterLateral;
}
/*!
......@@ -330,7 +330,7 @@ public:
* \param localSubFaceIdx The local index of the subface
* \param parallelDegreeIdx The index describing how many faces away from the self
*/
Scalar cellCenteredParallelDistance(const int localSubFaceIdx, const int parallelDegreeIdx) const
Scalar parallelDofsDistance(const int localSubFaceIdx, const int parallelDegreeIdx) const
{
if (parallelDegreeIdx == 0)
return (faceLength(localSubFaceIdx) + pairData(localSubFaceIdx).parallelCellWidths[0]) * 0.5;
......
......@@ -117,38 +117,37 @@ public:
const FVElementGeometry& fvGeometry,
const Element& element,
const SubControlVolumeFace& scvf,
const SubControlVolumeFace& normalFace,
const ElementVolumeVariables& elemVolVars,
const FaceVariables& faceVars,
const GridFluxVariablesCache& gridFluxVarsCache,
const int localSubFaceIdx,
const bool lateralFaceHasDirichletPressure,
const bool lateralFaceHasBJS)
const Dune::Std::optional<BoundaryTypes>& lateralFaceBoundaryTypes)
{
const auto eIdx = scvf.insideScvIdx();
const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
// Get the transporting velocity, located at the scvf perpendicular to the current scvf where the dof
// of interest is located.
const Scalar transportingVelocity = faceVars.velocityNormalInside(localSubFaceIdx);
const Scalar transportingVelocity = faceVars.velocityLateralInside(localSubFaceIdx);
// Check whether the own or the neighboring element is upstream.
const bool selfIsUpstream = ( normalFace.directionSign() == sign(transportingVelocity) );
const bool selfIsUpstream = ( lateralFace.directionSign() == sign(transportingVelocity) );
const bool canHigherOrder = canLateralSecondOrder_(scvf, selfIsUpstream, localSubFaceIdx, std::integral_constant<bool, useHigherOrder>{});
if (canHigherOrder)
{
const auto parallelUpwindingMomenta = getLateralUpwindingMomenta_(problem, fvGeometry, element, scvf, elemVolVars, faceVars,
transportingVelocity, localSubFaceIdx,
lateralFaceHasDirichletPressure, lateralFaceHasBJS,
transportingVelocity, localSubFaceIdx, lateralFaceBoundaryTypes,
std::integral_constant<bool, useHigherOrder>{});
return doLateralMomentumUpwinding_(scvf, normalFace, parallelUpwindingMomenta, transportingVelocity, localSubFaceIdx, gridFluxVarsCache);
return doLateralMomentumUpwinding_(fvGeometry, scvf, parallelUpwindingMomenta, transportingVelocity, localSubFaceIdx, gridFluxVarsCache);
}
else
{
const auto parallelUpwindingMomenta = getLateralUpwindingMomenta_(problem, fvGeometry, element, scvf, elemVolVars, faceVars,
transportingVelocity, localSubFaceIdx,
lateralFaceHasDirichletPressure, lateralFaceHasBJS,
transportingVelocity, localSubFaceIdx, lateralFaceBoundaryTypes,
std::integral_constant<bool, false>{});
return doLateralMomentumUpwinding_(scvf, normalFace, parallelUpwindingMomenta, transportingVelocity, localSubFaceIdx, gridFluxVarsCache);
return doLateralMomentumUpwinding_(fvGeometry, scvf, parallelUpwindingMomenta, transportingVelocity, localSubFaceIdx, gridFluxVarsCache);
}
}
......@@ -370,11 +369,10 @@ private:
const FaceVariables& faceVars,
const Scalar transportingVelocity,
const int localSubFaceIdx,
bool lateralFaceHasDirichletPressure,
bool lateralFaceHasBJS,
const Dune::Std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
std::true_type)
{
const SubControlVolumeFace& lateralFace = fvGeometry.scvf(ownScvf.insideScvIdx(), ownScvf.pairData(localSubFaceIdx).localNormalFaceIdx);
const SubControlVolumeFace& lateralFace = fvGeometry.scvf(ownScvf.insideScvIdx(), ownScvf.pairData(localSubFaceIdx).localLateralFaceIdx);
// Get the volume variables of the own and the neighboring element
const auto& insideVolVars = elemVolVars[lateralFace.insideScvIdx()];
......@@ -384,10 +382,9 @@ private:
// thus we always use this value for the computation of the transported momentum.
if (!ownScvf.hasParallelNeighbor(localSubFaceIdx, 0))
{
const Scalar boundaryMomentum = getParallelVelocityFromBoundary_(problem, element, fvGeometry, ownScvf, lateralFace,
const Scalar boundaryMomentum = getParallelVelocityFromBoundary_(problem, element, fvGeometry, ownScvf,
faceVars.velocitySelf(), localSubFaceIdx,
lateralFaceHasDirichletPressure,
lateralFaceHasBJS) * insideVolVars.density();
lateralFaceBoundaryTypes) * insideVolVars.density();
return std::array<Scalar, 3>{boundaryMomentum, boundaryMomentum, boundaryMomentum};
}
......@@ -404,10 +401,9 @@ private:
if (ownScvf.hasParallelNeighbor(localSubFaceIdx, 0))
momenta[0] = faceVars.velocityParallel(localSubFaceIdx, 0) * insideVolVars.density();
else
momenta[0] = getParallelVelocityFromBoundary_(problem, element, fvGeometry, ownScvf, lateralFace,
momenta[0] = getParallelVelocityFromBoundary_(problem, element, fvGeometry, ownScvf,
faceVars.velocitySelf(), localSubFaceIdx,
lateralFaceHasDirichletPressure,
lateralFaceHasBJS) * insideVolVars.density();
lateralFaceBoundaryTypes) * insideVolVars.density();
// The local index of the faces that is opposite to localSubFaceIdx
const int oppositeSubFaceIdx = localSubFaceIdx % 2 ? localSubFaceIdx - 1 : localSubFaceIdx + 1;
......@@ -454,12 +450,11 @@ private:
const FaceVariables& faceVars,
const Scalar transportingVelocity,
const int localSubFaceIdx,
bool lateralFaceHasDirichletPressure,
bool lateralFaceHasBJS,
const Dune::Std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
std::false_type)
{
// Check whether the own or the neighboring element is upstream.
const SubControlVolumeFace& lateralFace = fvGeometry.scvf(ownScvf.insideScvIdx(), ownScvf.pairData(localSubFaceIdx).localNormalFaceIdx);
const SubControlVolumeFace& lateralFace = fvGeometry.scvf(ownScvf.insideScvIdx(), ownScvf.pairData(localSubFaceIdx).localLateralFaceIdx);
// Get the volume variables of the own and the neighboring element
const auto& insideVolVars = elemVolVars[lateralFace.insideScvIdx()];
......@@ -467,9 +462,9 @@ private:
const Scalar momentumParallel = ownScvf.hasParallelNeighbor(localSubFaceIdx, 0)
? faceVars.velocityParallel(localSubFaceIdx, 0) * outsideVolVars.density()
: (getParallelVelocityFromBoundary_(problem, element, fvGeometry, ownScvf, lateralFace,
faceVars.velocitySelf(), localSubFaceIdx,
lateralFaceHasDirichletPressure, lateralFaceHasBJS)
: (getParallelVelocityFromBoundary_(problem, element, fvGeometry, ownScvf,
faceVars.velocitySelf(), localSubFaceIdx,
lateralFaceBoundaryTypes)
* insideVolVars.density());
// If the lateral face lies on a boundary, we assume that the parallel velocity on the boundary is actually known,
......@@ -489,8 +484,8 @@ private:
*
* Fowards to Frontal Momentum Upwinding method
*/
static Scalar doLateralMomentumUpwinding_(const SubControlVolumeFace& scvf,
const SubControlVolumeFace& normalScvf,
static Scalar doLateralMomentumUpwinding_(const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const std::array<Scalar, 2>& momenta,
const Scalar transportingVelocity,
const int localSubFaceIdx,
......@@ -503,20 +498,22 @@ private:
* \brief Returns the upwinded momentum for higher order upwind schemes
*
* \param scvf The sub control volume face
* \param normalScvf The normal sub control volume face
* \param momenta The momenta to be upwinded
* \param transportingVelocity The average of the self and opposite velocities.
* \param localSubFaceIdx The local index of the subface
* \param gridFluxVarsCache The grid flux variables cache
*/
static Scalar doLateralMomentumUpwinding_(const SubControlVolumeFace& scvf,
const SubControlVolumeFace& normalScvf,
static Scalar doLateralMomentumUpwinding_(const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const std::array<Scalar, 3>& momenta,
const Scalar transportingVelocity,
const int localSubFaceIdx,
const GridFluxVariablesCache& gridFluxVarsCache)
{
const bool selfIsUpstream = ( normalScvf.directionSign() == sign(transportingVelocity) );
const auto eIdx = scvf.insideScvIdx();
const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
const bool selfIsUpstream = ( lateralFace.directionSign() == sign(transportingVelocity) );
const auto& upwindScheme = gridFluxVarsCache.staggeredUpwindMethods();
std::array<Scalar,3> distances = getLateralDistances_(scvf, localSubFaceIdx, selfIsUpstream);
return upwindScheme.tvd(momenta, distances, selfIsUpstream, upwindScheme.tvdApproach());
......@@ -542,8 +539,8 @@ private:
// The local index of the faces that is opposite to localSubFaceIdx
const int oppositeSubFaceIdx = localSubFaceIdx % 2 ? localSubFaceIdx - 1 : localSubFaceIdx + 1;
distances[0] = ownScvf.cellCenteredParallelDistance(localSubFaceIdx, 0);
distances[1] = ownScvf.cellCenteredParallelDistance(oppositeSubFaceIdx, 0);
distances[0] = ownScvf.parallelDofsDistance(localSubFaceIdx, 0);
distances[1] = ownScvf.parallelDofsDistance(oppositeSubFaceIdx, 0);
if (ownScvf.hasParallelNeighbor(localSubFaceIdx, 0))
distances[2] = ownScvf.pairData(localSubFaceIdx).parallelCellWidths[0];
else
......@@ -551,8 +548,8 @@ private:
}
else
{
distances[0] = ownScvf.cellCenteredParallelDistance(localSubFaceIdx, 0);
distances[1] = ownScvf.cellCenteredParallelDistance(localSubFaceIdx, 1);
distances[0] = ownScvf.parallelDofsDistance(localSubFaceIdx, 0);
distances[1] = ownScvf.parallelDofsDistance(localSubFaceIdx, 1);
distances[2] = ownScvf.pairData(localSubFaceIdx).parallelCellWidths[0];
}
......@@ -566,23 +563,23 @@ private:
*
* \param problem The problem
* \param scvf The SubControlVolumeFace that is normal to the boundary
* \param normalFace The face at the boundary
* \param velocitySelf the velocity at scvf
* \param localSubFaceIdx The local index of the face that is on the boundary
* \param element The element that is on the boundary
* \param lateralFaceHasDirichletPressure @c true if there is a dirichlet condition for the pressure on the boundary
* \param lateralFaceHasBJS @c true if there is a BJS condition fot the velocity on the boundary
* \param lateralFaceBoundaryTypes stores the type of boundary at the lateral face
*/
static Scalar getParallelVelocityFromBoundary_(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const SubControlVolumeFace& normalFace,
const Scalar velocitySelf,
const int localSubFaceIdx,
const bool lateralFaceHasDirichletPressure,
const bool lateralFaceHasBJS)
const Dune::Std::optional<BoundaryTypes>& lateralFaceBoundaryTypes)
{
// Find out what boundary type is set on the lateral face
const bool lateralFaceHasDirichletPressure = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx);
const bool lateralFaceHasBJS = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isBJS(Indices::velocity(scvf.directionIndex()));
// If there is a Dirichlet condition for the pressure we assume zero gradient for the velocity,
// so the velocity at the boundary equal to that on the scvf.
if (lateralFaceHasDirichletPressure)
......@@ -590,8 +587,11 @@ private:
if (lateralFaceHasBJS)
{
const auto eIdx = scvf.insideScvIdx();
const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
const SubControlVolume& scv = fvGeometry.scv(scvf.insideScvIdx());
return problem.bjsVelocity(element, scv, normalFace, velocitySelf);
return problem.bjsVelocity(element, scv, lateralFace, velocitySelf);
}
else
{
......@@ -626,31 +626,31 @@ private:
const Scalar parallelVelocity)
{
// A ghost subface at the boundary is created, featuring the location of the sub face's center
const SubControlVolumeFace& boundaryNormalFace = fvGeometry.scvf(scvf.insideScvIdx(), scvf.pairData(localIdx).localNormalFaceIdx);
GlobalPosition boundarySubFaceCenter = scvf.pairData(localIdx).virtualFirstParallelFaceDofPos + boundaryNormalFace.center();
boundarySubFaceCenter *= 0.5;
const SubControlVolumeFace& lateralFace = fvGeometry.scvf(scvf.insideScvIdx(), scvf.pairData(localIdx).localLateralFaceIdx);
GlobalPosition lateralBoundaryFaceCenter = scvf.pairData(localIdx).virtualBoundaryFaceDofPos + lateralFace.center();
lateralBoundaryFaceCenter *= 0.5;
// ________________
// --------####o### || frontal face of staggered half-control-volume
// | || | current scvf # localSubFace of interest, lies on lateral boundary
// | || | current scvf # lateralBoundaryFace of interest, lies on lateral boundary
// | || | x dof position
// | || x~~~~> vel.Self -- element boundaries
// | || | __ domain boundary
// | || | o position at which the boundary conditions will be evaluated
// ---------------- (boundarySubFaceCenter)
// ---------------- (lateralBoundaryFaceCenter)
const SubControlVolumeFace boundarySubFace = boundaryNormalFace.makeBoundaryFace(boundarySubFaceCenter);
const SubControlVolumeFace lateralBoundaryFace = lateralFace.makeBoundaryFace(lateralBoundaryFaceCenter);
// The boundary condition is checked, in case of symmetry or Dirichlet for the pressure
// a gradient of zero is assumed in the direction normal to the bounadry, while if there is
// Dirichlet of BJS for the velocity the related values are exploited.
const auto bcTypes = problem.boundaryTypes(boundaryElement, boundarySubFace);
const auto bcTypes = problem.boundaryTypes(boundaryElement, lateralBoundaryFace);
if (bcTypes.isDirichlet(Indices::velocity(scvf.directionIndex())))
{
// ________________
// --------#######o || frontal face of staggered half-control-volume
// | || | current scvf # localSubFace of interest, lies on lateral boundary
// | || | current scvf # lateralBoundaryFace of interest, lies on lateral boundary
// | || | x dof position
// | || x~~~~> vel.Self -- element boundaries
// | || | __ domain boundary
......@@ -665,19 +665,19 @@ private:
else if (bcTypes.isBJS(Indices::velocity(scvf.directionIndex())))
{
const SubControlVolume& scv = fvGeometry.scv(scvf.insideScvIdx());
return problem.bjsVelocity(boundaryElement, scv,boundaryNormalFace, parallelVelocity);
return problem.bjsVelocity(boundaryElement, scv, lateralFace, parallelVelocity);
}
else
{
// Neumann conditions are not well implemented
DUNE_THROW(Dune::InvalidStateException, "Something went wrong with the boundary conditions for the momentum equations at global position " << boundarySubFaceCenter);
DUNE_THROW(Dune::InvalidStateException, "Something went wrong with the boundary conditions for the momentum equations at global position " << lateralBoundaryFaceCenter);
}
}