[staggeredGrid] Add direction index for faces and fix flux calculation

dumux/discretization/staggered/staggeredgeometryhelper.hh

@@ -29,6 +29,7 @@
 
 #include <dumux/common/math.hh>
 #include <type_traits>
+#include <algorithm>
 
 namespace Dumux
 {
@@ -122,6 +123,17 @@ public:
         return pairData_;
     }
 
+     /*!
+     * \brief Returns the dirction index of the facet (0 = x, 1 = y, 2 = z)
+     */
+    int directionIndex() const
+    {
+        const Scalar eps = 1e-8;
+        const auto& direction = intersection_.centerUnitOuterNormal();
+        const int idx = std::find_if(direction.begin(), direction.end(), [eps](const auto& x) { return std::abs(x) > eps; } ) - direction.begin();
+        return idx;
+    }
+
 private:
      /*!
      * \brief Fills all entries of the pair data

dumux/discretization/staggered/subcontrolvolumeface.hh

@@ -104,6 +104,7 @@ public:
 
           pairData_ = geometryHelper.pairData();
           localFaceIdx_ = is.indexInInside();
+          dirIdx_ = geometryHelper.directionIndex();
       }
 
     /*//! The copy constrcutor
@@ -209,6 +210,12 @@ public:
         return localFaceIdx_;
     }
 
+    //! Returns the dirction index of the facet (0 = x, 1 = y, 2 = z)
+    int directionIndex() const
+    {
+        return dirIdx_;
+    }
+
     //! The global index of this sub control volume face
     Scalar selfToOppositeDistance() const
     {
@@ -242,6 +249,7 @@ private:
     std::vector<StaggeredSubFace> subfaces_;
     std::array<PairData<Scalar, GlobalPosition>, numPairs> pairData_;
     int localFaceIdx_;
+    int dirIdx_;
 
 };
 

dumux/freeflow/staggered/localresidual.hh

@@ -102,21 +102,19 @@ public:
 
         CellCenterPrimaryVariables flux(0.0);
 
-        const Scalar eps = 1e-6;
-
-        for(int direction = 0; direction < dim; ++direction)
+        if(scvf.unitOuterNormal()[scvf.directionIndex()] > 0.0) // positive coordinate direction
         {
-            if(scvf.unitOuterNormal()[direction] > eps && velocity > 0)
-            {
+            if(velocity > 0.0)
                 flux[0] = insideVolVars.density(0) * velocity;
-                return flux;
-            }
-
-            if(scvf.unitOuterNormal()[direction] < eps && velocity < 0)
-            {
+            else
+                flux[0] = outsideVolVars.density(0) * velocity;
+        }
+        else // negative coordinate direction
+        {
+            if(velocity > 0.0)
                 flux[0] = outsideVolVars.density(0) * velocity;
-                return flux;
-            }
+            else
+                flux[0] = insideVolVars.density(0) * velocity;
         }
         return flux;
     }