From 2d6c691e750efe8ffba5eea2cd655e795d683a68 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Thu, 8 Oct 2009 18:50:46 +0100
Subject: [PATCH] Adding instrumentation Pout statements.

Using an additional test for best guess face, which is the nearest
face, although this is not fully thought out, and can be a problem if
the nearest face is on a patch and this is not appropriate.
---
 src/lagrangian/basic/Cloud/Cloud.C       |  15 +--
 src/lagrangian/basic/Particle/Particle.C | 150 ++++++++++++++++++-----
 2 files changed, 127 insertions(+), 38 deletions(-)

diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
index 5d6cc5705a..81bf07b378 100644
--- a/src/lagrangian/basic/Cloud/Cloud.C
+++ b/src/lagrangian/basic/Cloud/Cloud.C
@@ -39,10 +39,10 @@ template<class ParticleType>
 const Foam::scalar Foam::Cloud<ParticleType>::minValidTrackFraction = 1e-10;
 
 template<class ParticleType>
-const Foam::scalar Foam::Cloud<ParticleType>::trackingRescueTolerance = 1e-3;
+const Foam::scalar Foam::Cloud<ParticleType>::trackingRescueTolerance = 1e-6;
 
 template<class ParticleType>
-const Foam::scalar Foam::Cloud<ParticleType>::intersectionTolerance = SMALL;
+const Foam::scalar Foam::Cloud<ParticleType>::intersectionTolerance = 1e-3;
 
 template<class ParticleType>
 const Foam::scalar Foam::Cloud<ParticleType>::planarCosAngle = (1 - 1e-6);
@@ -189,13 +189,14 @@ void Foam::Cloud<ParticleType>::calcConcaveCells() const
     intrudesIntoOwnerPtr_.reset(new PackedBoolList(pMesh().nFaces()));
     intrudesIntoNeighbourPtr_.reset(new PackedBoolList(pMesh().nFaces()));
 
-
     forAll(cells, cellI)
     {
-        if (isConcaveCell(cellI))
-        {
-            concaveCell[cellI] = 1;
-        }
+        // if (isConcaveCell(cellI))
+        // {
+        //     concaveCell[cellI] = 1;
+        // }
+
+        concaveCell[cellI] = 1;
     }
 
     {
diff --git a/src/lagrangian/basic/Particle/Particle.C b/src/lagrangian/basic/Particle/Particle.C
index b2dc2d3299..cbc4a9ccc0 100644
--- a/src/lagrangian/basic/Particle/Particle.C
+++ b/src/lagrangian/basic/Particle/Particle.C
@@ -118,6 +118,14 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
         return;
     }
 
+    Pout<< nl << origProc_ << " "
+        << origId_ << " "
+        << position_ << " "
+        << endPosition << " "
+        << celli_ << nl
+        << potentiallyCrossedFaces
+        << endl;
+
     label nFaceCrossings = 0;
 
     forAll (potentiallyCrossedFaces, pCFI)
@@ -156,11 +164,13 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
     if (nFaceCrossings > 1)
     {
         Pout<< "In cell " << celli_ << " there were " << nFaceCrossings
-            << " face crossings detected tracking from cell centre to "
+            << " face crossings detected tracking from concave cell centre to "
             << " endPosition"
             << endl;
     }
 
+    Pout<< nFaceCrossings << " " << (nFaceCrossings % 2) << endl;
+
     if (nFaceCrossings % 2 == 0)
     {
         // Even number of face crossings, so the particle must end up
@@ -201,7 +211,7 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
                 // has been found
 
                 Pout<< "New face crossing " << facei
-                    << " of cell " << celli_ << " found"
+                    << " of concave cell " << celli_ << " found"
                     << endl;
 
                 potentiallyCrossedFaces.append(facei);
@@ -209,6 +219,8 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
         }
     }
 
+    Pout<< potentiallyCrossedFaces << endl;
+
     vector deltaPosition = endPosition - position_;
 
     scalar tmpLambda = GREAT;
@@ -261,9 +273,19 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
     {
         // No face has been identified to be crossed yet, but the cell
         // must have been left, so the best guess of which face to
-        // cross is required.  This face will not be one of those
-        // faces already tested, and could be any other face of the
-        // cell.
+        // cross is required.
+
+        vector delta =
+            mag(mesh.cellCentres()[celli_] - position_)
+           *deltaPosition/mag(deltaPosition);
+
+        Pout<< "Need a best guess " << nl
+            << position_ << nl
+            << position_ + delta << nl
+            << position_ - delta
+            << endl;
+
+        // ---- Best guess by projection --
 
         DynamicList<scalar> tmpLambdas;
         DynamicList<label> tmpLambdaFaceIs;
@@ -272,34 +294,59 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
         {
             label facei = faces[i];
 
-            if (findIndex(potentiallyCrossedFaces, facei) == -1)
+            // Use exact FULL_RAY intersection to allow negative
+            // values
+
+            // TODO: A correction is required for moving meshes to
+            // calculate the correct lambda value.
+
+            pointHit inter = mesh.faces()[facei].intersection
+            (
+                position_,
+                deltaPosition,
+                mesh.faceCentres()[facei],
+                mesh.points(),
+                intersection::FULL_RAY,
+                Cloud<ParticleType>::intersectionTolerance
+            );
+
+            if (inter.hit())
             {
-                // This face has not been assessed yet.
-
-                // Use exact FULL_RAY intersection to allow negative
-                // values
-
-                // TODO: A correction is required for moving meshes to
-                // calculate the correct lambda value.
-
-                pointHit inter = mesh.faces()[facei].intersection
-                (
-                    position_,
-                    deltaPosition,
-                    mesh.faceCentres()[facei],
-                    mesh.points(),
-                    intersection::FULL_RAY,
-                    Cloud<ParticleType>::intersectionTolerance
-                );
-
-                if (inter.hit())
-                {
-                    tmpLambdas.append(inter.distance());
-                    tmpLambdaFaceIs.append(facei);
-                }
+                tmpLambdas.append(inter.distance());
+                tmpLambdaFaceIs.append(facei);
             }
+
+            inter = mesh.faces()[facei].intersection
+            (
+                position_,
+                deltaPosition,
+                mesh.faceCentres()[facei],
+                mesh.points(),
+                intersection::HALF_RAY,
+                Cloud<ParticleType>::intersectionTolerance
+            );
+
+            Pout<< "Test repeat HALF "
+                << facei << " " << inter.distance() << endl;
+
+            inter = mesh.faces()[facei].intersection
+            (
+                position_,
+                deltaPosition,
+                mesh.faceCentres()[facei],
+                mesh.points(),
+                intersection::FULL_RAY,
+                Cloud<ParticleType>::intersectionTolerance
+            );
+
+            Pout<< "Test repeat FULL "
+                << facei << " " << inter.distance() << nl
+                << (inter.distance()*deltaPosition) + position_
+                << endl;
         }
 
+        Pout<< tmpLambdaFaceIs << " " << tmpLambdas << endl;
+
         // The closest lambda value (less than zero) or (greater than or
         // equal to one), and corresponding face.
         scalar closestLambda = GREAT;
@@ -343,12 +390,48 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
         }
         else
         {
-            Pout<< "Tracking failure in concave cell - didn't find a "
-                << "best guess face" << endl;
+            Pout<< "Tracking in concave cell - didn't find a "
+                << "best guess face" << nl
+                << origId_ << " "
+                << origProc_ << " "
+                << position_ << " "
+                << endPosition << " "
+                << stepFraction_ << " "
+                << celli_ << " "
+                << endl;
+
+            // ---- Best guess by nearest ----
+            scalar minDistance = GREAT;
+
+            forAll(faces, i)
+            {
+                label facei = faces[i];
+
+                if (cloud_.internalFace(facei))
+                {
+                    pointHit nearest = mesh.faces()[facei].nearestPoint
+                    (
+                        position_,
+                        mesh.points()
+                    );
+
+                    if (nearest.distance() < minDistance)
+                    {
+                        minDistance = nearest.distance();
+                        facei_ = facei;
+                    }
+                }
+            }
+
+            // ---- End of best guess by nearest ----
         }
     }
 
+    Pout<< facei_ << " " << celli_ << endl;
+
     faceAction(trackFraction, endPosition, td);
+
+    Pout<< facei_ << " " << celli_ << endl;
 }
 
 
@@ -722,6 +805,11 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
             // Use a more careful tracking algorithm if the cell is concave
             trackToFaceConcave(trackFraction, endPosition, td);
         }
+        else
+        {
+            // Use a more careful tracking algorithm if the cell is concave
+            trackToFaceConvex(trackFraction, endPosition, td);
+        }
     }
     else
     {