Merge branch 'feature/ff-get-boundary-vel-in-fluxvars' into 'master'

[staggered][faceVariables] Retrive boundary values in fluxVars

See merge request !808

dumux/discretization/staggered/freeflow/facevariables.hh

@@ -86,24 +86,6 @@ public:
         velocitySelf_ = faceSol[scvf.dofIndex()];
         velocityOpposite_ = faceSol[scvf.dofIndexOpposingFace()];
 
-        // lambda to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
-        auto makeGhostFace = [](const auto& pos)
-        {
-            return SubControlVolumeFace(pos, std::vector<unsigned int>{0,0});
-        };
-
-        // lambda to check whether there is a parallel face neighbor
-        auto hasParallelNeighbor = [](const auto& subFaceData)
-        {
-            return subFaceData.outerParallelFaceDofIdx >= 0;
-        };
-
-        // lambda to check whether there is a normal face neighbor
-        auto hasNormalNeighbor = [](const auto& subFaceData)
-        {
-            return subFaceData.normalPair.second >= 0;
-        };
-
         // handle all sub faces
         for(int i = 0; i < scvf.pairData().size(); ++i)
         {
@@ -112,21 +94,12 @@ public:
             // treat the velocities normal to the face
             velocityNormalInside_[i] = faceSol[subFaceData.normalPair.first];
 
-            if(hasNormalNeighbor(subFaceData))
-            {
+            if(scvf.hasFrontalNeighbor(i))
                 velocityNormalOutside_[i] = faceSol[subFaceData.normalPair.second];
-            }
-            else
-            {
-                const auto& normalFace = fvGeometry.scvf(scvf.insideScvIdx(), subFaceData.localNormalFaceIdx);
-                const auto normalDirIdx = normalFace.directionIndex();
-                velocityNormalOutside_[i] = problem.dirichlet(element, makeGhostFace(subFaceData.virtualOuterNormalFaceDofPos))[Indices::velocity(normalDirIdx)];
-            }
 
             // treat the velocity parallel to the face
-            velocityParallel_[i] = hasParallelNeighbor(subFaceData) ?
-                                   velocityParallel_[i] = faceSol[subFaceData.outerParallelFaceDofIdx] :
-                                   problem.dirichlet(element, makeGhostFace(subFaceData.virtualOuterParallelFaceDofPos))[Indices::velocity(scvf.directionIndex())];
+            if(scvf.hasParallelNeighbor(i))
+                velocityParallel_[i] = faceSol[subFaceData.outerParallelFaceDofIdx];
         }
     }
 

dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh

@@ -31,7 +31,6 @@
 #include <dumux/discretization/staggered/subcontrolvolumeface.hh>
 #include <dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh>
 #include <dumux/common/properties.hh>
-#include <dumux/common/optional.hh>
 
 #include <typeinfo>
 
@@ -82,34 +81,32 @@ public:
       unitOuterNormal_(is.centerUnitOuterNormal()),
       scvfIndex_(scvfIndex),
       scvIndices_(scvIndices),
-      boundary_(is.boundary())
+      boundary_(is.boundary()),
+      dofIdx_(geometryHelper.dofIndex()),
+      dofIdxOpposingFace_(geometryHelper.dofIndexOpposingFace()),
+      selfToOppositeDistance_(geometryHelper.selfToOppositeDistance()),
+      pairData_(geometryHelper.pairData()),
+      localFaceIdx_(geometryHelper.localFaceIndex()),
+      dirIdx_(geometryHelper.directionIndex()),
+      normalInPosCoordDir_(unitOuterNormal_[directionIndex()] > 0.0),
+      outerNormalScalar_(unitOuterNormal_[directionIndex()]),
+      isGhostFace_(false)
       {
           corners_.resize(isGeometry.corners());
           for (int i = 0; i < isGeometry.corners(); ++i)
               corners_[i] = isGeometry.corner(i);
-
-          dofIdx_ = geometryHelper.dofIndex();
-          oppositeIdx_ = geometryHelper.dofIndexOpposingFace();
-          selfToOppositeDistance_ = geometryHelper.selfToOppositeDistance();
-
-          pairData_ = geometryHelper.pairData();
-          localFaceIdx_ = geometryHelper.localFaceIndex();
-          dirIdx_ = geometryHelper.directionIndex();
-          normalInPosCoordDir_ = unitOuterNormal()[directionIndex()] > 0.0;
-          outerNormalScalar_ = unitOuterNormal()[directionIndex()];
-          isGhostFace_ = false;
       }
 
-      //! Constructor for a ghost face outside of the domain. Only needed to retrieve the center and scvIndices
-      FreeFlowStaggeredSubControlVolumeFace(const GlobalPosition& dofPosition,
-                                    const std::vector<GridIndexType>& scvIndices)
-      {
-          isGhostFace_ = true;
-          center_ = dofPosition;
-          scvIndices_ = scvIndices;
-          scvfIndex_ = -1;
-          dofIdx_ = -1;
-      }
+    //! Constructor for a ghost face outside of the domain. Only needed to retrieve the center and scvIndices
+    FreeFlowStaggeredSubControlVolumeFace(const GlobalPosition& dofPosition,
+                                          const std::vector<GridIndexType>& scvIndices,
+                                          const int dofIdx = -1)
+    : center_(dofPosition),
+      scvfIndex_(-1),
+      scvIndices_(scvIndices),
+      dofIdx_(dofIdx),
+      isGhostFace_(true)
+    {}
 
     //! The center of the sub control volume face
     const GlobalPosition& center() const
@@ -189,7 +186,7 @@ public:
     //! The global index of the dof living on the opposing face
     GridIndexType dofIndexOpposingFace() const
     {
-        return oppositeIdx_;
+        return dofIdxOpposingFace_;
     }
 
     //! The local index of this sub control volume face
@@ -234,6 +231,21 @@ public:
         return pairData_;
     }
 
+    bool isGhostFace() const
+    {
+        return isGhostFace_;
+    }
+
+    bool hasParallelNeighbor(const int localFaceIdx) const
+    {
+        return pairData_[localFaceIdx].outerParallelFaceDofIdx >= 0;
+    }
+
+    bool hasFrontalNeighbor(const int localFaceIdx) const
+    {
+        return pairData_[localFaceIdx].normalPair.second >= 0;
+    }
+
 private:
     Dune::GeometryType geomType_;
     std::vector<GlobalPosition> corners_;
@@ -245,7 +257,7 @@ private:
     bool boundary_;
 
     int dofIdx_;
-    int oppositeIdx_;
+    int dofIdxOpposingFace_;
     Scalar selfToOppositeDistance_;
     std::array<PairData<Scalar, GlobalPosition>, numPairs> pairData_;
     int localFaceIdx_;

dumux/freeflow/navierstokes/staggered/fluxvariables.hh

@@ -283,16 +283,10 @@ public:
             const auto& normalFace = fvGeometry.scvf(eIdx, scvf.pairData()[localSubFaceIdx].localNormalFaceIdx);
 
             // Check if we have a symmetry boundary condition. If yes, the tangental part of the momentum flux can be neglected.
-            if(scvf.pairData()[localSubFaceIdx].outerParallelFaceDofIdx < 0)
+            if(!scvf.hasParallelNeighbor(localSubFaceIdx))
             {
-                // Lambda to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest.
-                auto makeGhostFace = [eIdx] (const GlobalPosition& pos)
-                {
-                    return SubControlVolumeFace(pos, std::vector<unsigned int>{eIdx,eIdx});
-                };
-
                 // Use the ghost face to check if there is a symmetry boundary condition and skip any further steps if yes.
-                const auto bcTypes = problem.boundaryTypes(element, makeGhostFace(scvf.pairData()[localSubFaceIdx].virtualOuterParallelFaceDofPos));
+                const auto bcTypes = problem.boundaryTypes(element, makeParallelGhostFace_(scvf, localSubFaceIdx));
                 if(bcTypes.isSymmetry())
                     continue;
             }
@@ -346,7 +340,18 @@ private:
 
         // Get the velocities at the current (own) scvf and at the parallel one at the neighboring scvf.
         const Scalar velocitySelf = faceVars.velocitySelf();
-        const Scalar velocityParallel = faceVars.velocityParallel(localSubFaceIdx);
+
+        // Lambda to conveniently get the outer parallel velocity for normal faces that are on the boundary
+        // and therefore have no neighbor. Calls the problem to retrieve a fixed value set on the boundary.
+        auto getParallelVelocityFromBoundary = [&]()
+        {
+            const auto ghostFace = makeParallelGhostFace_(scvf, localSubFaceIdx);
+            return problem.dirichlet(element, ghostFace)[Indices::velocity(scvf.directionIndex())];
+        };
+
+        const Scalar velocityParallel = scvf.hasParallelNeighbor(localSubFaceIdx) ?
+                                        faceVars.velocityParallel(localSubFaceIdx)
+                                      : getParallelVelocityFromBoundary();
 
         // Get the volume variables of the own and the neighboring element
         const auto& insideVolVars = elemVolVars[normalFace.insideScvIdx()];
@@ -420,7 +425,19 @@ private:
         // the outer one at the respective staggered face of the element on the other side of the
         // current scvf.
         const Scalar innerNormalVelocity = faceVars.velocityNormalInside(localSubFaceIdx);
-        const Scalar outerNormalVelocity = faceVars.velocityNormalOutside(localSubFaceIdx);
+
+        // Lambda to conveniently get the outer normal velocity for faces that are on the boundary
+        // and therefore have no neighbor. Calls the problem to retrieve a fixed value set on the boundary.
+        auto getNormalVelocityFromBoundary = [&]()
+        {
+            const auto ghostFace = makeNormalGhostFace_(scvf, localSubFaceIdx);
+            return problem.dirichlet(element, ghostFace)[Indices::velocity(normalFace.directionIndex())];
+        };
+
+        const Scalar outerNormalVelocity = scvf.hasFrontalNeighbor(localSubFaceIdx) ?
+                                           faceVars.velocityNormalOutside(localSubFaceIdx)
+                                         : getNormalVelocityFromBoundary();
+
 
         // Calculate the velocity gradient in positive coordinate direction.
         const Scalar normalDeltaV = scvf.normalInPosCoordDir() ?
@@ -435,7 +452,18 @@ private:
         // For the parallel derivative, get the velocities at the current (own) scvf
         // and at the parallel one at the neighboring scvf.
         const Scalar innerParallelVelocity = faceVars.velocitySelf();
-        const Scalar outerParallelVelocity = faceVars.velocityParallel(localSubFaceIdx);
+
+        // Lambda to conveniently get the outer parallel velocity for normal faces that are on the boundary
+        // and therefore have no neighbor. Calls the problem to retrieve a fixed value set on the boundary.
+        auto getParallelVelocityFromBoundary = [&]()
+        {
+            const auto ghostFace = makeParallelGhostFace_(scvf, localSubFaceIdx);
+            return problem.dirichlet(element, ghostFace)[Indices::velocity(scvf.directionIndex())];
+        };
+
+        const Scalar outerParallelVelocity = scvf.hasParallelNeighbor(localSubFaceIdx) ?
+                                             faceVars.velocityParallel(localSubFaceIdx)
+                                           : getParallelVelocityFromBoundary();
 
         // The velocity gradient already accounts for the orientation
         // of the staggered face's outer normal vector.
@@ -494,6 +522,27 @@ private:
         // Account for the orientation of the face at the boundary,
         return outflow * scvf.directionSign();
     }
+
+private:
+
+    //! helper functiuob to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
+    SubControlVolumeFace makeGhostFace_(const SubControlVolumeFace& ownScvf, const GlobalPosition& pos) const
+    {
+        return SubControlVolumeFace(pos, std::vector<unsigned int>{ownScvf.insideScvIdx(), ownScvf.outsideScvIdx()}, ownScvf.dofIndex());
+    };
+
+    //! helper functiuob to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
+    SubControlVolumeFace makeNormalGhostFace_(const SubControlVolumeFace& ownScvf, const int localSubFaceIdx) const
+    {
+        return makeGhostFace_(ownScvf, ownScvf.pairData(localSubFaceIdx).virtualOuterNormalFaceDofPos);
+    };
+
+    //! helper functiuob to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
+    SubControlVolumeFace makeParallelGhostFace_(const SubControlVolumeFace& ownScvf, const int localSubFaceIdx) const
+    {
+        return makeGhostFace_(ownScvf, ownScvf.pairData(localSubFaceIdx).virtualOuterParallelFaceDofPos);
+    };
+
 };
 
 } // end namespace