From 13ebe79ddf89945235e84aefe1f77c8052dd63bb Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Mon, 5 Dec 2016 09:09:34 +0100
Subject: [PATCH] [staggeredGrid] Add direction index for faces and fix flux
 calculation

---
 .../staggered/staggeredgeometryhelper.hh      | 12 ++++++++++
 .../staggered/subcontrolvolumeface.hh         |  8 +++++++
 dumux/freeflow/staggered/localresidual.hh     | 22 +++++++++----------
 3 files changed, 30 insertions(+), 12 deletions(-)

diff --git a/dumux/discretization/staggered/staggeredgeometryhelper.hh b/dumux/discretization/staggered/staggeredgeometryhelper.hh
index 236fa8b9ff..0eb8233df9 100644
--- a/dumux/discretization/staggered/staggeredgeometryhelper.hh
+++ b/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
diff --git a/dumux/discretization/staggered/subcontrolvolumeface.hh b/dumux/discretization/staggered/subcontrolvolumeface.hh
index d513bc9e45..8d2c86023b 100644
--- a/dumux/discretization/staggered/subcontrolvolumeface.hh
+++ b/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_;
 
 };
 
diff --git a/dumux/freeflow/staggered/localresidual.hh b/dumux/freeflow/staggered/localresidual.hh
index 0d03535433..1a11198584 100644
--- a/dumux/freeflow/staggered/localresidual.hh
+++ b/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;
     }
-- 
GitLab