From af52a30e43a43d4e59274666ef5d2de7d8de9411 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Tue, 24 Jul 2018 17:23:29 +0200
Subject: [PATCH] [rans] Check isOnWall based on scvf, not on pos

* internally forward to isOnWallAtPos
* needed for coupling
---
 dumux/freeflow/rans/problem.hh                | 33 ++++++++++++-------
 dumux/freeflow/rans/twoeq/kepsilon/problem.hh |  2 +-
 .../twoeq/kepsilon/staggered/fluxvariables.hh |  2 +-
 test/freeflow/rans/pipelauferproblem.hh       | 14 ++++----
 test/freeflow/ransnc/flatplatetestproblem.hh  | 12 +++----
 5 files changed, 36 insertions(+), 27 deletions(-)

diff --git a/dumux/freeflow/rans/problem.hh b/dumux/freeflow/rans/problem.hh
index 0605d693a0..a886af7141 100644
--- a/dumux/freeflow/rans/problem.hh
+++ b/dumux/freeflow/rans/problem.hh
@@ -110,25 +110,24 @@ public:
         std::vector<unsigned int> wallElements;
         std::vector<GlobalPosition> wallPositions;
         std::vector<unsigned int> wallNormalAxisTemp;
-        auto& gridView(this->fvGridGeometry().gridView());
+
+        const auto gridView = this->fvGridGeometry().gridView();
+        auto fvGeometry = localView(this->fvGridGeometry());
+
         for (const auto& element : elements(gridView))
         {
-            for (const auto& intersection : intersections(gridView, element))
+            fvGeometry.bindElement(element);
+            for (const auto& scvf : scvfs(fvGeometry))
             {
                 // only search for walls at a global boundary
-                if (!intersection.boundary())
+                if (!scvf.boundary())
                     continue;
 
-                GlobalPosition global = intersection.geometry().center();
-                if (asImp_().isOnWall(global))
+                if (asImp_().isOnWall(scvf))
                 {
                     wallElements.push_back(this->fvGridGeometry().elementMapper().index(element));
-                    wallPositions.push_back(global);
-                    for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
-                    {
-                        if (abs(intersection.centerUnitOuterNormal()[dimIdx]) > 1e-8)
-                            wallNormalAxisTemp.push_back(dimIdx);
-                    }
+                    wallPositions.push_back(scvf.center());
+                    wallNormalAxisTemp.push_back(scvf.directionIndex());
                 }
             }
         }
@@ -417,12 +416,22 @@ public:
         }
     }
 
+    /*!
+     * \brief Returns whether a given sub control volume face is on a wall
+     *
+     * \param scvf The sub control volume face.
+     */
+    bool isOnWall(const SubControlVolumeFace& scvf) const
+    {
+        return asImp_().isOnWallAtPos(scvf.center());
+    }
+
     /*!
      * \brief Returns whether a given point is on a wall
      *
      * \param globalPos The position in global coordinates.
      */
-    bool isOnWall(const GlobalPosition &globalPos) const
+    bool isOnWallAtPos(const GlobalPosition &globalPos) const
     {
         // Throw an exception if no walls are implemented
         DUNE_THROW(Dune::InvalidStateException,
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
index 8333634b66..f22d5b088e 100644
--- a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
+++ b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
@@ -312,7 +312,7 @@ public:
     {
         unsigned int elementIdx = asImp_().fvGridGeometry().elementMapper().index(element);
         auto bcTypes = asImp_().boundaryTypes(element, localSubFace);
-        return asImp_().isOnWall(localSubFace.center())
+        return asImp_().isOnWall(localSubFace)
                && bcTypes.isDirichlet(eqIdx)
                && isMatchingPoint(elementIdx);
     }
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh
index 1ee5672027..621e11e3cb 100644
--- a/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh
+++ b/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh
@@ -154,7 +154,7 @@ public:
         const auto bcTypes = problem.boundaryTypes(element, scvf);
         if (!(scvf.boundary() && (bcTypes.isOutflow(Indices::turbulentKineticEnergyEqIdx)
                                   || bcTypes.isSymmetry()
-                                  || problem.isOnWall(scvf.center()))))
+                                  || problem.isOnWall(scvf))))
         {
             if (!(insideVolVars.isMatchingPoint() && outsideVolVars.isMatchingPoint())
                 || !(insideVolVars.isMatchingPoint() && outsideVolVars.inNearWallRegion())
diff --git a/test/freeflow/rans/pipelauferproblem.hh b/test/freeflow/rans/pipelauferproblem.hh
index 3441e79575..d7d1af40c3 100644
--- a/test/freeflow/rans/pipelauferproblem.hh
+++ b/test/freeflow/rans/pipelauferproblem.hh
@@ -185,7 +185,7 @@ public:
      */
     // \{
 
-    bool isOnWall(const GlobalPosition &globalPos) const
+    bool isOnWallAtPos(const GlobalPosition &globalPos) const
     {
         return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_
                || globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_;
@@ -260,7 +260,7 @@ public:
 #endif
 #if KOMEGA
             // set a fixed dissipation (omega) in one cell
-            if (isOnWall(globalPos))
+            if (isOnWallAtPos(globalPos))
                 values.setDirichletCell(Indices::dissipationIdx);
 #endif
 #if ONEEQ
@@ -285,7 +285,7 @@ public:
         if (time() > 10.0)
         {
             values[Indices::temperatureIdx] = inletTemperature_;
-            if (isOnWall(globalPos))
+            if (isOnWallAtPos(globalPos))
             {
                 values[Indices::temperatureIdx] = wallTemperature_;
             }
@@ -327,14 +327,14 @@ public:
         values[Indices::velocityXIdx] = time() > initializationTime_
                                         ? inletVelocity_
                                         : time() / initializationTime_ * inletVelocity_;
-        if (isOnWall(globalPos))
+        if (isOnWallAtPos(globalPos))
         {
             values[Indices::velocityXIdx] = 0.0;
         }
 
 #if NONISOTHERMAL
         values[Indices::temperatureIdx] = inletTemperature_;
-        if (isOnWall(globalPos))
+        if (isOnWallAtPos(globalPos))
         {
             values[Indices::temperatureIdx] = wallTemperature_;
         }
@@ -343,7 +343,7 @@ public:
 #if LOWREKEPSILON || KEPSILON || KOMEGA
         values[Indices::turbulentKineticEnergyIdx] = turbulentKineticEnergy_;
         values[Indices::dissipationIdx] = dissipation_;
-        if (isOnWall(globalPos))
+        if (isOnWallAtPos(globalPos))
         {
             values[Indices::turbulentKineticEnergyIdx] = 0.0;
             values[Indices::dissipationIdx] = 0.0;
@@ -352,7 +352,7 @@ public:
 
 #if ONEEQ
         values[Indices::viscosityTildeIdx] = viscosityTilde_;
-        if (isOnWall(globalPos))
+        if (isOnWallAtPos(globalPos))
         {
             values[Indices::viscosityTildeIdx] = 0.0;
         }
diff --git a/test/freeflow/ransnc/flatplatetestproblem.hh b/test/freeflow/ransnc/flatplatetestproblem.hh
index 39bffa6f7c..6b2a71660d 100644
--- a/test/freeflow/ransnc/flatplatetestproblem.hh
+++ b/test/freeflow/ransnc/flatplatetestproblem.hh
@@ -266,7 +266,7 @@ public:
             values.setOutflow(Indices::viscosityTildeIdx);
 #endif
         }
-        else if(isOnWall(globalPos))
+        else if(isOnWallAtPos(globalPos))
         {
             values.setDirichlet(Indices::velocityXIdx);
             values.setDirichlet(Indices::velocityYIdx);
@@ -316,7 +316,7 @@ public:
                 values[transportCompIdx] = 1e-3;
             }
 #if NONISOTHERMAL
-            if (isOnWall(globalPos))
+            if (isOnWallAtPos(globalPos))
             {
                 values[Indices::temperatureIdx] += 30.;
             }
@@ -370,14 +370,14 @@ public:
 
         // block velocity profile
         values[Indices::velocityXIdx] = 0.0;
-        if (!isOnWall(globalPos))
+        if (!isOnWallAtPos(globalPos))
             values[Indices::velocityXIdx] =  inletVelocity_;
         values[Indices::velocityYIdx] = 0.0;
 
 #if KEPSILON || KOMEGA || LOWREKEPSILON
         values[Indices::turbulentKineticEnergyEqIdx] = turbulentKineticEnergy_;
         values[Indices::dissipationEqIdx] = dissipation_;
-        if (isOnWall(globalPos))
+        if (isOnWallAtPos(globalPos))
         {
             values[Indices::turbulentKineticEnergyEqIdx] = 0.0;
             values[Indices::dissipationEqIdx] = 0.0;
@@ -386,7 +386,7 @@ public:
 
 #if ONEEQ
         values[Indices::viscosityTildeIdx] = viscosityTilde_;
-        if (isOnWall(globalPos))
+        if (isOnWallAtPos(globalPos))
         {
             values[Indices::viscosityTildeIdx] = 0.0;
         }
@@ -407,7 +407,7 @@ public:
         return timeLoop_->time();
     }
 
-    bool isOnWall(const GlobalPosition& globalPos) const
+    bool isOnWallAtPos(const GlobalPosition& globalPos) const
     {
         return globalPos[1] < eps_;
     }
-- 
GitLab