From adfc2de66e35bae5a7788cbed57b01ad72baafb9 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Thu, 7 Apr 2011 11:23:18 +0100
Subject: [PATCH] ENH: Moving the calculation of the lambda
 denominator/numerator upwards.

Expensive to look it up in the lambda function, passing from trackToFace as an
argument.
---
 src/lagrangian/basic/particle/particle.C      |  2 ++
 src/lagrangian/basic/particle/particle.H      | 13 ++++++---
 src/lagrangian/basic/particle/particleI.H     | 19 ++++++-------
 .../basic/particle/particleTemplates.C        | 27 ++++++++++++++++---
 4 files changed, 45 insertions(+), 16 deletions(-)

diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C
index 38c45bf57f..f3566c93e6 100644
--- a/src/lagrangian/basic/particle/particle.C
+++ b/src/lagrangian/basic/particle/particle.C
@@ -32,6 +32,8 @@ Foam::label Foam::particle::particleCount_ = 0;
 
 const Foam::scalar Foam::particle::trackingCorrectionTol = 1e-5;
 
+const Foam::scalar Foam::particle::lambdaDistanceToleranceCoeff = 1e3*SMALL;
+
 namespace Foam
 {
     defineTypeNameAndDebug(particle, 0);
diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H
index df64e7b9e7..06d32c321c 100644
--- a/src/lagrangian/basic/particle/particle.H
+++ b/src/lagrangian/basic/particle/particle.H
@@ -164,7 +164,8 @@ protected:
             DynamicList<label>& faceList,
             const tetPointRef& tet,
             const FixedList<vector, 4>& tetAreas,
-            const FixedList<label, 4>& tetPlaneBasePtIs
+            const FixedList<label, 4>& tetPlaneBasePtIs,
+            const scalar tol
         ) const;
 
         //- Find the lambda value for the line to-from across the
@@ -178,7 +179,8 @@ protected:
             const label tetPlaneBasePtI,
             const label cellI,
             const label tetFaceI,
-            const label tetPtI
+            const label tetPtI,
+            const scalar tol
         ) const;
 
         //- Find the lambda value for a moving tri face
@@ -191,7 +193,8 @@ protected:
             const label tetPlaneBasePtI,
             const label cellI,
             const label tetFaceI,
-            const label tetPtI
+            const label tetPtI,
+            const scalar tol
         ) const;
 
         //- Modify the tet owner data by crossing triI
@@ -291,6 +294,10 @@ public:
         // 'rescue' it from a tracking problem
         static const scalar trackingCorrectionTol;
 
+        //- Fraction of the cell volume to use in determining tolerance values
+        //  for the denominator and numerator of lambda
+        static const scalar lambdaDistanceToleranceCoeff;
+
 
     // Constructors
 
diff --git a/src/lagrangian/basic/particle/particleI.H b/src/lagrangian/basic/particle/particleI.H
index a11018a5d5..6107df0c4e 100644
--- a/src/lagrangian/basic/particle/particleI.H
+++ b/src/lagrangian/basic/particle/particleI.H
@@ -34,7 +34,8 @@ inline void Foam::particle::findTris
     DynamicList<label>& faceList,
     const tetPointRef& tet,
     const FixedList<vector, 4>& tetAreas,
-    const FixedList<label, 4>& tetPlaneBasePtIs
+    const FixedList<label, 4>& tetPlaneBasePtIs,
+    const scalar tol
 ) const
 {
     faceList.clear();
@@ -52,7 +53,8 @@ inline void Foam::particle::findTris
             tetPlaneBasePtIs[i],
             cellI_,
             tetFaceI_,
-            tetPtI_
+            tetPtI_,
+            tol
         );
 
         if ((lambda > 0.0) && (lambda < 1.0))
@@ -72,7 +74,8 @@ inline Foam::scalar Foam::particle::tetLambda
     const label tetPlaneBasePtI,
     const label cellI,
     const label tetFaceI,
-    const label tetPtI
+    const label tetPtI,
+    const scalar tol
 ) const
 {
     const pointField& pPts = mesh_.points();
@@ -88,7 +91,8 @@ inline Foam::scalar Foam::particle::tetLambda
             tetPlaneBasePtI,
             cellI,
             tetFaceI,
-            tetPtI
+            tetPtI,
+            tol
         );
     }
 
@@ -102,8 +106,6 @@ inline Foam::scalar Foam::particle::tetLambda
     // delta-length in the direction of n times the face area to a fraction of
     // the cell volume.
 
-    scalar tol = 1e3*SMALL*mesh_.cellVolumes()[cellI];
-
     if (mag(lambdaDenominator) < tol)
     {
         if (mag(lambdaNumerator) < tol)
@@ -147,7 +149,8 @@ inline Foam::scalar Foam::particle::movingTetLambda
     const label tetPlaneBasePtI,
     const label cellI,
     const label tetFaceI,
-    const label tetPtI
+    const label tetPtI,
+    const scalar tol
 ) const
 {
     const pointField& pPts = mesh_.points();
@@ -299,8 +302,6 @@ inline Foam::scalar Foam::particle::movingTetLambda
 
     }
 
-    scalar tol = 1e3*SMALL*mesh_.cellVolumes()[cellI];
-
     if (mag(lambdaDenominator) < tol)
     {
         if (mag(lambdaNumerator) < tol)
diff --git a/src/lagrangian/basic/particle/particleTemplates.C b/src/lagrangian/basic/particle/particleTemplates.C
index e77376150a..cb58c21653 100644
--- a/src/lagrangian/basic/particle/particleTemplates.C
+++ b/src/lagrangian/basic/particle/particleTemplates.C
@@ -286,6 +286,11 @@ Foam::scalar Foam::particle::trackToFace
     // be a different tet to the one that the particle occupies.
     tetIndices faceHitTetIs;
 
+    // What tolerance is appropriate the minimum lambda numerator and
+    // denominator for tracking in this cell.
+    scalar lambdaDistanceTolerance =
+        lambdaDistanceToleranceCoeff*mesh_.cellVolumes()[cellI_];
+
     do
     {
         if (triI != -1)
@@ -371,7 +376,15 @@ Foam::scalar Foam::particle::trackToFace
         tetPlaneBasePtIs[2] = basePtI;
         tetPlaneBasePtIs[3] = basePtI;
 
-        findTris(endPosition, tris, tet, tetAreas, tetPlaneBasePtIs);
+        findTris
+        (
+            endPosition,
+            tris,
+            tet,
+            tetAreas,
+            tetPlaneBasePtIs,
+            lambdaDistanceTolerance
+        );
 
         // Reset variables for new track
         triI = -1;
@@ -415,7 +428,8 @@ Foam::scalar Foam::particle::trackToFace
                     tetPlaneBasePtIs[tI],
                     cellI_,
                     tetFaceI_,
-                    tetPtI_
+                    tetPtI_,
+                    lambdaDistanceTolerance
                 );
 
                 if (lam < lambdaMin)
@@ -704,6 +718,9 @@ void Foam::particle::hitWallFaces
 
     const Foam::cell& thisCell = mesh_.cells()[cellI_];
 
+    scalar lambdaDistanceTolerance =
+        lambdaDistanceToleranceCoeff*mesh_.cellVolumes()[cellI_];
+
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
     forAll(thisCell, cFI)
@@ -755,7 +772,8 @@ void Foam::particle::hitWallFaces
                     f[tetIs.faceBasePt()],
                     cellI_,
                     fI,
-                    tetIs.tetPt()
+                    tetIs.tetPt(),
+                    lambdaDistanceTolerance
                 );
 
                 if ((tetClambda <= 0.0) || (tetClambda >= 1.0))
@@ -781,7 +799,8 @@ void Foam::particle::hitWallFaces
                     f[tetIs.faceBasePt()],
                     cellI_,
                     fI,
-                    tetIs.tetPt()
+                    tetIs.tetPt(),
+                    lambdaDistanceTolerance
                 );
 
                 pointHit hitInfo(vector::zero);