From 3e12a89a041f728471846cbc7828bbd4d837ea8b Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Mon, 2 Sep 2013 22:51:15 +0100
Subject: [PATCH] Check mesh is 2D and scale resulting steam-function by the
 thickness

---
 .../streamFunction/streamFunction.C           | 23 +++++++++++++++----
 1 file changed, 18 insertions(+), 5 deletions(-)

diff --git a/applications/utilities/postProcessing/velocityField/streamFunction/streamFunction.C b/applications/utilities/postProcessing/velocityField/streamFunction/streamFunction.C
index aa56086c90..de9af38b87 100644
--- a/applications/utilities/postProcessing/velocityField/streamFunction/streamFunction.C
+++ b/applications/utilities/postProcessing/velocityField/streamFunction/streamFunction.C
@@ -42,14 +42,26 @@ Description
 int main(int argc, char *argv[])
 {
     timeSelector::addOptions();
-#   include "addRegionOption.H"
+    #include "addRegionOption.H"
 
-#   include "setRootCase.H"
-#   include "createTime.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
 
     instantList timeDirs = timeSelector::select0(runTime, args);
 
-#   include "createNamedMesh.H"
+    #include "createNamedMesh.H"
+
+    label nD = mesh.nGeometricD();
+
+    if (nD != 2)
+    {
+        FatalErrorIn(args.executable())
+            << "Case is not 2D, stream-function cannot be computed"
+            << exit(FatalError);
+    }
+
+    const vector slabDir((Vector<label>::one - mesh.geometricD())/2);
+    scalar thickness = slabDir & mesh.bounds().span();
 
     const pointMesh& pMesh = pointMesh::New(mesh);
 
@@ -428,7 +440,6 @@ int main(int argc, char *argv[])
                      }
 
                      Info<< ".";
-//                     Info<< "One pass, n visited = " << nVisited << endl;
 
                      if (nVisited == nVisitedOld)
                      {
@@ -449,6 +460,8 @@ int main(int argc, char *argv[])
                 Info<< endl;
             } while (!finished);
 
+            // Normalise the stream-function by the 2D mesh thickness
+            streamFunction /= thickness;
             streamFunction.boundaryField() = 0.0;
             streamFunction.write();
         }