Merge branch 'feature/improve-staggered-scvf' into 'master'

[staggered][scvf] Replace virtualBoundaryFaceDofPos by laterFaceCenter

See merge request !1631

dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh

@@ -57,7 +57,7 @@ struct PairData
     std::pair<GridIndexType, GridIndexType> lateralPair;
     SmallLocalIndexType localLateralFaceIdx;
     Scalar lateralDistance;
-    GlobalPosition virtualBoundaryFaceDofPos;
+    GlobalPosition lateralStaggeredFaceCenter;
 };
 
 
@@ -370,6 +370,10 @@ private:
             {
                 const auto innerElementIntersectionIdx = innerElementIntersection.indexInInside();
                 setLateralPairFirstInfo_(innerElementIntersectionIdx, element_, numPairInnerLateralIdx);
+
+                const auto distance = innerElementIntersection.geometry().center() - selfElementCenter;
+                pairData_[numPairInnerLateralIdx].lateralStaggeredFaceCenter = std::move(selfFacetCenter + distance);
+
                 numPairInnerLateralIdx++;
             }
         }
@@ -399,8 +403,6 @@ private:
                 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());
@@ -419,8 +421,6 @@ private:
     void fillParallelPairData_(std::false_type)
     {
         // set basic global positions and stencil size definitions
-        const auto& elementCenter = element_.geometry().center();
-
         // get the parallel Dofs
         const auto parallelLocalIdx = intersection_.indexInInside();
         SmallLocalIndexType numPairParallelIdx = 0;
@@ -438,15 +438,7 @@ private:
                     pairData_[numPairParallelIdx].parallelDofs[0] = gridView_.indexSet().subIndex(outerElement, parallelLocalIdx, codimIntersection);
                     pairData_[numPairParallelIdx].parallelCellWidths[0] = setParallelPairCellWidths_(outerElement, parallelAxisIdx);
                 }
-                else  // No parallel neighbor available
-                {
-                    // 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 virtualBoundaryFaceDofPos = intersection_.geometry().center() + distance;
 
-                    pairData_[numPairParallelIdx].virtualBoundaryFaceDofPos = std::move(virtualBoundaryFaceDofPos);
-                }
                 numPairParallelIdx++;
             }
         }
@@ -460,7 +452,6 @@ private:
     {
         // set basic global positions and stencil size definitions
         const auto numParallelFaces = pairData_[0].parallelCellWidths.size();
-        const auto& elementCenter = element_.geometry().center();
 
         // get the parallel Dofs
         const auto parallelLocalIdx = intersection_.indexInInside();
@@ -506,16 +497,7 @@ private:
                     }
 
                 }
-                else
-                {
-                    // 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 virtualBoundaryFaceDofPos = intersection_.geometry().center() + distance;
-
-                    pairData_[numPairParallelIdx].virtualBoundaryFaceDofPos = std::move(virtualBoundaryFaceDofPos);
-                }
                 numPairParallelIdx++;
             }
         }

dumux/freeflow/navierstokes/staggered/fluxvariables.hh

@@ -287,7 +287,7 @@ public:
         {
             const auto eIdx = scvf.insideScvIdx();
             // Get the face normal to the face the dof lives on. The staggered sub face conincides with half of this lateral face.
-            const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
+            const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
 
             // create a boundaryTypes object (will be empty if not at a boundary)
             Dune::Std::optional<BoundaryTypes> lateralFaceBoundaryTypes;
@@ -298,7 +298,7 @@ public:
             {
                 // Construct a temporary scvf which corresponds to the staggered sub face, featuring the location
                 // the sub faces's center.
-                auto lateralBoundaryFaceCenter = scvf.pairData(localSubFaceIdx).virtualBoundaryFaceDofPos + lateralFace.center();
+                auto lateralBoundaryFaceCenter = scvf.pairData(localSubFaceIdx).lateralStaggeredFaceCenter + lateralScvf.center();
                 lateralBoundaryFaceCenter *= 0.5;
 
                 //    ________________
@@ -310,7 +310,7 @@ public:
                 //    |      ||      |                 o  position at which the boundary conditions will be evaluated
                 //    ----------------                    (lateralBoundaryFaceCenter)
 
-                const auto lateralBoundaryFace = lateralFace.makeBoundaryFace(lateralBoundaryFaceCenter);
+                const auto lateralBoundaryFace = lateralScvf.makeBoundaryFace(lateralBoundaryFaceCenter);
 
                 // Retrieve the boundary types that correspond to the sub face.
                 lateralFaceBoundaryTypes.emplace(problem.boundaryTypes(element, lateralBoundaryFace));
@@ -324,7 +324,7 @@ public:
                 if (lateralFaceBoundaryTypes->isNeumann(Indices::velocity(scvf.directionIndex())))
                 {
                     lateralFlux += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, lateralBoundaryFace)[Indices::velocity(scvf.directionIndex())]
-                                                  * elemVolVars[lateralFace.insideScvIdx()].extrusionFactor() * lateralFace.area() * 0.5;
+                                                  * elemVolVars[lateralScvf.insideScvIdx()].extrusionFactor() * lateralScvf.area() * 0.5;
                     continue;
                 }
 
@@ -457,7 +457,7 @@ private:
         {
             // Create a boundaryTypes object (will be empty if not at a boundary).
             Dune::Std::optional<BoundaryTypes> bcTypes;
-            const auto& boundaryFace = scvf.makeBoundaryFace(scvf.pairData(localSubFaceIdx).virtualBoundaryFaceDofPos);
+            const auto& boundaryFace = scvf.makeBoundaryFace(scvf.pairData(localSubFaceIdx).lateralStaggeredFaceCenter);
 
             // Get the boundary conditions if we are at a boundary.
             if (scvf.boundary())
@@ -612,7 +612,7 @@ private:
         //    |      ||      |                 o  position at which the boundary conditions will be evaluated
         //    ----------------
 
-        return ownScvf.makeBoundaryFace(ownScvf.pairData(localSubFaceIdx).virtualBoundaryFaceDofPos);
+        return ownScvf.makeBoundaryFace(ownScvf.pairData(localSubFaceIdx).lateralStaggeredFaceCenter);
     };
 
     //! helper function to get the averaged extrusion factor for a face

dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh

@@ -426,7 +426,7 @@ private:
                 momenta[2] = faceVars.velocityParallel(localSubFaceIdx, 1) * outsideVolVars.density();
             else
             {
-                const Element& elementParallel = fvGeometry.fvGridGeometry().element(fvGeometry.scv(lateralFace.outsideScvIdx()));
+                const Element& elementParallel = fvGeometry.fvGridGeometry().element(lateralFace.outsideScvIdx());
                 const SubControlVolumeFace& firstParallelScvf = fvGeometry.scvf(lateralFace.outsideScvIdx(), ownScvf.localFaceIdx());
                 momenta[2] = getParallelVelocityFromOtherBoundary_(problem, fvGeometry, firstParallelScvf,
                                                                    localSubFaceIdx, elementParallel,
@@ -626,8 +626,8 @@ private:
                                                         const Scalar parallelVelocity)
     {
         // A ghost subface at the boundary is created, featuring the location of the sub face's center
-        const SubControlVolumeFace& lateralFace = fvGeometry.scvf(scvf.insideScvIdx(), scvf.pairData(localIdx).localLateralFaceIdx);
-        GlobalPosition lateralBoundaryFaceCenter = scvf.pairData(localIdx).virtualBoundaryFaceDofPos + lateralFace.center();
+        const SubControlVolumeFace& lateralScvf = fvGeometry.scvf(scvf.insideScvIdx(), scvf.pairData(localIdx).localLateralFaceIdx);
+        GlobalPosition lateralBoundaryFaceCenter = scvf.pairData(localIdx).lateralStaggeredFaceCenter + lateralScvf.center();
         lateralBoundaryFaceCenter *= 0.5;
 
         //    ________________
@@ -639,7 +639,7 @@ private:
         //    |      ||      |                 o  position at which the boundary conditions will be evaluated
         //    ----------------                    (lateralBoundaryFaceCenter)
 
-        const SubControlVolumeFace lateralBoundaryFace = lateralFace.makeBoundaryFace(lateralBoundaryFaceCenter);
+        const SubControlVolumeFace lateralBoundaryFace = lateralScvf.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
@@ -665,7 +665,7 @@ private:
         else if (bcTypes.isBJS(Indices::velocity(scvf.directionIndex())))
         {
             const SubControlVolume& scv = fvGeometry.scv(scvf.insideScvIdx());
-            return problem.bjsVelocity(boundaryElement, scv, lateralFace, parallelVelocity);
+            return problem.bjsVelocity(boundaryElement, scv, lateralScvf, parallelVelocity);
         }
         else
         {
@@ -677,7 +677,7 @@ private:
     //! helper function to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
     static SubControlVolumeFace makeParallelGhostFace_(const SubControlVolumeFace& ownScvf, const int localSubFaceIdx)
     {
-        return ownScvf.makeBoundaryFace(ownScvf.pairData(localSubFaceIdx).virtualBoundaryFaceDofPos);
+        return ownScvf.makeBoundaryFace(ownScvf.pairData(localSubFaceIdx).lateralStaggeredFaceCenter);
     };
 };