From 9ef96ef0a8e6c0959fe8db6199f60432cfdd6c4e Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Tue, 6 Oct 2009 10:21:43 +0100
Subject: [PATCH 01/16] Modified version of Mattijs' concave tracking
 modifications.  Changes to Mattijs' work are:

+ Correct use of normals and cosine (dot product) to identify planar
  faces.

+ Correct use for lambda in tracking - it is a fraction, not a distance.

+ Not doing a reduce on demand driven construction data - not all
  processors call it at the same time, so crashes.

This implementation contains an attempt at making the calculation of
lambdaC (from the cell centre) use decomposed triangles for faces, but
this is a bad approach, the concept of using lambdaC relies on the
definition of a convex polyhedron, i.e. none of the planes of the
faces of the polyhedron are inside the volume.  Can't use this method
and will need to treat a convex cell completely differently, not just
some of its faces.
---
 src/lagrangian/basic/Cloud/Cloud.C       | 285 ++++++++++++++++++++++-
 src/lagrangian/basic/Cloud/Cloud.H       |  65 +++++-
 src/lagrangian/basic/Cloud/CloudIO.C     |  13 +-
 src/lagrangian/basic/Particle/Particle.C | 193 +++++++++++++--
 4 files changed, 533 insertions(+), 23 deletions(-)

diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
index e3e0b4f63f..7c2b383b83 100644
--- a/src/lagrangian/basic/Cloud/Cloud.C
+++ b/src/lagrangian/basic/Cloud/Cloud.C
@@ -31,6 +31,250 @@ License
 #include "mapPolyMesh.H"
 #include "Time.H"
 #include "OFstream.H"
+#include "triPointRef.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<class ParticleType>
+const Foam::scalar Foam::Cloud<ParticleType>::minValidTrackFraction = 1e-10;
+
+template<class ParticleType>
+const Foam::scalar Foam::Cloud<ParticleType>::trackingRescueTolerance = 1e-3;
+
+template<class ParticleType>
+const Foam::scalar Foam::Cloud<ParticleType>::intersectionTolerance = 1e-6;
+
+template<class ParticleType>
+const Foam::scalar Foam::Cloud<ParticleType>::planarCosAngle = (1 - 1e-6);
+
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+// Calculate using face properties only
+template<class ParticleType>
+Foam::label Foam::Cloud<ParticleType>::isIntrudingFace
+(
+    const label cellI,
+    const label f0I,
+    const label f1I
+) const
+{
+    const vectorField& fAreas = pMesh().faceAreas();
+    const pointField& fCentres = pMesh().faceCentres();
+    const labelList& fOwner = pMesh().faceOwner();
+    const pointField& points = pMesh().points();
+
+    // Get f0 centre
+    const point& fc0 = fCentres[f0I];
+
+    const face& f0 = pMesh().faces()[f0I];
+
+    const face& f1 = pMesh().faces()[f1I];
+
+    vector fn0 = fAreas[f0I];
+    fn0 /= (mag(fn0) + VSMALL);
+
+    // Flip normal if required so that it is always pointing out of
+    // the cell
+    if (fOwner[f0I] != cellI)
+    {
+        fn0 *= -1;
+    }
+
+    // Is any vertex of f1 on wrong side of the plane of f0?
+
+    forAll(f1, f1pI)
+    {
+        label ptI = f1[f1pI];
+
+        // Skip points that are shared between f1 and f0
+        if (findIndex(f0, ptI) > -1)
+        {
+            continue;
+        }
+
+        const point& pt = points[ptI];
+
+        // Normal distance from fc0 to pt
+        scalar d = ((fc0 - pt) & fn0);
+
+        if (d < 0)
+        {
+            // Proper concave: f1 on wrong side of f0
+            return 2;
+        }
+    }
+
+    // Could be a convex cell, but check angle for planar faces.
+
+    vector fn1 = fAreas[f1I];
+    fn1 /= (mag(fn1) + VSMALL);
+
+    // Flip normal if required so that it is always pointing out of
+    // the cell
+    if (fOwner[f1I] != cellI)
+    {
+        fn1 *= -1;
+    }
+
+    if ((fn0 & fn1) > planarCosAngle)
+    {
+        // Planar face
+        return 1;
+    }
+
+    return 0;
+}
+
+
+template<class ParticleType>
+void Foam::Cloud<ParticleType>::calcIntrudingFaces() const
+{
+    const labelList& fOwner = pMesh().faceOwner();
+    const cellList& cells = pMesh().cells();
+
+    intrudesIntoOwnerPtr_.reset(new PackedBoolList(pMesh().nFaces()));
+    PackedBoolList& intrudesIntoOwner = intrudesIntoOwnerPtr_();
+    intrudesIntoNeighbourPtr_.reset(new PackedBoolList(pMesh().nFaces()));
+    PackedBoolList& intrudesIntoNeighbour = intrudesIntoNeighbourPtr_();
+
+    PackedBoolList cellsWithProblemFaces(cells.size());
+    DynamicList<label> planarFaces;
+    DynamicList<label> intrudingFaces;
+
+
+    label nIntruded = 0;
+    label nPlanar = 0;
+
+    forAll(cells, cellI)
+    {
+        const cell& cFaces = cells[cellI];
+
+        forAll(cFaces, i)
+        {
+            label f0 = cFaces[i];
+
+            bool own0 = (fOwner[f0] == cellI);
+
+            // Now check that all other faces of the cell are on the 'inside'
+            // of the face.
+
+            bool intrudes = false;
+
+            forAll(cFaces, j)
+            {
+                if (j != i)
+                {
+                    label state = isIntrudingFace
+                    (
+                        cellI,
+                        f0,
+                        cFaces[j]
+                    );
+
+                    if (state == 1)
+                    {
+                        // planar
+                        planarFaces.append(f0);
+
+                        intrudes = true;
+                        nPlanar++;
+                        break;
+                    }
+                    else if (state == 2)
+                    {
+                        // concave
+                        intrudingFaces.append(f0);
+
+                        intrudes = true;
+                        nIntruded++;
+                        break;
+                    }
+                }
+            }
+
+            if (intrudes)
+            {
+                cellsWithProblemFaces[cellI] = 1;
+
+                if (own0)
+                {
+                    intrudesIntoOwner[f0] = 1;
+                }
+                else
+                {
+                    intrudesIntoNeighbour[f0] = 1;
+                }
+            }
+
+            // intrudesIntoOwner[f0] = 1;
+            // intrudesIntoNeighbour[f0] = 1;
+        }
+    }
+
+    // if (debug)
+    {
+        Pout<< "Cloud<ParticleType>::calcIntrudingFaces() :"
+            << " overall faces in mesh : "
+            << pMesh().nFaces() << nl
+            << "    of these planar : "
+            << nPlanar << nl
+            << "    of these intruding into their cell (concave) : "
+            << nIntruded << endl;
+
+        // Info<< "Hard coded picking up every face." << endl;
+
+        // Write the faces and cells that are a problem to file to be
+        // compared and made into sets
+
+        {
+            DynamicList<label> cellsWithIntrudingFaces;
+
+            forAll(cells, cellI)
+            {
+                if (cellsWithProblemFaces[cellI])
+                {
+                    cellsWithIntrudingFaces.append(cellI);
+                }
+            }
+
+            fileName fName = pMesh().time().path()/"cellsWithProblemFaces";
+
+            Pout<< "    Writing " << fName.name() << endl;
+
+            OFstream file(fName);
+
+            file << cellsWithIntrudingFaces;
+
+            file.flush();
+        }
+
+        {
+            fileName fName = pMesh().time().path()/"planarFaces";
+
+            Pout<< "    Writing " << fName.name() << endl;
+
+            OFstream file(fName);
+
+            file << planarFaces;
+
+            file.flush();
+        }
+
+        {
+            fileName fName = pMesh().time().path()/"intrudingFaces";
+
+            Pout<< "    Writing " << fName.name() << endl;
+
+            OFstream file(fName);
+
+            file << intrudingFaces;
+
+            file.flush();
+        }
+    }
+}
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -38,12 +282,14 @@ template<class ParticleType>
 Foam::Cloud<ParticleType>::Cloud
 (
     const polyMesh& pMesh,
-    const IDLList<ParticleType>& particles
+    const IDLList<ParticleType>& particles,
+    const bool concaveCheck
 )
 :
     cloud(pMesh),
     IDLList<ParticleType>(),
     polyMesh_(pMesh),
+    concaveCheck_(concaveCheck),
     particleCount_(0)
 {
     IDLList<ParticleType>::operator=(particles);
@@ -55,20 +301,55 @@ Foam::Cloud<ParticleType>::Cloud
 (
     const polyMesh& pMesh,
     const word& cloudName,
-    const IDLList<ParticleType>& particles
+    const IDLList<ParticleType>& particles,
+    const bool concaveCheck
 )
 :
     cloud(pMesh, cloudName),
     IDLList<ParticleType>(),
     polyMesh_(pMesh),
+    concaveCheck_(concaveCheck),
     particleCount_(0)
 {
     IDLList<ParticleType>::operator=(particles);
 }
 
 
+template<class ParticleType>
+void Foam::Cloud<ParticleType>::clearOut()
+{
+    intrudesIntoOwnerPtr_.clear();
+    intrudesIntoNeighbourPtr_.clear();
+    labels_.clearStorage();
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class ParticleType>
+const Foam::PackedBoolList& Foam::Cloud<ParticleType>::intrudesIntoOwner()
+const
+{
+    if (!intrudesIntoOwnerPtr_.valid())
+    {
+        calcIntrudingFaces();
+    }
+    return intrudesIntoOwnerPtr_();
+}
+
+
+template<class ParticleType>
+const Foam::PackedBoolList& Foam::Cloud<ParticleType>::intrudesIntoNeighbour()
+const
+{
+    if (!intrudesIntoNeighbourPtr_.valid())
+    {
+        calcIntrudingFaces();
+    }
+    return intrudesIntoNeighbourPtr_();
+}
+
+
 template<class ParticleType>
 Foam::label Foam::Cloud<ParticleType>::getNewParticleID() const
 {
diff --git a/src/lagrangian/basic/Cloud/Cloud.H b/src/lagrangian/basic/Cloud/Cloud.H
index 626042e68a..11a0809c8f 100644
--- a/src/lagrangian/basic/Cloud/Cloud.H
+++ b/src/lagrangian/basic/Cloud/Cloud.H
@@ -40,6 +40,7 @@ SourceFiles
 #include "IDLList.H"
 #include "IOField.H"
 #include "polyMesh.H"
+#include "PackedBoolList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -75,6 +76,13 @@ class Cloud
 
         const polyMesh& polyMesh_;
 
+        const bool concaveCheck_;
+
+        //- Does face intrude into owner cell
+        mutable autoPtr<PackedBoolList> intrudesIntoOwnerPtr_;
+        //- Does face intrude into neighbour cell
+        mutable autoPtr<PackedBoolList> intrudesIntoNeighbourPtr_;
+
         //- Overall count of particles ever created. Never decreases.
         mutable label particleCount_;
 
@@ -84,6 +92,21 @@ class Cloud
 
     // Private member functions
 
+        //- Helper: is face f1 on wrong side of f0 (both faces are on same cell)
+        //  0 : convex
+        //  1 : planar to within planarCosAngle
+        //  2 : concave
+        label isIntrudingFace
+        (
+            const label cellI,
+            const label f0,
+            const label f1
+        ) const;
+
+        //- Analyse mesh for faces not being convex.
+        //  Sets intrudesIntoOwner_,intrudesIntoNeighbour_.
+        void calcIntrudingFaces() const;
+
         //- Initialise cloud on IO constructor
         void initCloud(const bool checkClass);
 
@@ -104,13 +127,31 @@ public:
     TypeName("Cloud");
 
 
+    // Static data members
+
+        //- The smallest allowable trackFraction in trackToFace before
+        //  an intervention is required
+        static const scalar minValidTrackFraction;
+
+        //- Fraction of distance to cell centre to move a particle to
+        //  'rescue' it from a tracking problem
+        static const scalar trackingRescueTolerance;
+
+        //- Overlap tolerance for face intersections.
+        static const scalar intersectionTolerance;
+
+        //- cosine of angle for faces to be planar.
+        static const scalar planarCosAngle;
+
+
     // Constructors
 
         //- Construct from mesh and a list of particles
         Cloud
         (
             const polyMesh& mesh,
-            const IDLList<ParticleType>& particles
+            const IDLList<ParticleType>& particles,
+            const bool concaveCheck = true
         );
 
         //- Construct from mesh, cloud name, and a list of particles
@@ -118,7 +159,8 @@ public:
         (
             const polyMesh& mesh,
             const word& cloudName,
-            const IDLList<ParticleType>& particles
+            const IDLList<ParticleType>& particles,
+            const bool concaveCheck = true
         );
 
         //- Construct from mesh by reading from file
@@ -126,7 +168,8 @@ public:
         Cloud
         (
             const polyMesh& mesh,
-            const bool checkClass = true
+            const bool checkClass = true,
+            const bool concaveCheck = true
         );
 
 
@@ -136,10 +179,17 @@ public:
         (
             const polyMesh& pMesh,
             const word& cloudName,
-            const bool checkClass = true
+            const bool checkClass = true,
+            const bool concaveCheck = true
         );
 
 
+    // Destructor
+
+        //- Destroy any demand-driven data
+        void clearOut();
+
+
     // Member Functions
 
         // Access
@@ -179,6 +229,13 @@ public:
                 return IDLList<ParticleType>::size();
             };
 
+            //- Per face whether it intrudes into the owner cell.(demand driven)
+            const PackedBoolList& intrudesIntoOwner() const;
+
+            //- Per face whether it intrudes into the neighbour cell. (demand
+            //  driven.
+            const PackedBoolList& intrudesIntoNeighbour() const;
+
 
             // Iterators
 
diff --git a/src/lagrangian/basic/Cloud/CloudIO.C b/src/lagrangian/basic/Cloud/CloudIO.C
index ac23db92e4..1eb35e09d0 100644
--- a/src/lagrangian/basic/Cloud/CloudIO.C
+++ b/src/lagrangian/basic/Cloud/CloudIO.C
@@ -62,11 +62,13 @@ template<class ParticleType>
 Foam::Cloud<ParticleType>::Cloud
 (
     const polyMesh& pMesh,
-    const bool checkClass
+    const bool checkClass,
+    const bool concaveCheck
 )
 :
     cloud(pMesh),
     polyMesh_(pMesh),
+    concaveCheck_(concaveCheck),
     particleCount_(0)
 {
     initCloud(checkClass);
@@ -78,14 +80,21 @@ Foam::Cloud<ParticleType>::Cloud
 (
     const polyMesh& pMesh,
     const word& cloudName,
-    const bool checkClass
+    const bool checkClass,
+    const bool concaveCheck
 )
 :
     cloud(pMesh, cloudName),
     polyMesh_(pMesh),
+    concaveCheck_(concaveCheck),
     particleCount_(0)
 {
     initCloud(checkClass);
+
+    // if (concaveCheck_)
+    // {
+    //     intrudesIntoOwner();
+    // }
 }
 
 
diff --git a/src/lagrangian/basic/Particle/Particle.C b/src/lagrangian/basic/Particle/Particle.C
index d99dfb476c..dc0c383ea4 100644
--- a/src/lagrangian/basic/Particle/Particle.C
+++ b/src/lagrangian/basic/Particle/Particle.C
@@ -47,10 +47,68 @@ void Foam::Particle<ParticleType>::findFaces
     const vector& C = mesh.cellCentres()[celli_];
 
     faceList.clear();
+
+    // forAll(faces, i)
+    // {
+    //     label facei = faces[i];
+    //     scalar lam = lambda(C, position, facei);
+
+    //     if ((lam > 0) && (lam < 1.0))
+    //     {
+    //         faceList.append(facei);
+    //     }
+    // }
+
+    vector deltaPosition = position - C;
+    scalar deltaPositionMag = mag(deltaPosition);
+    scalar lam;
+
     forAll(faces, i)
     {
         label facei = faces[i];
-        scalar lam = lambda(C, position, facei);
+
+        lam = GREAT;
+
+        if (cloud_.concaveCheck_)
+        {
+            // Check for intruding face
+            // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+            label isOwner = (celli_ == mesh.faceOwner()[facei]);
+
+            if
+            (
+                (isOwner && cloud_.intrudesIntoOwner()[facei])
+             || (!isOwner && cloud_.intrudesIntoNeighbour()[facei])
+            )
+            {
+                // Concave face. Use exact intersection.
+
+                pointHit inter = mesh.faces()[facei].intersection
+                (
+                    position_,
+                    deltaPosition,
+                    mesh.faceCentres()[facei],
+                    mesh.points(),
+                    intersection::HALF_RAY,
+                    Cloud<ParticleType>::intersectionTolerance
+                );
+
+                if (inter.hit())
+                {
+                    lam = inter.distance()/deltaPositionMag;
+                }
+            }
+            else
+            {
+                // Convex face. Can use lambda.
+                lam = lambda(C, position, facei);
+            }
+        }
+        else
+        {
+            lam = lambda(C, position, facei);
+        }
 
         if ((lam > 0) && (lam < 1.0))
         {
@@ -241,6 +299,17 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     facei_ = -1;
     scalar trackFraction = 0.0;
 
+    scalar lambdaMin = GREAT;
+
+    // Temporary addition to report on infinite loops ->
+    List<scalar> lambdas(faces.size());
+
+    forAll(faces, i)
+    {
+        lambdas[i] = lambda(position_, endPosition, faces[i], stepFraction_);
+    }
+    // <- Temporary addition to report on infinite loops
+
     if (faces.empty())  // inside cell
     {
         trackFraction = 1.0;
@@ -248,7 +317,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     }
     else // hit face
     {
-        scalar lambdaMin = GREAT;
+        // scalar lambdaMin = GREAT;
 
         if (faces.size() == 1)
         {
@@ -257,21 +326,86 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
         }
         else
         {
-            // If the particle has to cross more than one cell to reach the
-            // endPosition, we check which way to go.
-            // If one of the faces is a boundary face and the particle is
-            // outside, we choose the boundary face.
+            // // If the particle has to cross more than one cell to reach the
+            // // endPosition, we check which way to go.
+            // // If one of the faces is a boundary face and the particle is
+            // // outside, we choose the boundary face.
             // The particle is outside if one of the lambda's is > 1 or < 0
+
+            scalar fallBackLambdaMin = GREAT;
+            label fallBackFacei = -1;
+
+            vector deltaPosition = endPosition - position_;
+            scalar deltaPositionMag = mag(deltaPosition);
+            scalar hitLambda = GREAT;
+
             forAll(faces, i)
             {
                 scalar lam =
                     lambda(position_, endPosition, faces[i], stepFraction_);
 
-                if (lam < lambdaMin)
+                // Calculate fallback
+                if (lam < fallBackLambdaMin)
                 {
-                    lambdaMin = lam;
-                    facei_ = faces[i];
+                    fallBackLambdaMin = lam;
+                    fallBackFacei = faces[i];
                 }
+
+                if (cloud_.concaveCheck_)
+                {
+                    // Check for intruding face
+                    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+                    label isOwner = (celli_ == mesh.faceOwner()[faces[i]]);
+
+                    if
+                    (
+                        (isOwner && cloud_.intrudesIntoOwner()[faces[i]])
+                     || (!isOwner && cloud_.intrudesIntoNeighbour()[faces[i]])
+                    )
+                    {
+                        // Concave face. Use exact intersection.
+
+                        pointHit inter = mesh.faces()[faces[i]].intersection
+                        (
+                            position_,
+                            deltaPosition,
+                            mesh.faceCentres()[faces[i]],
+                            mesh.points(),
+                            intersection::HALF_RAY,
+                            Cloud<ParticleType>::intersectionTolerance
+                        );
+
+                        if (inter.hit())
+                        {
+                            hitLambda = inter.distance()/deltaPositionMag;
+
+                            if (hitLambda < lambdaMin)
+                            {
+                                lambdaMin = hitLambda;
+                                facei_ = faces[i];
+                            }
+                        }
+                    }
+                    else
+                    {
+                        // Convex face. Can use lambda.
+
+                        if (lam < lambdaMin)
+                        {
+                            lambdaMin = lam;
+                            facei_ = faces[i];
+                        }
+                    }
+                }
+            }
+
+            // Fall back to the original algorithm if no
+            // intersection was found
+            if (facei_ == -1)
+            {
+                lambdaMin = fallBackLambdaMin;
+                facei_ = fallBackFacei;
             }
         }
 
@@ -320,9 +454,8 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
                 FatalErrorIn
                 (
                     "Particle::trackToFace(const vector&, TrackData&)"
-                )
-                    << "addressing failure" << nl
-                    << abort(FatalError);
+                )<< "addressing failure" << nl
+                 << abort(FatalError);
             }
         }
         else
@@ -381,7 +514,6 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
                     p.hitPatch(patch, td);
                 }
             }
-
         }
     }
 
@@ -391,14 +523,44 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     // caused by face warpage.
     // In both cases resolve the positional ambiguity by moving the particle
     // slightly towards the cell-centre.
-    if (trackFraction < SMALL)
+
+    if (trackFraction < Cloud<ParticleType>::minValidTrackFraction)
     {
-        position_ += 1.0e-3*(mesh.cellCentres()[celli_] - position_);
+        Pout<< "tracking error "
+            << origId_ << " "
+            << origProc_ << " "
+            << position_ << " "
+            << endPosition << " "
+            << position_ + 1000*(endPosition - position_) << " "
+            << position_ - 1000*(endPosition - position_) << " "
+            << lambdaMin << " "
+            << trackFraction << " "
+            << stepFraction_ << " "
+            << celli_ << " "
+            << facei_ << " "
+            << cloud_.intrudesIntoOwner()[facei_] << " "
+            << cloud_.intrudesIntoNeighbour()[facei_] << " "
+            << faces << " "
+            << lambdas
+            << endl;
+
+        // if (stepFraction_  < SMALL)
+        // {
+        //     Pout<< "Intervene to mark face" << endl;
+
+        //     cloud_.intrudesIntoOwnerPtr_()[facei_] = 1;
+        //     cloud_.intrudesIntoNeighbourPtr_()[facei_] = 1;
+        // }
+
+        const point& cc = mesh.cellCentres()[celli_];
+        position_ +=
+            Cloud<ParticleType>::trackingRescueTolerance*(cc - position_);
     }
 
     return trackFraction;
 }
 
+
 template<class ParticleType>
 Foam::scalar Foam::Particle<ParticleType>::trackToFace
 (
@@ -409,6 +571,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     return trackToFace(endPosition, dummyTd);
 }
 
+
 template<class ParticleType>
 void Foam::Particle<ParticleType>::transformPosition(const tensor& T)
 {

From 8fd827aeb3ed7fa5510383d5fb0abbb9ba3b423f Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Tue, 6 Oct 2009 16:59:46 +0100
Subject: [PATCH 02/16] Moving to a cell approach for cells that are concave. 
 WIP - still problems.

---
 src/lagrangian/basic/Cloud/Cloud.C       | 297 ++++++----------
 src/lagrangian/basic/Cloud/Cloud.H       |  29 +-
 src/lagrangian/basic/Cloud/CloudIO.C     |   5 -
 src/lagrangian/basic/Particle/Particle.C | 409 ++++++++++-------------
 src/lagrangian/basic/Particle/Particle.H |  17 +
 5 files changed, 307 insertions(+), 450 deletions(-)

diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
index 7c2b383b83..b4a38f7ba0 100644
--- a/src/lagrangian/basic/Cloud/Cloud.C
+++ b/src/lagrangian/basic/Cloud/Cloud.C
@@ -42,7 +42,7 @@ template<class ParticleType>
 const Foam::scalar Foam::Cloud<ParticleType>::trackingRescueTolerance = 1e-3;
 
 template<class ParticleType>
-const Foam::scalar Foam::Cloud<ParticleType>::intersectionTolerance = 1e-6;
+const Foam::scalar Foam::Cloud<ParticleType>::intersectionTolerance = 1e-9;
 
 template<class ParticleType>
 const Foam::scalar Foam::Cloud<ParticleType>::planarCosAngle = (1 - 1e-6);
@@ -52,226 +52,143 @@ const Foam::scalar Foam::Cloud<ParticleType>::planarCosAngle = (1 - 1e-6);
 
 // Calculate using face properties only
 template<class ParticleType>
-Foam::label Foam::Cloud<ParticleType>::isIntrudingFace
-(
-    const label cellI,
-    const label f0I,
-    const label f1I
-) const
+bool Foam::Cloud<ParticleType>::isConcaveCell(const label cellI) const
 {
+    const cellList& cells = pMesh().cells();
+
+    const labelList& fOwner = pMesh().faceOwner();
     const vectorField& fAreas = pMesh().faceAreas();
     const pointField& fCentres = pMesh().faceCentres();
-    const labelList& fOwner = pMesh().faceOwner();
     const pointField& points = pMesh().points();
 
-    // Get f0 centre
-    const point& fc0 = fCentres[f0I];
+    const cell& cFaces = cells[cellI];
 
-    const face& f0 = pMesh().faces()[f0I];
-
-    const face& f1 = pMesh().faces()[f1I];
-
-    vector fn0 = fAreas[f0I];
-    fn0 /= (mag(fn0) + VSMALL);
-
-    // Flip normal if required so that it is always pointing out of
-    // the cell
-    if (fOwner[f0I] != cellI)
+    forAll(cFaces, i)
     {
-        fn0 *= -1;
-    }
+        label f0I = cFaces[i];
 
-    // Is any vertex of f1 on wrong side of the plane of f0?
-
-    forAll(f1, f1pI)
-    {
-        label ptI = f1[f1pI];
-
-        // Skip points that are shared between f1 and f0
-        if (findIndex(f0, ptI) > -1)
+        forAll(cFaces, j)
         {
-            continue;
-        }
+            if (j != i)
+            {
+                label f1I = cFaces[j];
 
-        const point& pt = points[ptI];
+                // Get f0 centre
+                const point& fc0 = fCentres[f0I];
 
-        // Normal distance from fc0 to pt
-        scalar d = ((fc0 - pt) & fn0);
+                const face& f0 = pMesh().faces()[f0I];
 
-        if (d < 0)
-        {
-            // Proper concave: f1 on wrong side of f0
-            return 2;
+                const face& f1 = pMesh().faces()[f1I];
+
+                vector fn0 = fAreas[f0I];
+                fn0 /= (mag(fn0) + VSMALL);
+
+                // Flip normal if required so that it is always pointing out of
+                // the cell
+                if (fOwner[f0I] != cellI)
+                {
+                    fn0 *= -1;
+                }
+
+                // Is any vertex of f1 on wrong side of the plane of f0?
+
+                forAll(f1, f1pI)
+                {
+                    label ptI = f1[f1pI];
+
+                    // Skip points that are shared between f1 and f0
+                    if (findIndex(f0, ptI) > -1)
+                    {
+                        continue;
+                    }
+
+                    const point& pt = points[ptI];
+
+                    // If the cell is concave, the point on f1 will be
+                    // on the outside of the plane of f0, defined by
+                    // its centre and normal, and the angle between
+                    // (fc0 -pt) and fn0 will be greater than 90
+                    // degrees, so the dot product will be negative.
+
+                    scalar d = ((fc0 - pt) & fn0);
+
+                    if (d < 0)
+                    {
+                        return true;
+                    }
+                }
+
+                // Check for co-planar faces, which are also treated
+                // as concave, as they are not strictly convex.
+
+                vector fn1 = fAreas[f1I];
+                fn1 /= (mag(fn1) + VSMALL);
+
+                // Flip normal if required so that it is always pointing out of
+                // the cell
+                if (fOwner[f1I] != cellI)
+                {
+                    fn1 *= -1;
+                }
+
+                if ((fn0 & fn1) > planarCosAngle)
+                {
+                    // Planar face
+                    return true;
+                }
+            }
         }
     }
 
-    // Could be a convex cell, but check angle for planar faces.
-
-    vector fn1 = fAreas[f1I];
-    fn1 /= (mag(fn1) + VSMALL);
-
-    // Flip normal if required so that it is always pointing out of
-    // the cell
-    if (fOwner[f1I] != cellI)
-    {
-        fn1 *= -1;
-    }
-
-    if ((fn0 & fn1) > planarCosAngle)
-    {
-        // Planar face
-        return 1;
-    }
-
-    return 0;
+    return false;
 }
 
 
 template<class ParticleType>
-void Foam::Cloud<ParticleType>::calcIntrudingFaces() const
+void Foam::Cloud<ParticleType>::calcConcaveCells() const
 {
-    const labelList& fOwner = pMesh().faceOwner();
     const cellList& cells = pMesh().cells();
 
-    intrudesIntoOwnerPtr_.reset(new PackedBoolList(pMesh().nFaces()));
-    PackedBoolList& intrudesIntoOwner = intrudesIntoOwnerPtr_();
-    intrudesIntoNeighbourPtr_.reset(new PackedBoolList(pMesh().nFaces()));
-    PackedBoolList& intrudesIntoNeighbour = intrudesIntoNeighbourPtr_();
-
-    PackedBoolList cellsWithProblemFaces(cells.size());
-    DynamicList<label> planarFaces;
-    DynamicList<label> intrudingFaces;
-
-
-    label nIntruded = 0;
-    label nPlanar = 0;
+    concaveCellPtr_.reset(new PackedBoolList(pMesh().nCells()));
+    PackedBoolList& concaveCell = concaveCellPtr_();
 
     forAll(cells, cellI)
     {
-        const cell& cFaces = cells[cellI];
-
-        forAll(cFaces, i)
+        if (isConcaveCell(cellI))
         {
-            label f0 = cFaces[i];
-
-            bool own0 = (fOwner[f0] == cellI);
-
-            // Now check that all other faces of the cell are on the 'inside'
-            // of the face.
-
-            bool intrudes = false;
-
-            forAll(cFaces, j)
-            {
-                if (j != i)
-                {
-                    label state = isIntrudingFace
-                    (
-                        cellI,
-                        f0,
-                        cFaces[j]
-                    );
-
-                    if (state == 1)
-                    {
-                        // planar
-                        planarFaces.append(f0);
-
-                        intrudes = true;
-                        nPlanar++;
-                        break;
-                    }
-                    else if (state == 2)
-                    {
-                        // concave
-                        intrudingFaces.append(f0);
-
-                        intrudes = true;
-                        nIntruded++;
-                        break;
-                    }
-                }
-            }
-
-            if (intrudes)
-            {
-                cellsWithProblemFaces[cellI] = 1;
-
-                if (own0)
-                {
-                    intrudesIntoOwner[f0] = 1;
-                }
-                else
-                {
-                    intrudesIntoNeighbour[f0] = 1;
-                }
-            }
-
-            // intrudesIntoOwner[f0] = 1;
-            // intrudesIntoNeighbour[f0] = 1;
+            concaveCell[cellI] = 1;
         }
+
+        // Force all cells to be decomposed
+        concaveCell[cellI] = 1;
     }
 
-    // if (debug)
     {
-        Pout<< "Cloud<ParticleType>::calcIntrudingFaces() :"
-            << " overall faces in mesh : "
-            << pMesh().nFaces() << nl
-            << "    of these planar : "
-            << nPlanar << nl
-            << "    of these intruding into their cell (concave) : "
-            << nIntruded << endl;
+        // Write cells that are a problem to file
 
-        // Info<< "Hard coded picking up every face." << endl;
-
-        // Write the faces and cells that are a problem to file to be
-        // compared and made into sets
+        DynamicList<label> tmpConcaveCells;
 
+        forAll(cells, cellI)
         {
-            DynamicList<label> cellsWithIntrudingFaces;
-
-            forAll(cells, cellI)
+            if (concaveCell[cellI])
             {
-                if (cellsWithProblemFaces[cellI])
-                {
-                    cellsWithIntrudingFaces.append(cellI);
-                }
+                tmpConcaveCells.append(cellI);
             }
-
-            fileName fName = pMesh().time().path()/"cellsWithProblemFaces";
-
-            Pout<< "    Writing " << fName.name() << endl;
-
-            OFstream file(fName);
-
-            file << cellsWithIntrudingFaces;
-
-            file.flush();
         }
 
-        {
-            fileName fName = pMesh().time().path()/"planarFaces";
+        Pout<< "Cloud<ParticleType>::calcConcaveCells() :"
+            << " overall cells in mesh : " << pMesh().nCells() << nl
+            << "    of these concave : " << tmpConcaveCells.size() << endl;
 
-            Pout<< "    Writing " << fName.name() << endl;
+        fileName fName = pMesh().time().path()/"concaveCells";
 
-            OFstream file(fName);
+        Pout<< "    Writing " << fName.name() << endl;
 
-            file << planarFaces;
+        OFstream file(fName);
 
-            file.flush();
-        }
+        file << tmpConcaveCells;
 
-        {
-            fileName fName = pMesh().time().path()/"intrudingFaces";
-
-            Pout<< "    Writing " << fName.name() << endl;
-
-            OFstream file(fName);
-
-            file << intrudingFaces;
-
-            file.flush();
-        }
+        file.flush();
     }
 }
 
@@ -318,8 +235,7 @@ Foam::Cloud<ParticleType>::Cloud
 template<class ParticleType>
 void Foam::Cloud<ParticleType>::clearOut()
 {
-    intrudesIntoOwnerPtr_.clear();
-    intrudesIntoNeighbourPtr_.clear();
+    concaveCellPtr_.clear();
     labels_.clearStorage();
 }
 
@@ -327,26 +243,15 @@ void Foam::Cloud<ParticleType>::clearOut()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class ParticleType>
-const Foam::PackedBoolList& Foam::Cloud<ParticleType>::intrudesIntoOwner()
+const Foam::PackedBoolList& Foam::Cloud<ParticleType>::concaveCell()
 const
 {
-    if (!intrudesIntoOwnerPtr_.valid())
+    if (!concaveCellPtr_.valid())
     {
-        calcIntrudingFaces();
+        calcConcaveCells();
     }
-    return intrudesIntoOwnerPtr_();
-}
 
-
-template<class ParticleType>
-const Foam::PackedBoolList& Foam::Cloud<ParticleType>::intrudesIntoNeighbour()
-const
-{
-    if (!intrudesIntoNeighbourPtr_.valid())
-    {
-        calcIntrudingFaces();
-    }
-    return intrudesIntoNeighbourPtr_();
+    return concaveCellPtr_();
 }
 
 
diff --git a/src/lagrangian/basic/Cloud/Cloud.H b/src/lagrangian/basic/Cloud/Cloud.H
index 11a0809c8f..e30204abd6 100644
--- a/src/lagrangian/basic/Cloud/Cloud.H
+++ b/src/lagrangian/basic/Cloud/Cloud.H
@@ -78,10 +78,8 @@ class Cloud
 
         const bool concaveCheck_;
 
-        //- Does face intrude into owner cell
-        mutable autoPtr<PackedBoolList> intrudesIntoOwnerPtr_;
-        //- Does face intrude into neighbour cell
-        mutable autoPtr<PackedBoolList> intrudesIntoNeighbourPtr_;
+        //- Is the cell concave
+        mutable autoPtr<PackedBoolList> concaveCellPtr_;
 
         //- Overall count of particles ever created. Never decreases.
         mutable label particleCount_;
@@ -92,20 +90,11 @@ class Cloud
 
     // Private member functions
 
-        //- Helper: is face f1 on wrong side of f0 (both faces are on same cell)
-        //  0 : convex
-        //  1 : planar to within planarCosAngle
-        //  2 : concave
-        label isIntrudingFace
-        (
-            const label cellI,
-            const label f0,
-            const label f1
-        ) const;
+        //- Determine if a particular cell is concave
+        bool isConcaveCell(const label cellI) const;
 
         //- Analyse mesh for faces not being convex.
-        //  Sets intrudesIntoOwner_,intrudesIntoNeighbour_.
-        void calcIntrudingFaces() const;
+        void calcConcaveCells() const;
 
         //- Initialise cloud on IO constructor
         void initCloud(const bool checkClass);
@@ -229,12 +218,8 @@ public:
                 return IDLList<ParticleType>::size();
             };
 
-            //- Per face whether it intrudes into the owner cell.(demand driven)
-            const PackedBoolList& intrudesIntoOwner() const;
-
-            //- Per face whether it intrudes into the neighbour cell. (demand
-            //  driven.
-            const PackedBoolList& intrudesIntoNeighbour() const;
+            //- Whether each cell is concave (demand driven data)
+            const PackedBoolList& concaveCell() const;
 
 
             // Iterators
diff --git a/src/lagrangian/basic/Cloud/CloudIO.C b/src/lagrangian/basic/Cloud/CloudIO.C
index 1eb35e09d0..91e5808e6a 100644
--- a/src/lagrangian/basic/Cloud/CloudIO.C
+++ b/src/lagrangian/basic/Cloud/CloudIO.C
@@ -90,11 +90,6 @@ Foam::Cloud<ParticleType>::Cloud
     particleCount_(0)
 {
     initCloud(checkClass);
-
-    // if (concaveCheck_)
-    // {
-    //     intrudesIntoOwner();
-    // }
 }
 
 
diff --git a/src/lagrangian/basic/Particle/Particle.C b/src/lagrangian/basic/Particle/Particle.C
index dc0c383ea4..81b0532c97 100644
--- a/src/lagrangian/basic/Particle/Particle.C
+++ b/src/lagrangian/basic/Particle/Particle.C
@@ -48,67 +48,10 @@ void Foam::Particle<ParticleType>::findFaces
 
     faceList.clear();
 
-    // forAll(faces, i)
-    // {
-    //     label facei = faces[i];
-    //     scalar lam = lambda(C, position, facei);
-
-    //     if ((lam > 0) && (lam < 1.0))
-    //     {
-    //         faceList.append(facei);
-    //     }
-    // }
-
-    vector deltaPosition = position - C;
-    scalar deltaPositionMag = mag(deltaPosition);
-    scalar lam;
-
     forAll(faces, i)
     {
         label facei = faces[i];
-
-        lam = GREAT;
-
-        if (cloud_.concaveCheck_)
-        {
-            // Check for intruding face
-            // ~~~~~~~~~~~~~~~~~~~~~~~~
-
-            label isOwner = (celli_ == mesh.faceOwner()[facei]);
-
-            if
-            (
-                (isOwner && cloud_.intrudesIntoOwner()[facei])
-             || (!isOwner && cloud_.intrudesIntoNeighbour()[facei])
-            )
-            {
-                // Concave face. Use exact intersection.
-
-                pointHit inter = mesh.faces()[facei].intersection
-                (
-                    position_,
-                    deltaPosition,
-                    mesh.faceCentres()[facei],
-                    mesh.points(),
-                    intersection::HALF_RAY,
-                    Cloud<ParticleType>::intersectionTolerance
-                );
-
-                if (inter.hit())
-                {
-                    lam = inter.distance()/deltaPositionMag;
-                }
-            }
-            else
-            {
-                // Convex face. Can use lambda.
-                lam = lambda(C, position, facei);
-            }
-        }
-        else
-        {
-            lam = lambda(C, position, facei);
-        }
+        scalar lam = lambda(C, position, facei);
 
         if ((lam > 0) && (lam < 1.0))
         {
@@ -145,6 +88,163 @@ void Foam::Particle<ParticleType>::findFaces
 }
 
 
+template<class ParticleType>
+void Foam::Particle<ParticleType>::trackToFaceByDecomposition
+(
+    scalar& trackFraction,
+    const vector& endPosition
+)
+{
+    const polyMesh& mesh = cloud_.polyMesh_;
+    const labelList& faces = mesh.cells()[celli_];
+
+    vector deltaPosition = endPosition - position_;
+
+    scalar lambdaMin = GREAT;
+    label workingFace = -1;
+    scalar hitLambda = GREAT;
+
+    forAll(faces, i)
+    {
+        label facei = faces[i];
+
+        pointHit inter = mesh.faces()[facei].intersection
+        (
+            position_,
+            deltaPosition,
+            mesh.faceCentres()[facei],
+            mesh.points(),
+            intersection::HALF_RAY,
+            Cloud<ParticleType>::intersectionTolerance
+        );
+
+        if (inter.hit())
+        {
+            hitLambda = inter.distance();
+
+            // There is a problem if the particle is essentially on
+            // the face, it will flip.
+
+            if (hitLambda <= 1.0 && hitLambda < lambdaMin)
+            {
+                Info<< origProc_ << " "
+                    << origId_ << " "
+                    << position_ << " "
+                    << endPosition << " "
+                    << trackFraction << " "
+                    << stepFraction_ << " "
+                    << hitLambda << " "
+                    << inter.hitPoint() << " "
+                    << celli_ << " "
+                    << facei << endl;
+
+                lambdaMin = hitLambda;
+                workingFace = facei;
+            }
+        }
+    }
+
+    if (workingFace == -1)
+    {
+        // inside cell
+
+        facei_ = workingFace;
+        trackFraction = 1.0;
+        position_ = endPosition;
+    }
+    else
+    {
+        if (lambdaMin < 0.0)
+        {
+            FatalErrorIn("Particle::trackToFaceByDecomposition")
+                << "Decomposed tracking failure, lambdaMin = " << lambdaMin
+                << nl << abort(FatalError);
+        }
+
+        trackFraction = lambdaMin;
+
+        position_ += trackFraction*deltaPosition;
+
+        // What about softImpact?
+    }
+}
+
+
+template<class ParticleType>
+void Foam::Particle<ParticleType>::trackToFaceByPlanes
+(
+    scalar& trackFraction,
+    const vector& endPosition
+)
+{
+    facei_ = -1;
+
+    DynamicList<label>& faces = cloud_.labels_;
+    findFaces(endPosition, faces);
+
+    if (faces.empty())
+    {
+        // inside cell
+
+        trackFraction = 1.0;
+        position_ = endPosition;
+    }
+    else
+    {
+        // hit face
+
+        scalar lambdaMin = GREAT;
+
+        if (faces.size() == 1)
+        {
+            lambdaMin = lambda(position_, endPosition, faces[0], stepFraction_);
+            facei_ = faces[0];
+        }
+        else
+        {
+            forAll(faces, i)
+            {
+                scalar lam =
+                    lambda(position_, endPosition, faces[i], stepFraction_);
+
+                if (lam < lambdaMin)
+                {
+                    lambdaMin = lam;
+                    facei_ = faces[i];
+                }
+            }
+        }
+
+        // For warped faces the particle can be 'outside' the cell.
+        // This will yield a lambda larger than 1, or smaller than 0
+        // For values < 0, the particle travels away from the cell
+        // and we don't move the particle, only change cell.
+        // For values larger than 1, we move the particle to endPosition only.
+        if (lambdaMin > 0.0)
+        {
+            if (lambdaMin <= 1.0)
+            {
+                trackFraction = lambdaMin;
+                position_ += trackFraction*(endPosition - position_);
+            }
+            else
+            {
+                trackFraction = 1.0;
+                position_ = endPosition;
+            }
+        }
+        else if (static_cast<ParticleType&>(*this).softImpact())
+        {
+            // Soft-sphere particles can travel outside the domain
+            // but we don't use lambda since this the particle
+            // is going away from face
+            trackFraction = 1.0;
+            position_ = endPosition;
+        }
+    }
+}
+
+
 template<class ParticleType>
 template<class TrackData>
 void Foam::Particle<ParticleType>::prepareForParallelTransfer
@@ -275,7 +375,6 @@ Foam::label Foam::Particle<ParticleType>::track
 }
 
 
-
 template<class ParticleType>
 Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
 {
@@ -283,6 +382,7 @@ Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
     return track(endPosition, dummyTd);
 }
 
+
 template<class ParticleType>
 template<class TrackData>
 Foam::scalar Foam::Particle<ParticleType>::trackToFace
@@ -293,154 +393,25 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
 {
     const polyMesh& mesh = cloud_.polyMesh_;
 
-    DynamicList<label>& faces = cloud_.labels_;
-    findFaces(endPosition, faces);
-
-    facei_ = -1;
     scalar trackFraction = 0.0;
 
-    scalar lambdaMin = GREAT;
-
-    // Temporary addition to report on infinite loops ->
-    List<scalar> lambdas(faces.size());
-
-    forAll(faces, i)
+    if (cloud_.concaveCell()[celli_])
     {
-        lambdas[i] = lambda(position_, endPosition, faces[i], stepFraction_);
+        trackToFaceByDecomposition(trackFraction, endPosition);
     }
-    // <- Temporary addition to report on infinite loops
-
-    if (faces.empty())  // inside cell
+    else
     {
-        trackFraction = 1.0;
-        position_ = endPosition;
+        trackToFaceByPlanes(trackFraction, endPosition);
     }
-    else // hit face
+
+    if (facei_ > -1)
     {
-        // scalar lambdaMin = GREAT;
+        // A face was hit
 
-        if (faces.size() == 1)
+        if (cloud_.internalFace(facei_))
         {
-            lambdaMin = lambda(position_, endPosition, faces[0], stepFraction_);
-            facei_ = faces[0];
-        }
-        else
-        {
-            // // If the particle has to cross more than one cell to reach the
-            // // endPosition, we check which way to go.
-            // // If one of the faces is a boundary face and the particle is
-            // // outside, we choose the boundary face.
-            // The particle is outside if one of the lambda's is > 1 or < 0
+            // Internal face, change cell
 
-            scalar fallBackLambdaMin = GREAT;
-            label fallBackFacei = -1;
-
-            vector deltaPosition = endPosition - position_;
-            scalar deltaPositionMag = mag(deltaPosition);
-            scalar hitLambda = GREAT;
-
-            forAll(faces, i)
-            {
-                scalar lam =
-                    lambda(position_, endPosition, faces[i], stepFraction_);
-
-                // Calculate fallback
-                if (lam < fallBackLambdaMin)
-                {
-                    fallBackLambdaMin = lam;
-                    fallBackFacei = faces[i];
-                }
-
-                if (cloud_.concaveCheck_)
-                {
-                    // Check for intruding face
-                    // ~~~~~~~~~~~~~~~~~~~~~~~~
-
-                    label isOwner = (celli_ == mesh.faceOwner()[faces[i]]);
-
-                    if
-                    (
-                        (isOwner && cloud_.intrudesIntoOwner()[faces[i]])
-                     || (!isOwner && cloud_.intrudesIntoNeighbour()[faces[i]])
-                    )
-                    {
-                        // Concave face. Use exact intersection.
-
-                        pointHit inter = mesh.faces()[faces[i]].intersection
-                        (
-                            position_,
-                            deltaPosition,
-                            mesh.faceCentres()[faces[i]],
-                            mesh.points(),
-                            intersection::HALF_RAY,
-                            Cloud<ParticleType>::intersectionTolerance
-                        );
-
-                        if (inter.hit())
-                        {
-                            hitLambda = inter.distance()/deltaPositionMag;
-
-                            if (hitLambda < lambdaMin)
-                            {
-                                lambdaMin = hitLambda;
-                                facei_ = faces[i];
-                            }
-                        }
-                    }
-                    else
-                    {
-                        // Convex face. Can use lambda.
-
-                        if (lam < lambdaMin)
-                        {
-                            lambdaMin = lam;
-                            facei_ = faces[i];
-                        }
-                    }
-                }
-            }
-
-            // Fall back to the original algorithm if no
-            // intersection was found
-            if (facei_ == -1)
-            {
-                lambdaMin = fallBackLambdaMin;
-                facei_ = fallBackFacei;
-            }
-        }
-
-        bool internalFace = cloud_.internalFace(facei_);
-
-        // For warped faces the particle can be 'outside' the cell.
-        // This will yield a lambda larger than 1, or smaller than 0
-        // For values < 0, the particle travels away from the cell
-        // and we don't move the particle, only change cell.
-        // For values larger than 1, we move the particle to endPosition only.
-        if (lambdaMin > 0.0)
-        {
-            if (lambdaMin <= 1.0)
-            {
-                trackFraction = lambdaMin;
-                position_ += trackFraction*(endPosition - position_);
-            }
-            else
-            {
-                trackFraction = 1.0;
-                position_ = endPosition;
-            }
-        }
-        else if (static_cast<ParticleType&>(*this).softImpact())
-        {
-            // Soft-sphere particles can travel outside the domain
-            // but we don't use lambda since this the particle
-            // is going away from face
-            trackFraction = 1.0;
-            position_ = endPosition;
-        }
-
-        // change cell
-        if (internalFace) // Internal face
-        {
             if (celli_ == mesh.faceOwner()[facei_])
             {
                 celli_ = mesh.faceNeighbour()[facei_];
@@ -451,11 +422,9 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
             }
             else
             {
-                FatalErrorIn
-                (
-                    "Particle::trackToFace(const vector&, TrackData&)"
-                )<< "addressing failure" << nl
-                 << abort(FatalError);
+                FatalErrorIn("Particle::trackToFace(const vector&, TrackData&)")
+                    << "addressing failure" << nl
+                    << abort(FatalError);
             }
         }
         else
@@ -526,31 +495,17 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
 
     if (trackFraction < Cloud<ParticleType>::minValidTrackFraction)
     {
-        Pout<< "tracking error "
-            << origId_ << " "
-            << origProc_ << " "
-            << position_ << " "
-            << endPosition << " "
-            << position_ + 1000*(endPosition - position_) << " "
-            << position_ - 1000*(endPosition - position_) << " "
-            << lambdaMin << " "
-            << trackFraction << " "
-            << stepFraction_ << " "
-            << celli_ << " "
-            << facei_ << " "
-            << cloud_.intrudesIntoOwner()[facei_] << " "
-            << cloud_.intrudesIntoNeighbour()[facei_] << " "
-            << faces << " "
-            << lambdas
-            << endl;
-
-        // if (stepFraction_  < SMALL)
-        // {
-        //     Pout<< "Intervene to mark face" << endl;
-
-        //     cloud_.intrudesIntoOwnerPtr_()[facei_] = 1;
-        //     cloud_.intrudesIntoNeighbourPtr_()[facei_] = 1;
-        // }
+        // Pout<< "tracking error "
+        //     << origId_ << " "
+        //     << origProc_ << " "
+        //     << position_ << " "
+        //     << endPosition << " "
+        //     << trackFraction << " "
+        //     << stepFraction_ << " "
+        //     << celli_ << " "
+        //     << facei_ << " "
+        //     << cloud_.concaveCell()[celli_]
+        //     << endl;
 
         const point& cc = mesh.cellCentres()[celli_];
         position_ +=
diff --git a/src/lagrangian/basic/Particle/Particle.H b/src/lagrangian/basic/Particle/Particle.H
index a1077a604e..cf44cc1198 100644
--- a/src/lagrangian/basic/Particle/Particle.H
+++ b/src/lagrangian/basic/Particle/Particle.H
@@ -173,6 +173,23 @@ protected:
             DynamicList<label>& faceList
         ) const;
 
+        //- Track to cell face using face decomposition, used for
+        //  concave cells.
+        void trackToFaceByDecomposition
+        (
+            scalar& trackFraction,
+            const vector& endPosition
+        );
+
+        //- Track to cell face using the infinite planes of the faces.
+        //  The cell *must* be convex.
+        void trackToFaceByPlanes
+        (
+            scalar& trackFraction,
+            const vector& endPosition
+        );
+
+
 
     // Patch interactions
 

From 46b240ba86346c8833546b21132e96848ed70004 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Thu, 8 Oct 2009 14:28:30 +0100
Subject: [PATCH 03/16] Using a new, more careful method when a cell is
 concave, not relying on plane crossings for any actual motion or face
 intersection decisions.

---
 src/lagrangian/basic/Cloud/Cloud.C       | 105 +++-
 src/lagrangian/basic/Cloud/Cloud.H       |  13 +
 src/lagrangian/basic/Particle/Particle.C | 617 +++++++++++++++--------
 src/lagrangian/basic/Particle/Particle.H |  28 +-
 4 files changed, 533 insertions(+), 230 deletions(-)

diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
index b4a38f7ba0..5d6cc5705a 100644
--- a/src/lagrangian/basic/Cloud/Cloud.C
+++ b/src/lagrangian/basic/Cloud/Cloud.C
@@ -42,7 +42,7 @@ template<class ParticleType>
 const Foam::scalar Foam::Cloud<ParticleType>::trackingRescueTolerance = 1e-3;
 
 template<class ParticleType>
-const Foam::scalar Foam::Cloud<ParticleType>::intersectionTolerance = 1e-9;
+const Foam::scalar Foam::Cloud<ParticleType>::intersectionTolerance = SMALL;
 
 template<class ParticleType>
 const Foam::scalar Foam::Cloud<ParticleType>::planarCosAngle = (1 - 1e-6);
@@ -59,37 +59,52 @@ bool Foam::Cloud<ParticleType>::isConcaveCell(const label cellI) const
     const labelList& fOwner = pMesh().faceOwner();
     const vectorField& fAreas = pMesh().faceAreas();
     const pointField& fCentres = pMesh().faceCentres();
+    const pointField& cCentres = pMesh().cellCentres();
     const pointField& points = pMesh().points();
 
+    PackedBoolList& intrudesIntoOwner = intrudesIntoOwnerPtr_();
+    PackedBoolList& intrudesIntoNeighbour = intrudesIntoNeighbourPtr_();
+
     const cell& cFaces = cells[cellI];
 
+    bool concave = false;
+
     forAll(cFaces, i)
     {
         label f0I = cFaces[i];
 
+        // Get f0 centre
+        const point& fc0 = fCentres[f0I];
+
+        const face& f0 = pMesh().faces()[f0I];
+
+        vector fn0 = fAreas[f0I];
+        fn0 /= (mag(fn0) + VSMALL);
+
+        // Flip normal if required so that it is always pointing out of
+        // the cell
+        if (fOwner[f0I] != cellI)
+        {
+            fn0 *= -1;
+        }
+
+        // Check if the cell centre is visible from the plane of the face
+
+        if (((fc0 - cCentres[cellI]) & fn0) < 0)
+        {
+            Pout<< "Face " << f0I
+                << " is not visible from the centre of cell " << cellI
+                << endl;
+        }
+
         forAll(cFaces, j)
         {
             if (j != i)
             {
                 label f1I = cFaces[j];
 
-                // Get f0 centre
-                const point& fc0 = fCentres[f0I];
-
-                const face& f0 = pMesh().faces()[f0I];
-
                 const face& f1 = pMesh().faces()[f1I];
 
-                vector fn0 = fAreas[f0I];
-                fn0 /= (mag(fn0) + VSMALL);
-
-                // Flip normal if required so that it is always pointing out of
-                // the cell
-                if (fOwner[f0I] != cellI)
-                {
-                    fn0 *= -1;
-                }
-
                 // Is any vertex of f1 on wrong side of the plane of f0?
 
                 forAll(f1, f1pI)
@@ -114,7 +129,17 @@ bool Foam::Cloud<ParticleType>::isConcaveCell(const label cellI) const
 
                     if (d < 0)
                     {
-                        return true;
+                        // Concave face
+                        if (fOwner[f0I] == cellI)
+                        {
+                            intrudesIntoOwner[f0I] = 1;
+                        }
+                        else
+                        {
+                            intrudesIntoNeighbour[f0I] = 1;
+                        }
+
+                        concave = true;
                     }
                 }
 
@@ -134,13 +159,22 @@ bool Foam::Cloud<ParticleType>::isConcaveCell(const label cellI) const
                 if ((fn0 & fn1) > planarCosAngle)
                 {
                     // Planar face
-                    return true;
+                    if (fOwner[f0I] == cellI)
+                    {
+                        intrudesIntoOwner[f0I] = 1;
+                    }
+                    else
+                    {
+                        intrudesIntoNeighbour[f0I] = 1;
+                    }
+
+                    concave = true;
                 }
             }
         }
     }
 
-    return false;
+    return concave;
 }
 
 
@@ -152,15 +186,16 @@ void Foam::Cloud<ParticleType>::calcConcaveCells() const
     concaveCellPtr_.reset(new PackedBoolList(pMesh().nCells()));
     PackedBoolList& concaveCell = concaveCellPtr_();
 
+    intrudesIntoOwnerPtr_.reset(new PackedBoolList(pMesh().nFaces()));
+    intrudesIntoNeighbourPtr_.reset(new PackedBoolList(pMesh().nFaces()));
+
+
     forAll(cells, cellI)
     {
         if (isConcaveCell(cellI))
         {
             concaveCell[cellI] = 1;
         }
-
-        // Force all cells to be decomposed
-        concaveCell[cellI] = 1;
     }
 
     {
@@ -236,6 +271,8 @@ template<class ParticleType>
 void Foam::Cloud<ParticleType>::clearOut()
 {
     concaveCellPtr_.clear();
+    intrudesIntoOwnerPtr_.clear();
+    intrudesIntoNeighbourPtr_.clear();
     labels_.clearStorage();
 }
 
@@ -255,6 +292,30 @@ const
 }
 
 
+template<class ParticleType>
+const Foam::PackedBoolList& Foam::Cloud<ParticleType>::intrudesIntoOwner()
+const
+{
+    if (!intrudesIntoOwnerPtr_.valid())
+    {
+        calcConcaveCells();
+    }
+    return intrudesIntoOwnerPtr_();
+}
+
+
+template<class ParticleType>
+const Foam::PackedBoolList& Foam::Cloud<ParticleType>::intrudesIntoNeighbour()
+const
+{
+    if (!intrudesIntoNeighbourPtr_.valid())
+    {
+        calcConcaveCells();
+    }
+    return intrudesIntoNeighbourPtr_();
+}
+
+
 template<class ParticleType>
 Foam::label Foam::Cloud<ParticleType>::getNewParticleID() const
 {
diff --git a/src/lagrangian/basic/Cloud/Cloud.H b/src/lagrangian/basic/Cloud/Cloud.H
index e30204abd6..b151827d0a 100644
--- a/src/lagrangian/basic/Cloud/Cloud.H
+++ b/src/lagrangian/basic/Cloud/Cloud.H
@@ -81,6 +81,12 @@ class Cloud
         //- Is the cell concave
         mutable autoPtr<PackedBoolList> concaveCellPtr_;
 
+        //- Does face intrude into owner cell
+        mutable autoPtr<PackedBoolList> intrudesIntoOwnerPtr_;
+
+        //- Does face intrude into neighbour cell
+        mutable autoPtr<PackedBoolList> intrudesIntoNeighbourPtr_;
+
         //- Overall count of particles ever created. Never decreases.
         mutable label particleCount_;
 
@@ -221,6 +227,13 @@ public:
             //- Whether each cell is concave (demand driven data)
             const PackedBoolList& concaveCell() const;
 
+            //- Per face whether it intrudes into the owner cell.(demand driven)
+            const PackedBoolList& intrudesIntoOwner() const;
+
+            //- Per face whether it intrudes into the neighbour cell. (demand
+            //  driven.
+            const PackedBoolList& intrudesIntoNeighbour() const;
+
 
             // Iterators
 
diff --git a/src/lagrangian/basic/Particle/Particle.C b/src/lagrangian/basic/Particle/Particle.C
index 81b0532c97..b2dc2d3299 100644
--- a/src/lagrangian/basic/Particle/Particle.C
+++ b/src/lagrangian/basic/Particle/Particle.C
@@ -38,7 +38,7 @@ License
 template<class ParticleType>
 void Foam::Particle<ParticleType>::findFaces
 (
-    const vector& position,
+    const vector& endPosition,
     DynamicList<label>& faceList
 ) const
 {
@@ -51,7 +51,7 @@ void Foam::Particle<ParticleType>::findFaces
     forAll(faces, i)
     {
         label facei = faces[i];
-        scalar lam = lambda(C, position, facei);
+        scalar lam = lambda(C, endPosition, facei);
 
         if ((lam > 0) && (lam < 1.0))
         {
@@ -64,7 +64,7 @@ void Foam::Particle<ParticleType>::findFaces
 template<class ParticleType>
 void Foam::Particle<ParticleType>::findFaces
 (
-    const vector& position,
+    const vector& endPosition,
     const label celli,
     const scalar stepFraction,
     DynamicList<label>& faceList
@@ -78,7 +78,7 @@ void Foam::Particle<ParticleType>::findFaces
     forAll(faces, i)
     {
         label facei = faces[i];
-        scalar lam = lambda(C, position, facei, stepFraction);
+        scalar lam = lambda(C, endPosition, facei, stepFraction);
 
         if ((lam > 0) && (lam < 1.0))
         {
@@ -89,25 +89,140 @@ void Foam::Particle<ParticleType>::findFaces
 
 
 template<class ParticleType>
-void Foam::Particle<ParticleType>::trackToFaceByDecomposition
+template<class TrackData>
+void Foam::Particle<ParticleType>::trackToFaceConcave
 (
     scalar& trackFraction,
-    const vector& endPosition
+    const vector& endPosition,
+    TrackData& td
 )
 {
-    const polyMesh& mesh = cloud_.polyMesh_;
+    facei_ = -1;
+
+    const polyMesh& mesh = cloud_.pMesh();
     const labelList& faces = mesh.cells()[celli_];
+    const vector& C = mesh.cellCentres()[celli_];
 
-    vector deltaPosition = endPosition - position_;
+    DynamicList<label>& potentiallyCrossedFaces = cloud_.labels_;
+    findFaces(endPosition, potentiallyCrossedFaces);
 
-    scalar lambdaMin = GREAT;
-    label workingFace = -1;
-    scalar hitLambda = GREAT;
+    // Check all possible face crossings to see if they are actually
+    // crossed, determining if endPosition is outside the current cell.
+
+    if (potentiallyCrossedFaces.empty())
+    {
+        // endPosition is inside the cell
+
+        position_ = endPosition;
+        trackFraction = 1.0;
+        return;
+    }
+
+    label nFaceCrossings = 0;
+
+    forAll (potentiallyCrossedFaces, pCFI)
+    {
+        label facei = potentiallyCrossedFaces[pCFI];
+
+        pointHit inter = mesh.faces()[facei].intersection
+        (
+            C,
+            endPosition - C,
+            mesh.faceCentres()[facei],
+            mesh.points(),
+            intersection::HALF_RAY,
+            Cloud<ParticleType>::intersectionTolerance
+        );
+
+        if (inter.hit())
+        {
+            if (inter.distance() < 1.0)
+            {
+                // This face was actually crossed.
+
+                nFaceCrossings++;
+
+                if (inter.distance() < 0.0)
+                {
+                    FatalErrorIn("Particle::trackToFaceConcave")
+                        << "Negative lambda " << inter.distance()
+                        << "in potentiallyCrossedFaces test"
+                        << nl << abort(FatalError);
+                }
+            }
+        }
+    }
+
+    if (nFaceCrossings > 1)
+    {
+        Pout<< "In cell " << celli_ << " there were " << nFaceCrossings
+            << " face crossings detected tracking from cell centre to "
+            << " endPosition"
+            << endl;
+    }
+
+    if (nFaceCrossings % 2 == 0)
+    {
+        // Even number of face crossings, so the particle must end up
+        // still in the cell.
+
+        position_ = endPosition;
+        trackFraction = 1.0;
+        return;
+    }
+
+    // The particle *must* have left the cell.
+
+    // a) It may have crossed a face not yet identified by findFaces using the
+    //    cell centre to endPosition line, so the potentially crossed faces of
+    //    the position to endPosition line must be assessed.
+
+    // b) It may have been outside the cell in the first place, and, despite
+    //    trying to pick up more faces using a) the correct face to be crossed
+    //    is not knowable.  A best guess will be used, with the expectation that
+    //    the tracking in the destination cell will be able to recover form a
+    //    bad guess.
+
+    // For all face assessments, a full intersection test is required,
+    // as nothing can be assumed about the order of crossing the
+    // planes of faces.
 
     forAll(faces, i)
     {
         label facei = faces[i];
 
+        scalar lam = lambda(position_, endPosition, facei);
+
+        if ((lam > 0) && (lam < 1.0))
+        {
+            if (findIndex(potentiallyCrossedFaces, facei) == -1)
+            {
+                // A new potential face crossing, previously missed,
+                // has been found
+
+                Pout<< "New face crossing " << facei
+                    << " of cell " << celli_ << " found"
+                    << endl;
+
+                potentiallyCrossedFaces.append(facei);
+            }
+        }
+    }
+
+    vector deltaPosition = endPosition - position_;
+
+    scalar tmpLambda = GREAT;
+    scalar correctLambda = GREAT;
+
+    forAll(potentiallyCrossedFaces, pCFI)
+    {
+        label facei = potentiallyCrossedFaces[pCFI];
+
+        // Use exact intersection.
+
+        // TODO: A correction is required for moving meshes to
+        // calculate the correct lambda value.
+
         pointHit inter = mesh.faces()[facei].intersection
         (
             position_,
@@ -120,61 +235,130 @@ void Foam::Particle<ParticleType>::trackToFaceByDecomposition
 
         if (inter.hit())
         {
-            hitLambda = inter.distance();
+            tmpLambda = inter.distance();
 
-            // There is a problem if the particle is essentially on
-            // the face, it will flip.
-
-            if (hitLambda <= 1.0 && hitLambda < lambdaMin)
+            if
+            (
+                tmpLambda < 1.0
+             && tmpLambda < correctLambda
+            )
             {
-                Info<< origProc_ << " "
-                    << origId_ << " "
-                    << position_ << " "
-                    << endPosition << " "
-                    << trackFraction << " "
-                    << stepFraction_ << " "
-                    << hitLambda << " "
-                    << inter.hitPoint() << " "
-                    << celli_ << " "
-                    << facei << endl;
+                // This face is crossed before any other that has
+                // been found so far
 
-                lambdaMin = hitLambda;
-                workingFace = facei;
+                correctLambda = tmpLambda;
+                facei_ = facei;
             }
         }
     }
 
-    if (workingFace == -1)
+    if (facei_ > -1)
     {
-        // inside cell
-
-        facei_ = workingFace;
-        trackFraction = 1.0;
-        position_ = endPosition;
+        trackFraction = correctLambda;
+        position_ += trackFraction*(endPosition - position_);
     }
     else
     {
-        if (lambdaMin < 0.0)
+        // 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.
+
+        DynamicList<scalar> tmpLambdas;
+        DynamicList<label> tmpLambdaFaceIs;
+
+        forAll(faces, i)
         {
-            FatalErrorIn("Particle::trackToFaceByDecomposition")
-                << "Decomposed tracking failure, lambdaMin = " << lambdaMin
-                << nl << abort(FatalError);
+            label facei = faces[i];
+
+            if (findIndex(potentiallyCrossedFaces, facei) == -1)
+            {
+                // 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);
+                }
+            }
         }
 
-        trackFraction = lambdaMin;
+        // The closest lambda value (less than zero) or (greater than or
+        // equal to one), and corresponding face.
+        scalar closestLambda = GREAT;
+        scalar closestLambdaDifference = GREAT;
+        label closestLambdaDifferenceFaceI = -1;
 
-        position_ += trackFraction*deltaPosition;
+        forAll(tmpLambdaFaceIs, i)
+        {
+            label facei = tmpLambdaFaceIs[i];
+            scalar tmpLambda = tmpLambdas[i];
 
-        // What about softImpact?
+            if (tmpLambda < 0)
+            {
+                if (mag(tmpLambda) < closestLambdaDifference)
+                {
+                    closestLambdaDifference = mag(tmpLambda);
+                    closestLambda = tmpLambda;
+                    closestLambdaDifferenceFaceI = facei;
+                }
+            }
+            else if (tmpLambda >= 1.0)
+            {
+                if ((tmpLambda - 1.0) < closestLambdaDifference)
+                {
+                    closestLambdaDifference = tmpLambda - 1.0;
+                    closestLambda = tmpLambda;
+                    closestLambdaDifferenceFaceI = facei;
+                }
+            }
+        }
+
+        if (closestLambdaDifferenceFaceI > -1)
+        {
+            // A face has been found to cross.  Not moving the
+            // particle at all, but the option is there to decide do
+            // it, as the sign of closestLambda is available to tell
+            // if the particle is heading towards, or away from the
+            // face.
+
+            facei_ = closestLambdaDifferenceFaceI;
+        }
+        else
+        {
+            Pout<< "Tracking failure in concave cell - didn't find a "
+                << "best guess face" << endl;
+        }
     }
+
+    faceAction(trackFraction, endPosition, td);
 }
 
 
 template<class ParticleType>
-void Foam::Particle<ParticleType>::trackToFaceByPlanes
+template<class TrackData>
+void Foam::Particle<ParticleType>::trackToFaceConvex
 (
     scalar& trackFraction,
-    const vector& endPosition
+    const vector& endPosition,
+    TrackData& td
 )
 {
     facei_ = -1;
@@ -184,66 +368,204 @@ void Foam::Particle<ParticleType>::trackToFaceByPlanes
 
     if (faces.empty())
     {
-        // inside cell
+        // endPosition is inside the cell
+
+        position_ = endPosition;
+        trackFraction = 1.0;
+        return;
+    }
+
+    // A face has been hit
+
+    scalar lambdaMin = GREAT;
+
+    if (faces.size() == 1)
+    {
+        lambdaMin = lambda(position_, endPosition, faces[0], stepFraction_);
+        facei_ = faces[0];
+    }
+    else
+    {
+        forAll(faces, i)
+        {
+            scalar lam =
+                lambda(position_, endPosition, faces[i], stepFraction_);
+
+            if (lam < lambdaMin)
+            {
+                lambdaMin = lam;
+                facei_ = faces[i];
+            }
+        }
+    }
+
+    // For warped faces the particle can be 'outside' the cell.
+    // This will yield a lambda larger than 1, or smaller than 0.
+
+    // For values < 0, the particle travels away from the cell and we
+    // don't move the particle (except by a small value to move it off
+    // the face), only change cell.
+
+    if (static_cast<ParticleType&>(*this).softImpact())
+    {
+        // Soft-sphere particles can travel outside the domain
+        // but we don't use lambda since this the particle
+        // is going away from face
 
         trackFraction = 1.0;
         position_ = endPosition;
     }
+    else if (lambdaMin <= 0.0)
+    {
+        Pout<< "convex tracking recovery "
+            << origId_ << " "
+            << origProc_ << " "
+            << position_ << " "
+            << endPosition << " "
+            << stepFraction_ << " "
+            << lambdaMin << " "
+            << celli_ << " "
+            << facei_ << " "
+            << endl;
+
+        trackFraction = Cloud<ParticleType>::minValidTrackFraction;
+        position_ += trackFraction*(endPosition - position_);
+    }
     else
     {
-        // hit face
-
-        scalar lambdaMin = GREAT;
-
-        if (faces.size() == 1)
+        if (lambdaMin <= 1.0)
         {
-            lambdaMin = lambda(position_, endPosition, faces[0], stepFraction_);
-            facei_ = faces[0];
+            trackFraction = lambdaMin;
+            position_ += trackFraction*(endPosition - position_);
         }
         else
         {
-            forAll(faces, i)
-            {
-                scalar lam =
-                    lambda(position_, endPosition, faces[i], stepFraction_);
-
-                if (lam < lambdaMin)
-                {
-                    lambdaMin = lam;
-                    facei_ = faces[i];
-                }
-            }
-        }
-
-        // For warped faces the particle can be 'outside' the cell.
-        // This will yield a lambda larger than 1, or smaller than 0
-        // For values < 0, the particle travels away from the cell
-        // and we don't move the particle, only change cell.
-        // For values larger than 1, we move the particle to endPosition only.
-        if (lambdaMin > 0.0)
-        {
-            if (lambdaMin <= 1.0)
-            {
-                trackFraction = lambdaMin;
-                position_ += trackFraction*(endPosition - position_);
-            }
-            else
-            {
-                trackFraction = 1.0;
-                position_ = endPosition;
-            }
-        }
-        else if (static_cast<ParticleType&>(*this).softImpact())
-        {
-            // Soft-sphere particles can travel outside the domain
-            // but we don't use lambda since this the particle
-            // is going away from face
+            // For values larger than 1, we move the particle to endPosition
+            // only.
             trackFraction = 1.0;
             position_ = endPosition;
         }
     }
+
+    faceAction(trackFraction, endPosition, td);
+
+    // If the trackFraction = 0 something went wrong.
+    // Either the particle is flipping back and forth across a face perhaps
+    // due to velocity interpolation errors or it is in a "hole" in the mesh
+    // caused by face warpage.
+    // In both cases resolve the positional ambiguity by moving the particle
+    // slightly towards the cell-centre.
+
+    if (trackFraction < Cloud<ParticleType>::minValidTrackFraction)
+    {
+        Pout<< "convex tracking error "
+            << origId_ << " "
+            << origProc_ << " "
+            << position_ << " "
+            << endPosition << " "
+            << trackFraction << " "
+            << stepFraction_ << " "
+            << lambdaMin << " "
+            << celli_ << " "
+            << facei_ << " "
+            << endl;
+
+        const polyMesh& mesh = cloud_.pMesh();
+
+        const point& cc = mesh.cellCentres()[celli_];
+        position_ +=
+            Cloud<ParticleType>::trackingRescueTolerance*(cc - position_);
+    }
 }
 
+template<class ParticleType>
+template<class TrackData>
+void Foam::Particle<ParticleType>::faceAction
+(
+    scalar& trackFraction,
+    const vector& endPosition,
+    TrackData& td
+)
+{
+    const polyMesh& mesh = cloud_.pMesh();
+
+    if (cloud_.internalFace(facei_))
+    {
+        // Internal face, change cell
+
+        if (celli_ == mesh.faceOwner()[facei_])
+        {
+            celli_ = mesh.faceNeighbour()[facei_];
+        }
+        else if (celli_ == mesh.faceNeighbour()[facei_])
+        {
+            celli_ = mesh.faceOwner()[facei_];
+        }
+        else
+        {
+            FatalErrorIn("Particle::trackToFace(const vector&, TrackData&)")
+            << "addressing failure" << nl
+                << abort(FatalError);
+        }
+    }
+    else
+    {
+        ParticleType& p = static_cast<ParticleType&>(*this);
+
+        // Soft-sphere algorithm ignores the boundary
+        if (p.softImpact())
+        {
+            trackFraction = 1.0;
+            position_ = endPosition;
+        }
+
+        label patchi = patch(facei_);
+        const polyPatch& patch = mesh.boundaryMesh()[patchi];
+
+        if (!p.hitPatch(patch, td, patchi))
+        {
+            if (isA<wedgePolyPatch>(patch))
+            {
+                p.hitWedgePatch
+                (
+                    static_cast<const wedgePolyPatch&>(patch), td
+                );
+            }
+            else if (isA<symmetryPolyPatch>(patch))
+            {
+                p.hitSymmetryPatch
+                (
+                    static_cast<const symmetryPolyPatch&>(patch), td
+                );
+            }
+            else if (isA<cyclicPolyPatch>(patch))
+            {
+                p.hitCyclicPatch
+                (
+                    static_cast<const cyclicPolyPatch&>(patch), td
+                );
+            }
+            else if (isA<processorPolyPatch>(patch))
+            {
+                p.hitProcessorPatch
+                (
+                    static_cast<const processorPolyPatch&>(patch), td
+                );
+            }
+            else if (isA<wallPolyPatch>(patch))
+            {
+                p.hitWallPatch
+                (
+                    static_cast<const wallPolyPatch&>(patch), td
+                );
+            }
+            else
+            {
+                p.hitPatch(patch, td);
+            }
+        }
+    }
+}
 
 template<class ParticleType>
 template<class TrackData>
@@ -391,125 +713,20 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     TrackData& td
 )
 {
-    const polyMesh& mesh = cloud_.polyMesh_;
-
     scalar trackFraction = 0.0;
 
-    if (cloud_.concaveCell()[celli_])
+    if (cloud_.concaveCheck_)
     {
-        trackToFaceByDecomposition(trackFraction, endPosition);
+        if (cloud_.concaveCell()[celli_])
+        {
+            // Use a more careful tracking algorithm if the cell is concave
+            trackToFaceConcave(trackFraction, endPosition, td);
+        }
     }
     else
     {
-        trackToFaceByPlanes(trackFraction, endPosition);
-    }
-
-    if (facei_ > -1)
-    {
-        // A face was hit
-
-        if (cloud_.internalFace(facei_))
-        {
-            // Internal face, change cell
-
-            if (celli_ == mesh.faceOwner()[facei_])
-            {
-                celli_ = mesh.faceNeighbour()[facei_];
-            }
-            else if (celli_ == mesh.faceNeighbour()[facei_])
-            {
-                celli_ = mesh.faceOwner()[facei_];
-            }
-            else
-            {
-                FatalErrorIn("Particle::trackToFace(const vector&, TrackData&)")
-                    << "addressing failure" << nl
-                    << abort(FatalError);
-            }
-        }
-        else
-        {
-            ParticleType& p = static_cast<ParticleType&>(*this);
-
-            // Soft-sphere algorithm ignores the boundary
-            if (p.softImpact())
-            {
-                trackFraction = 1.0;
-                position_ = endPosition;
-            }
-
-            label patchi = patch(facei_);
-            const polyPatch& patch = mesh.boundaryMesh()[patchi];
-
-            if (!p.hitPatch(patch, td, patchi))
-            {
-                if (isA<wedgePolyPatch>(patch))
-                {
-                    p.hitWedgePatch
-                    (
-                        static_cast<const wedgePolyPatch&>(patch), td
-                    );
-                }
-                else if (isA<symmetryPolyPatch>(patch))
-                {
-                    p.hitSymmetryPatch
-                    (
-                        static_cast<const symmetryPolyPatch&>(patch), td
-                    );
-                }
-                else if (isA<cyclicPolyPatch>(patch))
-                {
-                    p.hitCyclicPatch
-                    (
-                        static_cast<const cyclicPolyPatch&>(patch), td
-                    );
-                }
-                else if (isA<processorPolyPatch>(patch))
-                {
-                    p.hitProcessorPatch
-                    (
-                        static_cast<const processorPolyPatch&>(patch), td
-                    );
-                }
-                else if (isA<wallPolyPatch>(patch))
-                {
-                    p.hitWallPatch
-                    (
-                        static_cast<const wallPolyPatch&>(patch), td
-                    );
-                }
-                else
-                {
-                    p.hitPatch(patch, td);
-                }
-            }
-        }
-    }
-
-    // If the trackFraction = 0 something went wrong.
-    // Either the particle is flipping back and forth across a face perhaps
-    // due to velocity interpolation errors or it is in a "hole" in the mesh
-    // caused by face warpage.
-    // In both cases resolve the positional ambiguity by moving the particle
-    // slightly towards the cell-centre.
-
-    if (trackFraction < Cloud<ParticleType>::minValidTrackFraction)
-    {
-        // Pout<< "tracking error "
-        //     << origId_ << " "
-        //     << origProc_ << " "
-        //     << position_ << " "
-        //     << endPosition << " "
-        //     << trackFraction << " "
-        //     << stepFraction_ << " "
-        //     << celli_ << " "
-        //     << facei_ << " "
-        //     << cloud_.concaveCell()[celli_]
-        //     << endl;
-
-        const point& cc = mesh.cellCentres()[celli_];
-        position_ +=
-            Cloud<ParticleType>::trackingRescueTolerance*(cc - position_);
+        // Use the original tracking algorithm if the cell is convex
+        trackToFaceConvex(trackFraction, endPosition, td);
     }
 
     return trackFraction;
diff --git a/src/lagrangian/basic/Particle/Particle.H b/src/lagrangian/basic/Particle/Particle.H
index cf44cc1198..0aea781b55 100644
--- a/src/lagrangian/basic/Particle/Particle.H
+++ b/src/lagrangian/basic/Particle/Particle.H
@@ -157,17 +157,17 @@ protected:
             const label facei
         ) const;
 
-        //- Find the faces between position and cell centre
+        //- Find the faces between endPosition and cell centre
         void findFaces
         (
-            const vector& position,
+            const vector& endPosition,
             DynamicList<label>& faceList
         ) const;
 
-        //- Find the faces between position and cell centre
+        //- Find the faces between endPosition and cell centre
         void findFaces
         (
-            const vector& position,
+            const vector& endPosition,
             const label celli,
             const scalar stepFraction,
             DynamicList<label>& faceList
@@ -175,20 +175,32 @@ protected:
 
         //- Track to cell face using face decomposition, used for
         //  concave cells.
-        void trackToFaceByDecomposition
+        template<class TrackData>
+        void trackToFaceConcave
         (
             scalar& trackFraction,
-            const vector& endPosition
+            const vector& endPosition,
+            TrackData& td
         );
 
         //- Track to cell face using the infinite planes of the faces.
         //  The cell *must* be convex.
-        void trackToFaceByPlanes
+        template<class TrackData>
+        void trackToFaceConvex
         (
             scalar& trackFraction,
-            const vector& endPosition
+            const vector& endPosition,
+            TrackData& td
         );
 
+        //- Change cell or interact with a patch as required
+        template<class TrackData>
+        void faceAction
+        (
+            scalar& trackFraction,
+            const vector& endPosition,
+            TrackData& td
+        );
 
 
     // Patch interactions

From 2d8e739837a29c04fde143076c01907e84e505e5 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Wed, 7 Oct 2009 17:08:26 +0100
Subject: [PATCH 04/16] Adding incompressible solver where the value of density
 is read from the transportDictionary.

---
 .../Make/files                                |  3 +
 .../Make/options                              | 25 ++++++
 .../createFields.H                            | 85 +++++++++++++++++++
 ...mpressibleInteractingKinematicParcelFoam.C | 78 +++++++++++++++++
 4 files changed, 191 insertions(+)
 create mode 100644 applications/solvers/lagrangian/incompressibleInteractingKinematicParcelFoam/Make/files
 create mode 100644 applications/solvers/lagrangian/incompressibleInteractingKinematicParcelFoam/Make/options
 create mode 100644 applications/solvers/lagrangian/incompressibleInteractingKinematicParcelFoam/createFields.H
 create mode 100644 applications/solvers/lagrangian/incompressibleInteractingKinematicParcelFoam/incompressibleInteractingKinematicParcelFoam.C

diff --git a/applications/solvers/lagrangian/incompressibleInteractingKinematicParcelFoam/Make/files b/applications/solvers/lagrangian/incompressibleInteractingKinematicParcelFoam/Make/files
new file mode 100644
index 0000000000..918868e1ce
--- /dev/null
+++ b/applications/solvers/lagrangian/incompressibleInteractingKinematicParcelFoam/Make/files
@@ -0,0 +1,3 @@
+incompressibleInteractingKinematicParcelFoam.C
+
+EXE = $(FOAM_APPBIN)/incompressibleInteractingKinematicParcelFoam
diff --git a/applications/solvers/lagrangian/incompressibleInteractingKinematicParcelFoam/Make/options b/applications/solvers/lagrangian/incompressibleInteractingKinematicParcelFoam/Make/options
new file mode 100644
index 0000000000..845dbb7b45
--- /dev/null
+++ b/applications/solvers/lagrangian/incompressibleInteractingKinematicParcelFoam/Make/options
@@ -0,0 +1,25 @@
+EXE_INC = \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+
+EXE_LIBS = \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -lthermophysicalFunctions \
+    -lbasicThermophysicalModels \
+    -lspecie \
+    -lradiation \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools
diff --git a/applications/solvers/lagrangian/incompressibleInteractingKinematicParcelFoam/createFields.H b/applications/solvers/lagrangian/incompressibleInteractingKinematicParcelFoam/createFields.H
new file mode 100644
index 0000000000..f816d38607
--- /dev/null
+++ b/applications/solvers/lagrangian/incompressibleInteractingKinematicParcelFoam/createFields.H
@@ -0,0 +1,85 @@
+    Info<< "\nReading transportProperties\n" << endl;
+
+    IOdictionary transportProperties
+    (
+        IOobject
+        (
+            "transportProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    dimensionedScalar rhoValue
+    (
+        transportProperties.lookup("rho")
+    );
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        rhoValue
+    );
+
+    Info<< "\nReading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "createPhi.H"
+
+    Info<< "Creating turbulence model\n" << endl;
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
+    );
+
+    volScalarField mu
+    (
+        IOobject
+        (
+            "mu",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        laminarTransport.nu()*rhoValue
+    );
+
+    word interactingKinematicCloudName("interactingKinematicCloud");
+    args.optionReadIfPresent("cloudName", interactingKinematicCloudName);
+
+    Info<< "Constructing interactingKinematicCloud "
+        << interactingKinematicCloudName << endl;
+
+    basicInteractingKinematicCloud interactingKinematicCloud
+    (
+        interactingKinematicCloudName,
+        rho,
+        U,
+        mu,
+        g
+    );
diff --git a/applications/solvers/lagrangian/incompressibleInteractingKinematicParcelFoam/incompressibleInteractingKinematicParcelFoam.C b/applications/solvers/lagrangian/incompressibleInteractingKinematicParcelFoam/incompressibleInteractingKinematicParcelFoam.C
new file mode 100644
index 0000000000..dd70f50d49
--- /dev/null
+++ b/applications/solvers/lagrangian/incompressibleInteractingKinematicParcelFoam/incompressibleInteractingKinematicParcelFoam.C
@@ -0,0 +1,78 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    interactingKinematicParcelFoam
+
+Description
+    Transient solver for the passive transport of a single interactingKinematic
+    particle could, where the particles interact.
+
+    Uses a pre-calculated velocity field to evolve the cloud.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulenceModel.H"
+#include "basicInteractingKinematicCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::validOptions.insert("cloudName", "cloud name");
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readGravitationalAcceleration.H"
+    #include "createFields.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        Info<< "Evolving " << interactingKinematicCloud.name() << endl;
+        interactingKinematicCloud.evolve();
+        interactingKinematicCloud.info();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

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 05/16] 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
     {

From 37d84af4530aba1819b15f4974e473e489e8fdbf Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Thu, 8 Oct 2009 18:53:03 +0100
Subject: [PATCH 06/16] Reducing compilation time by removing items from
 Make/files, temporary

---
 src/lagrangian/intermediate/Make/files | 38 +++++++++++++-------------
 1 file changed, 19 insertions(+), 19 deletions(-)

diff --git a/src/lagrangian/intermediate/Make/files b/src/lagrangian/intermediate/Make/files
index 7c95a699aa..4f4e4bd6a1 100644
--- a/src/lagrangian/intermediate/Make/files
+++ b/src/lagrangian/intermediate/Make/files
@@ -8,14 +8,14 @@ DERIVEDCLOUDS=$(CLOUDS)/derived
 
 
 /* Parcels */
-$(BASEPARCELS)/reactingParcel/reactingParcel.C
+// $(BASEPARCELS)/reactingParcel/reactingParcel.C
 
 
 /* Cloud base classes */
 $(BASECLOUDS)/kinematicCloud/kinematicCloud.C
-$(BASECLOUDS)/thermoCloud/thermoCloud.C
-$(BASECLOUDS)/reactingCloud/reactingCloud.C
-$(BASECLOUDS)/reactingMultiphaseCloud/reactingMultiphaseCloud.C
+// $(BASECLOUDS)/thermoCloud/thermoCloud.C
+// $(BASECLOUDS)/reactingCloud/reactingCloud.C
+// $(BASECLOUDS)/reactingMultiphaseCloud/reactingMultiphaseCloud.C
 
 
 /* kinematic parcel sub-models */
@@ -33,27 +33,27 @@ $(INTERACTINGKINEMATICPARCEL)/makeBasicInteractingKinematicParcelSubmodels.C
 
 
 /* thermo parcel sub-models */
-THERMOPARCEL=$(DERIVEDPARCELS)/basicThermoParcel
-$(THERMOPARCEL)/basicThermoParcel.C
-$(THERMOPARCEL)/defineBasicThermoParcel.C
-$(THERMOPARCEL)/makeBasicThermoParcelSubmodels.C
+// THERMOPARCEL=$(DERIVEDPARCELS)/basicThermoParcel
+// $(THERMOPARCEL)/basicThermoParcel.C
+// $(THERMOPARCEL)/defineBasicThermoParcel.C
+// $(THERMOPARCEL)/makeBasicThermoParcelSubmodels.C
 
 
 /* reacting parcel sub-models */
-REACTINGPARCEL=$(DERIVEDPARCELS)/BasicReactingParcel
-$(REACTINGPARCEL)/defineBasicReactingParcel.C
-$(REACTINGPARCEL)/makeBasicReactingParcelSubmodels.C
+// REACTINGPARCEL=$(DERIVEDPARCELS)/BasicReactingParcel
+// $(REACTINGPARCEL)/defineBasicReactingParcel.C
+// $(REACTINGPARCEL)/makeBasicReactingParcelSubmodels.C
 
 
 /* reacting multiphase parcel sub-models */
-REACTINGMPPARCEL=$(DERIVEDPARCELS)/BasicReactingMultiphaseParcel
-$(REACTINGMPPARCEL)/defineBasicReactingMultiphaseParcel.C
-$(REACTINGMPPARCEL)/makeBasicReactingMultiphaseParcelSubmodels.C
+// REACTINGMPPARCEL=$(DERIVEDPARCELS)/BasicReactingMultiphaseParcel
+// $(REACTINGMPPARCEL)/defineBasicReactingMultiphaseParcel.C
+// $(REACTINGMPPARCEL)/makeBasicReactingMultiphaseParcelSubmodels.C
 
 
 /* bolt-on models */
-submodels/addOns/radiation/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C
-submodels/addOns/radiation/scatter/cloudScatter/cloudScatter.C
+// submodels/addOns/radiation/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C
+// submodels/addOns/radiation/scatter/cloudScatter/cloudScatter.C
 
 
 /* data entries */
@@ -71,9 +71,9 @@ particleForces/particleForces.C
 
 
 /* phase properties */
-phaseProperties/phaseProperties/phaseProperties.C
-phaseProperties/phaseProperties/phasePropertiesIO.C
-phaseProperties/phasePropertiesList/phasePropertiesList.C
+// phaseProperties/phaseProperties/phaseProperties.C
+// phaseProperties/phaseProperties/phasePropertiesIO.C
+// phaseProperties/phasePropertiesList/phasePropertiesList.C
 
 
 LIB = $(FOAM_LIBBIN)/liblagrangianIntermediate

From 3ef3638c44ffcad4a9a75528a61df20dd733a812 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Fri, 9 Oct 2009 10:39:06 +0100
Subject: [PATCH 07/16] Fixing bug in faceIntersection where the sign of the
 distance was being discarded on return from the triPointRef intersecion test.
  Was only appearing when using FULL_RAY tests.

---
 src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C b/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C
index c5f9a80250..cf7dfc9400 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C
+++ b/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C
@@ -159,9 +159,9 @@ Foam::pointHit Foam::face::intersection
 
         if (curHit.hit())
         {
-            if (Foam::mag(curHit.distance()) < nearestHitDist)
+            if (Foam::mag(curHit.distance()) < Foam::mag(nearestHitDist))
             {
-                nearestHitDist = Foam::mag(curHit.distance());
+                nearestHitDist = curHit.distance();
                 nearest.setHit();
                 nearest.setPoint(curHit.hitPoint());
             }

From 5488acac272911f60d2085bd9598283171328d5d Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Fri, 9 Oct 2009 11:01:26 +0100
Subject: [PATCH 08/16] Small changes to output during best guess
 determination.

Making intersection test tolerance SMALL.
---
 src/lagrangian/basic/Cloud/Cloud.C       |  2 +-
 src/lagrangian/basic/Particle/Particle.C | 29 +-----------------------
 2 files changed, 2 insertions(+), 29 deletions(-)

diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
index 81bf07b378..c8f59f2904 100644
--- a/src/lagrangian/basic/Cloud/Cloud.C
+++ b/src/lagrangian/basic/Cloud/Cloud.C
@@ -42,7 +42,7 @@ template<class ParticleType>
 const Foam::scalar Foam::Cloud<ParticleType>::trackingRescueTolerance = 1e-6;
 
 template<class ParticleType>
-const Foam::scalar Foam::Cloud<ParticleType>::intersectionTolerance = 1e-3;
+const Foam::scalar Foam::Cloud<ParticleType>::intersectionTolerance = SMALL;
 
 template<class ParticleType>
 const Foam::scalar Foam::Cloud<ParticleType>::planarCosAngle = (1 - 1e-6);
diff --git a/src/lagrangian/basic/Particle/Particle.C b/src/lagrangian/basic/Particle/Particle.C
index cbc4a9ccc0..437cabd7b8 100644
--- a/src/lagrangian/basic/Particle/Particle.C
+++ b/src/lagrangian/basic/Particle/Particle.C
@@ -281,6 +281,7 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
 
         Pout<< "Need a best guess " << nl
             << position_ << nl
+            << stepFraction_ << nl
             << position_ + delta << nl
             << position_ - delta
             << endl;
@@ -315,34 +316,6 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
                 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;

From 20ab51ca50a3be1d932dd73b8918e6dbb5a91e17 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Fri, 9 Oct 2009 16:36:11 +0100
Subject: [PATCH 09/16] Trying various methods to decide what to do.  Best
 guess of which cell the particle is in is off, going to try a different way.

---
 src/lagrangian/basic/Cloud/Cloud.C       |   2 +-
 src/lagrangian/basic/Particle/Particle.C | 391 +++++++++++++++--------
 src/lagrangian/basic/Particle/Particle.H |   9 +
 3 files changed, 263 insertions(+), 139 deletions(-)

diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
index c8f59f2904..4ee3a1e38d 100644
--- a/src/lagrangian/basic/Cloud/Cloud.C
+++ b/src/lagrangian/basic/Cloud/Cloud.C
@@ -42,7 +42,7 @@ template<class ParticleType>
 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 = 0;
 
 template<class ParticleType>
 const Foam::scalar Foam::Cloud<ParticleType>::planarCosAngle = (1 - 1e-6);
diff --git a/src/lagrangian/basic/Particle/Particle.C b/src/lagrangian/basic/Particle/Particle.C
index 437cabd7b8..dfc9f6c5c1 100644
--- a/src/lagrangian/basic/Particle/Particle.C
+++ b/src/lagrangian/basic/Particle/Particle.C
@@ -89,53 +89,30 @@ void Foam::Particle<ParticleType>::findFaces
 
 
 template<class ParticleType>
-template<class TrackData>
-void Foam::Particle<ParticleType>::trackToFaceConcave
+bool Foam::Particle<ParticleType>::insideCellExact
 (
-    scalar& trackFraction,
-    const vector& endPosition,
-    TrackData& td
-)
+    const vector& testPt,
+    const label celli,
+    bool includeOnFace
+) const
 {
-    facei_ = -1;
-
     const polyMesh& mesh = cloud_.pMesh();
-    const labelList& faces = mesh.cells()[celli_];
-    const vector& C = mesh.cellCentres()[celli_];
-
-    DynamicList<label>& potentiallyCrossedFaces = cloud_.labels_;
-    findFaces(endPosition, potentiallyCrossedFaces);
-
-    // Check all possible face crossings to see if they are actually
-    // crossed, determining if endPosition is outside the current cell.
-
-    if (potentiallyCrossedFaces.empty())
-    {
-        // endPosition is inside the cell
-
-        position_ = endPosition;
-        trackFraction = 1.0;
-        return;
-    }
-
-    Pout<< nl << origProc_ << " "
-        << origId_ << " "
-        << position_ << " "
-        << endPosition << " "
-        << celli_ << nl
-        << potentiallyCrossedFaces
-        << endl;
+    const labelList& faces = mesh.cells()[celli];
+    const vector& C = mesh.cellCentres()[celli];
 
     label nFaceCrossings = 0;
 
-    forAll (potentiallyCrossedFaces, pCFI)
+    // The vector from the cell centre to the end point
+    vector delta = testPt - C;
+
+    forAll (faces, i)
     {
-        label facei = potentiallyCrossedFaces[pCFI];
+        label facei = faces[i];
 
         pointHit inter = mesh.faces()[facei].intersection
         (
             C,
-            endPosition - C,
+            delta,
             mesh.faceCentres()[facei],
             mesh.points(),
             intersection::HALF_RAY,
@@ -144,18 +121,25 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
 
         if (inter.hit())
         {
-            if (inter.distance() < 1.0)
+            // Pout<< "insideCellExact " << facei << " "
+            //     << inter.distance() << endl;
+
+            if (includeOnFace)
             {
-                // This face was actually crossed.
-
-                nFaceCrossings++;
-
-                if (inter.distance() < 0.0)
+                if (inter.distance() <= 1.0)
                 {
-                    FatalErrorIn("Particle::trackToFaceConcave")
-                        << "Negative lambda " << inter.distance()
-                        << "in potentiallyCrossedFaces test"
-                        << nl << abort(FatalError);
+                    // This face was actually crossed.
+
+                    nFaceCrossings++;
+                }
+            }
+            else
+            {
+                if (inter.distance() < 1.0)
+                {
+                    // This face was actually crossed.
+
+                    nFaceCrossings++;
                 }
             }
         }
@@ -169,9 +153,42 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
             << endl;
     }
 
-    Pout<< nFaceCrossings << " " << (nFaceCrossings % 2) << endl;
-
     if (nFaceCrossings % 2 == 0)
+    {
+        // Even number of face crossings, so the testPt must be in the
+        // cell.
+
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class ParticleType>
+template<class TrackData>
+void Foam::Particle<ParticleType>::trackToFaceConcave
+(
+    scalar& trackFraction,
+    const vector& endPosition,
+    TrackData& td
+)
+{
+    facei_ = -1;
+
+    const polyMesh& mesh = cloud_.pMesh();
+    const labelList& faces = mesh.cells()[celli_];
+    // const vector& C = mesh.cellCentres()[celli_];
+
+    // Check all possible face crossings to see if they are actually
+    // crossed, determining if endPosition is outside the current
+    // cell.  This allows situations where the cell is outside the
+    // cell to start with and enters the cell at the end of the track
+    // to be identified.
+
+    // Pout<< nl << "Outside test:" << endl;
+
+    if (insideCellExact(endPosition, celli_, false))
     {
         // Even number of face crossings, so the particle must end up
         // still in the cell.
@@ -181,11 +198,20 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
         return;
     }
 
+    Pout<< nl << origProc_ << " "
+        << origId_ << " "
+        << position_ << " "
+        << endPosition << " "
+        << stepFraction_ << " "
+        << celli_
+        << endl;
+
     // The particle *must* have left the cell.
 
-    // a) It may have crossed a face not yet identified by findFaces using the
-    //    cell centre to endPosition line, so the potentially crossed faces of
-    //    the position to endPosition line must be assessed.
+    // a) It may have crossed a face not yet identified by testing
+    //    faces using the cell centre to endPosition line, so the
+    //    potentially crossed faces of the position to endPosition
+    //    line must be assessed.
 
     // b) It may have been outside the cell in the first place, and, despite
     //    trying to pick up more faces using a) the correct face to be crossed
@@ -197,96 +223,144 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
     // as nothing can be assumed about the order of crossing the
     // planes of faces.
 
-    forAll(faces, i)
+    const vector deltaPosition = endPosition - position_;
+    vector deltaTrack =
+        mag(mesh.cellCentres()[celli_] - position_)
+       *deltaPosition/(mag(deltaPosition) + VSMALL);
+
+    // Pout<< "Inside test:" << endl;
+
+    if (insideCellExact(position_, celli_, false))
     {
-        label facei = faces[i];
+        Pout<< "The particle starts inside the cell and ends up outside of it"
+            << nl << position_ << " " << position_ + deltaTrack
+            << endl;
 
-        scalar lam = lambda(position_, endPosition, facei);
+        // The particle started inside the cell and finished outside
+        // of it, find which face to cross
 
-        if ((lam > 0) && (lam < 1.0))
+        scalar tmpLambda = GREAT;
+        scalar correctLambda = GREAT;
+
+        forAll(faces, i)
         {
-            if (findIndex(potentiallyCrossedFaces, facei) == -1)
-            {
-                // A new potential face crossing, previously missed,
-                // has been found
+            label facei = faces[i];
 
-                Pout<< "New face crossing " << facei
-                    << " of concave cell " << celli_ << " found"
-                    << endl;
+            // Use exact intersection.
 
-                potentiallyCrossedFaces.append(facei);
-            }
-        }
-    }
+            // TODO: A correction is required for moving meshes to
+            // calculate the correct lambda value.
 
-    Pout<< potentiallyCrossedFaces << endl;
-
-    vector deltaPosition = endPosition - position_;
-
-    scalar tmpLambda = GREAT;
-    scalar correctLambda = GREAT;
-
-    forAll(potentiallyCrossedFaces, pCFI)
-    {
-        label facei = potentiallyCrossedFaces[pCFI];
-
-        // Use exact intersection.
-
-        // 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::HALF_RAY,
-            Cloud<ParticleType>::intersectionTolerance
-        );
-
-        if (inter.hit())
-        {
-            tmpLambda = inter.distance();
-
-            if
+            pointHit inter = mesh.faces()[facei].intersection
             (
-                tmpLambda < 1.0
-             && tmpLambda < correctLambda
-            )
-            {
-                // This face is crossed before any other that has
-                // been found so far
+                position_,
+                deltaPosition,
+                mesh.faceCentres()[facei],
+                mesh.points(),
+                intersection::HALF_RAY,
+                Cloud<ParticleType>::intersectionTolerance
+            );
 
-                correctLambda = tmpLambda;
-                facei_ = facei;
+            if (inter.hit())
+            {
+                tmpLambda = inter.distance();
+
+                Pout<< facei << " " << tmpLambda << endl;
+
+                if
+                (
+                    tmpLambda <= 1.0
+                 && tmpLambda < correctLambda
+                )
+                {
+                    // This face is crossed before any other that has
+                    // been found so far
+
+                    correctLambda = tmpLambda;
+                    facei_ = facei;
+                }
             }
         }
-    }
 
-    if (facei_ > -1)
-    {
-        trackFraction = correctLambda;
-        position_ += trackFraction*(endPosition - position_);
+        if (facei_ > -1)
+        {
+            if (!cloud_.internalFace(facei_))
+            {
+                label patchi = patch(facei_);
+                const polyPatch& patch = mesh.boundaryMesh()[patchi];
+
+                if (isA<wallPolyPatch>(patch))
+                {
+                    if ((mesh.faceAreas()[facei_] & deltaPosition) <= 0)
+                    {
+                        // The particle has hit a wall face but it is
+                        // heading in the wrong direction with respect to
+                        // the face normal
+
+                        // Do not trigger a face hit and move the position
+                        // towards the cell centre
+
+                        Pout<< "Hit a wall face heading the wrong way"
+                            << endl;
+
+                        const point& cc = mesh.cellCentres()[celli_];
+                        position_ += 1e-2*(cc - position_);
+
+                        facei_ = -1;
+                    }
+                }
+            }
+            else if (correctLambda < Cloud<ParticleType>::minValidTrackFraction)
+            {
+                // The particle is not far enough away from the face
+                // to decide if it is valid crossing. Let it move a
+                // little without crossing the face to resolve the
+                // ambiguity.
+
+                Pout<< "Ambiguous face crossing, correcting towards cell "
+                    << "centre and not crossing face" << endl;
+
+                const point& cc = mesh.cellCentres()[celli_];
+                position_ +=
+                    Cloud<ParticleType>::trackingRescueTolerance
+                   *(cc - position_);
+
+                // Pout<< "Ambiguous face crossing. " << endl;
+                // correctLambda += Cloud<ParticleType>::minValidTrackFraction;
+
+                facei_ = -1;
+            }
+
+            trackFraction = correctLambda;
+            position_ += trackFraction*(endPosition - position_);
+        }
+        else
+        {
+            Pout<< "Particle " << origProc_ << " " << origId_
+                << " started inside cell " << celli_ << " and finished outside"
+                << " of it, but did not find a face to cross"
+                << endl;
+
+            const point& cc = mesh.cellCentres()[celli_];
+            position_ += 1e-2*(cc - position_);
+        }
     }
     else
     {
+        Pout<< "The particle started outside of the cell" << endl;
+
         // 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.
 
-        vector delta =
-            mag(mesh.cellCentres()[celli_] - position_)
-           *deltaPosition/mag(deltaPosition);
-
         Pout<< "Need a best guess " << nl
             << position_ << nl
             << stepFraction_ << nl
-            << position_ + delta << nl
-            << position_ - delta
+            << position_ + deltaTrack << nl
+            << position_ - deltaTrack
             << endl;
 
-        // ---- Best guess by projection --
+        // Best guess by projection
 
         DynamicList<scalar> tmpLambdas;
         DynamicList<label> tmpLambdaFaceIs;
@@ -320,38 +394,73 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
 
         Pout<< tmpLambdaFaceIs << " " << tmpLambdas << endl;
 
-        // The closest lambda value (less than zero) or (greater than or
-        // equal to one), and corresponding face.
         scalar closestLambda = GREAT;
-        scalar closestLambdaDifference = GREAT;
-        label closestLambdaDifferenceFaceI = -1;
+        label closestLambdaFaceI = -1;
 
+        // Is there a face crossing that has been missed?
         forAll(tmpLambdaFaceIs, i)
         {
             label facei = tmpLambdaFaceIs[i];
             scalar tmpLambda = tmpLambdas[i];
 
-            if (tmpLambda < 0)
+            if (!cloud_.internalFace(facei))
             {
-                if (mag(tmpLambda) < closestLambdaDifference)
-                {
-                    closestLambdaDifference = mag(tmpLambda);
-                    closestLambda = tmpLambda;
-                    closestLambdaDifferenceFaceI = facei;
-                }
+                // Boundary faces trump internal faces
+                closestLambda = tmpLambda;
+                closestLambdaFaceI = facei;
+
+                break;
             }
-            else if (tmpLambda >= 1.0)
+            else if
+            (
+                tmpLambda >= 0.0
+             && tmpLambda <= 1.0
+             && tmpLambda < closestLambda
+            )
             {
-                if ((tmpLambda - 1.0) < closestLambdaDifference)
+                closestLambda = tmpLambda;
+                closestLambdaFaceI = facei;
+            }
+        }
+
+        // Find the closest lambda by projection.
+        if (closestLambdaFaceI == -1)
+        {
+            Pout<< "Test by projection" << endl;
+
+            // The closest lambda value (less than zero) or (greater than or
+            // equal to one), and corresponding face.
+
+            // Limit how far away from the track a projected face can be
+            scalar closestLambdaDifference = 100.0;
+
+            forAll(tmpLambdaFaceIs, i)
+            {
+                label facei = tmpLambdaFaceIs[i];
+                scalar tmpLambda = tmpLambdas[i];
+
+                if (tmpLambda < 0)
                 {
-                    closestLambdaDifference = tmpLambda - 1.0;
-                    closestLambda = tmpLambda;
-                    closestLambdaDifferenceFaceI = facei;
+                    if (mag(tmpLambda) < closestLambdaDifference)
+                    {
+                        closestLambdaDifference = mag(tmpLambda);
+                        closestLambda = tmpLambda;
+                        closestLambdaFaceI = facei;
+                    }
+                }
+                else if (tmpLambda >= 1.0)
+                {
+                    if ((tmpLambda - 1.0) < closestLambdaDifference)
+                    {
+                        closestLambdaDifference = tmpLambda - 1.0;
+                        closestLambda = tmpLambda;
+                        closestLambdaFaceI = facei;
+                    }
                 }
             }
         }
 
-        if (closestLambdaDifferenceFaceI > -1)
+        if (closestLambdaFaceI > -1)
         {
             // A face has been found to cross.  Not moving the
             // particle at all, but the option is there to decide do
@@ -359,7 +468,10 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
             // if the particle is heading towards, or away from the
             // face.
 
-            facei_ = closestLambdaDifferenceFaceI;
+            facei_ = closestLambdaFaceI;
+
+            trackFraction = Cloud<ParticleType>::minValidTrackFraction;
+            position_ += trackFraction*(endPosition - position_);
         }
         else
         {
@@ -373,7 +485,6 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
                 << celli_ << " "
                 << endl;
 
-            // ---- Best guess by nearest ----
             scalar minDistance = GREAT;
 
             forAll(faces, i)
@@ -396,13 +507,17 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
                 }
             }
 
-            // ---- End of best guess by nearest ----
+            trackFraction = Cloud<ParticleType>::minValidTrackFraction;
+            position_ += trackFraction*(endPosition - position_);
         }
     }
 
     Pout<< facei_ << " " << celli_ << endl;
 
-    faceAction(trackFraction, endPosition, td);
+    if (facei_ > -1)
+    {
+        faceAction(trackFraction, endPosition, td);
+    }
 
     Pout<< facei_ << " " << celli_ << endl;
 }
diff --git a/src/lagrangian/basic/Particle/Particle.H b/src/lagrangian/basic/Particle/Particle.H
index 0aea781b55..42cc5483ad 100644
--- a/src/lagrangian/basic/Particle/Particle.H
+++ b/src/lagrangian/basic/Particle/Particle.H
@@ -173,6 +173,15 @@ protected:
             DynamicList<label>& faceList
         ) const;
 
+        //- Determine if the testPt is inside the cell by an exact
+        //  intersection test
+        bool insideCellExact
+        (
+            const vector& testPt,
+            const label celli,
+            bool includeOnFace
+        ) const;
+
         //- Track to cell face using face decomposition, used for
         //  concave cells.
         template<class TrackData>

From 9f6335164947373965a56b764f79eabb5520f20d Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Fri, 9 Oct 2009 17:45:04 +0100
Subject: [PATCH 10/16] New method for finding which cell a particle should be
 in if it is not in the cell it thinks that it is when it starts and the
 endPosition is also outside - an exact search off all point, edge and face
 connected cells.  If it still isn't found, a recovering move to cell centre
 is applied.

---
 src/lagrangian/basic/Cloud/Cloud.C       |   2 +-
 src/lagrangian/basic/Particle/Particle.C | 180 +++++------------------
 src/lagrangian/basic/Particle/Particle.H |   2 +-
 3 files changed, 38 insertions(+), 146 deletions(-)

diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
index 4ee3a1e38d..e0679c2d45 100644
--- a/src/lagrangian/basic/Cloud/Cloud.C
+++ b/src/lagrangian/basic/Cloud/Cloud.C
@@ -36,7 +36,7 @@ License
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 template<class ParticleType>
-const Foam::scalar Foam::Cloud<ParticleType>::minValidTrackFraction = 1e-10;
+const Foam::scalar Foam::Cloud<ParticleType>::minValidTrackFraction = 1e-6;
 
 template<class ParticleType>
 const Foam::scalar Foam::Cloud<ParticleType>::trackingRescueTolerance = 1e-6;
diff --git a/src/lagrangian/basic/Particle/Particle.C b/src/lagrangian/basic/Particle/Particle.C
index dfc9f6c5c1..c107a38e13 100644
--- a/src/lagrangian/basic/Particle/Particle.C
+++ b/src/lagrangian/basic/Particle/Particle.C
@@ -93,7 +93,7 @@ bool Foam::Particle<ParticleType>::insideCellExact
 (
     const vector& testPt,
     const label celli,
-    bool includeOnFace
+    bool beingOnAFaceMeansOutside
 ) const
 {
     const polyMesh& mesh = cloud_.pMesh();
@@ -121,10 +121,11 @@ bool Foam::Particle<ParticleType>::insideCellExact
 
         if (inter.hit())
         {
-            // Pout<< "insideCellExact " << facei << " "
+            // Pout<< "insideCellExact cell " << celli
+            //     << " face " << facei << " "
             //     << inter.distance() << endl;
 
-            if (includeOnFace)
+            if (beingOnAFaceMeansOutside)
             {
                 if (inter.distance() <= 1.0)
                 {
@@ -178,7 +179,6 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
 
     const polyMesh& mesh = cloud_.pMesh();
     const labelList& faces = mesh.cells()[celli_];
-    // const vector& C = mesh.cellCentres()[celli_];
 
     // Check all possible face crossings to see if they are actually
     // crossed, determining if endPosition is outside the current
@@ -349,166 +349,58 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
     {
         Pout<< "The particle started outside of the cell" << endl;
 
-        // 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.
+        // Find which cell the particle should be in.
 
-        Pout<< "Need a best guess " << nl
-            << position_ << nl
-            << stepFraction_ << nl
-            << position_ + deltaTrack << nl
-            << position_ - deltaTrack
-            << endl;
+        const labelList& cPts = mesh.cellPoints(celli_);
 
-        // Best guess by projection
+        DynamicList<label> checkedCells;
 
-        DynamicList<scalar> tmpLambdas;
-        DynamicList<label> tmpLambdaFaceIs;
+        bool found = false;
 
-        forAll(faces, i)
+        forAll(cPts, cPtI)
         {
-            label facei = faces[i];
+            label ptI = cPts[cPtI];
 
-            // Use exact FULL_RAY intersection to allow negative
-            // values
+            const labelList& pCs = mesh.pointCells(ptI);
 
-            // 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())
+            forAll(pCs, pCI)
             {
-                tmpLambdas.append(inter.distance());
-                tmpLambdaFaceIs.append(facei);
+                label cellI = pCs[pCI];
+
+                if (findIndex(checkedCells, cellI) == -1)
+                {
+                    checkedCells.append(cellI);
+
+                    if (insideCellExact(position_, cellI, false))
+                    {
+                        found = true;
+
+                        celli_ = cellI;
+
+                        break;
+                    }
+                }
             }
-        }
 
-        Pout<< tmpLambdaFaceIs << " " << tmpLambdas << endl;
-
-        scalar closestLambda = GREAT;
-        label closestLambdaFaceI = -1;
-
-        // Is there a face crossing that has been missed?
-        forAll(tmpLambdaFaceIs, i)
-        {
-            label facei = tmpLambdaFaceIs[i];
-            scalar tmpLambda = tmpLambdas[i];
-
-            if (!cloud_.internalFace(facei))
+            if (found)
             {
-                // Boundary faces trump internal faces
-                closestLambda = tmpLambda;
-                closestLambdaFaceI = facei;
-
                 break;
             }
-            else if
-            (
-                tmpLambda >= 0.0
-             && tmpLambda <= 1.0
-             && tmpLambda < closestLambda
-            )
-            {
-                closestLambda = tmpLambda;
-                closestLambdaFaceI = facei;
-            }
         }
 
-        // Find the closest lambda by projection.
-        if (closestLambdaFaceI == -1)
+        if (!found)
         {
-            Pout<< "Test by projection" << endl;
+            Pout<< "Didn't find a new cell after searching "
+                << checkedCells << endl;
 
-            // The closest lambda value (less than zero) or (greater than or
-            // equal to one), and corresponding face.
-
-            // Limit how far away from the track a projected face can be
-            scalar closestLambdaDifference = 100.0;
-
-            forAll(tmpLambdaFaceIs, i)
-            {
-                label facei = tmpLambdaFaceIs[i];
-                scalar tmpLambda = tmpLambdas[i];
-
-                if (tmpLambda < 0)
-                {
-                    if (mag(tmpLambda) < closestLambdaDifference)
-                    {
-                        closestLambdaDifference = mag(tmpLambda);
-                        closestLambda = tmpLambda;
-                        closestLambdaFaceI = facei;
-                    }
-                }
-                else if (tmpLambda >= 1.0)
-                {
-                    if ((tmpLambda - 1.0) < closestLambdaDifference)
-                    {
-                        closestLambdaDifference = tmpLambda - 1.0;
-                        closestLambda = tmpLambda;
-                        closestLambdaFaceI = facei;
-                    }
-                }
-            }
-        }
-
-        if (closestLambdaFaceI > -1)
-        {
-            // A face has been found to cross.  Not moving the
-            // particle at all, but the option is there to decide do
-            // it, as the sign of closestLambda is available to tell
-            // if the particle is heading towards, or away from the
-            // face.
-
-            facei_ = closestLambdaFaceI;
-
-            trackFraction = Cloud<ParticleType>::minValidTrackFraction;
-            position_ += trackFraction*(endPosition - position_);
+            const point& cc = mesh.cellCentres()[celli_];
+            position_ += 1e-2*(cc - position_);
         }
         else
         {
-            Pout<< "Tracking in concave cell - didn't find a "
-                << "best guess face" << nl
-                << origId_ << " "
-                << origProc_ << " "
-                << position_ << " "
-                << endPosition << " "
-                << stepFraction_ << " "
-                << celli_ << " "
+            Pout<< "Found new cell " << celli_
+                << " by searching " << checkedCells
                 << endl;
-
-            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;
-                    }
-                }
-            }
-
-            trackFraction = Cloud<ParticleType>::minValidTrackFraction;
-            position_ += trackFraction*(endPosition - position_);
         }
     }
 
@@ -599,7 +491,7 @@ void Foam::Particle<ParticleType>::trackToFaceConvex
             << facei_ << " "
             << endl;
 
-        trackFraction = Cloud<ParticleType>::minValidTrackFraction;
+        trackFraction = Cloud<ParticleType>::trackingRescueTolerance;
         position_ += trackFraction*(endPosition - position_);
     }
     else
diff --git a/src/lagrangian/basic/Particle/Particle.H b/src/lagrangian/basic/Particle/Particle.H
index 42cc5483ad..d5eb584429 100644
--- a/src/lagrangian/basic/Particle/Particle.H
+++ b/src/lagrangian/basic/Particle/Particle.H
@@ -179,7 +179,7 @@ protected:
         (
             const vector& testPt,
             const label celli,
-            bool includeOnFace
+            bool beingOnAFaceMeansOutside
         ) const;
 
         //- Track to cell face using face decomposition, used for

From bc119c4f2822e6c45931fd1a18d942acd91c7009 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Fri, 9 Oct 2009 18:22:07 +0100
Subject: [PATCH 11/16] Removing instrumentation statements

---
 src/lagrangian/basic/Particle/Particle.C | 102 +++++++++++------------
 1 file changed, 51 insertions(+), 51 deletions(-)

diff --git a/src/lagrangian/basic/Particle/Particle.C b/src/lagrangian/basic/Particle/Particle.C
index c107a38e13..18ec42bc1b 100644
--- a/src/lagrangian/basic/Particle/Particle.C
+++ b/src/lagrangian/basic/Particle/Particle.C
@@ -198,13 +198,13 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
         return;
     }
 
-    Pout<< nl << origProc_ << " "
-        << origId_ << " "
-        << position_ << " "
-        << endPosition << " "
-        << stepFraction_ << " "
-        << celli_
-        << endl;
+    // Pout<< nl << origProc_ << " "
+    //     << origId_ << " "
+    //     << position_ << " "
+    //     << endPosition << " "
+    //     << stepFraction_ << " "
+    //     << celli_
+    //     << endl;
 
     // The particle *must* have left the cell.
 
@@ -232,9 +232,9 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
 
     if (insideCellExact(position_, celli_, false))
     {
-        Pout<< "The particle starts inside the cell and ends up outside of it"
-            << nl << position_ << " " << position_ + deltaTrack
-            << endl;
+        // Pout<< "The particle starts inside the cell and ends up outside of it"
+        //     << nl << position_ << " " << position_ + deltaTrack
+        //     << endl;
 
         // The particle started inside the cell and finished outside
         // of it, find which face to cross
@@ -265,7 +265,7 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
             {
                 tmpLambda = inter.distance();
 
-                Pout<< facei << " " << tmpLambda << endl;
+                // Pout<< facei << " " << tmpLambda << endl;
 
                 if
                 (
@@ -300,8 +300,8 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
                         // Do not trigger a face hit and move the position
                         // towards the cell centre
 
-                        Pout<< "Hit a wall face heading the wrong way"
-                            << endl;
+                        // Pout<< "Hit a wall face heading the wrong way"
+                        //     << endl;
 
                         const point& cc = mesh.cellCentres()[celli_];
                         position_ += 1e-2*(cc - position_);
@@ -317,8 +317,8 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
                 // little without crossing the face to resolve the
                 // ambiguity.
 
-                Pout<< "Ambiguous face crossing, correcting towards cell "
-                    << "centre and not crossing face" << endl;
+                // Pout<< "Ambiguous face crossing, correcting towards cell "
+                //     << "centre and not crossing face" << endl;
 
                 const point& cc = mesh.cellCentres()[celli_];
                 position_ +=
@@ -336,10 +336,10 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
         }
         else
         {
-            Pout<< "Particle " << origProc_ << " " << origId_
-                << " started inside cell " << celli_ << " and finished outside"
-                << " of it, but did not find a face to cross"
-                << endl;
+            // Pout<< "Particle " << origProc_ << " " << origId_
+            //     << " started inside cell " << celli_ << " and finished outside"
+            //     << " of it, but did not find a face to cross"
+            //     << endl;
 
             const point& cc = mesh.cellCentres()[celli_];
             position_ += 1e-2*(cc - position_);
@@ -347,7 +347,7 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
     }
     else
     {
-        Pout<< "The particle started outside of the cell" << endl;
+        // Pout<< "The particle started outside of the cell" << endl;
 
         // Find which cell the particle should be in.
 
@@ -390,28 +390,28 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
 
         if (!found)
         {
-            Pout<< "Didn't find a new cell after searching "
-                << checkedCells << endl;
+            // Pout<< "Didn't find a new cell after searching "
+            //     << checkedCells << endl;
 
             const point& cc = mesh.cellCentres()[celli_];
             position_ += 1e-2*(cc - position_);
         }
-        else
-        {
-            Pout<< "Found new cell " << celli_
-                << " by searching " << checkedCells
-                << endl;
-        }
+        // else
+        // {
+        //     Pout<< "Found new cell " << celli_
+        //         << " by searching " << checkedCells
+        //         << endl;
+        // }
     }
 
-    Pout<< facei_ << " " << celli_ << endl;
+    // Pout<< facei_ << " " << celli_ << endl;
 
     if (facei_ > -1)
     {
         faceAction(trackFraction, endPosition, td);
     }
 
-    Pout<< facei_ << " " << celli_ << endl;
+    // Pout<< facei_ << " " << celli_ << endl;
 }
 
 
@@ -480,16 +480,16 @@ void Foam::Particle<ParticleType>::trackToFaceConvex
     }
     else if (lambdaMin <= 0.0)
     {
-        Pout<< "convex tracking recovery "
-            << origId_ << " "
-            << origProc_ << " "
-            << position_ << " "
-            << endPosition << " "
-            << stepFraction_ << " "
-            << lambdaMin << " "
-            << celli_ << " "
-            << facei_ << " "
-            << endl;
+        // Pout<< "convex tracking recovery "
+        //     << origId_ << " "
+        //     << origProc_ << " "
+        //     << position_ << " "
+        //     << endPosition << " "
+        //     << stepFraction_ << " "
+        //     << lambdaMin << " "
+        //     << celli_ << " "
+        //     << facei_ << " "
+        //     << endl;
 
         trackFraction = Cloud<ParticleType>::trackingRescueTolerance;
         position_ += trackFraction*(endPosition - position_);
@@ -521,17 +521,17 @@ void Foam::Particle<ParticleType>::trackToFaceConvex
 
     if (trackFraction < Cloud<ParticleType>::minValidTrackFraction)
     {
-        Pout<< "convex tracking error "
-            << origId_ << " "
-            << origProc_ << " "
-            << position_ << " "
-            << endPosition << " "
-            << trackFraction << " "
-            << stepFraction_ << " "
-            << lambdaMin << " "
-            << celli_ << " "
-            << facei_ << " "
-            << endl;
+        // Pout<< "convex tracking error "
+        //     << origId_ << " "
+        //     << origProc_ << " "
+        //     << position_ << " "
+        //     << endPosition << " "
+        //     << trackFraction << " "
+        //     << stepFraction_ << " "
+        //     << lambdaMin << " "
+        //     << celli_ << " "
+        //     << facei_ << " "
+        //     << endl;
 
         const polyMesh& mesh = cloud_.pMesh();
 

From b46cfdda5e904553f92c91b400de55435afc7d5f Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Fri, 9 Oct 2009 19:01:52 +0100
Subject: [PATCH 12/16] Not moving towards cell centre on ambiguous track, only
 a little down-range on the track.

---
 src/lagrangian/basic/Particle/Particle.C | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/lagrangian/basic/Particle/Particle.C b/src/lagrangian/basic/Particle/Particle.C
index 18ec42bc1b..45c14e536b 100644
--- a/src/lagrangian/basic/Particle/Particle.C
+++ b/src/lagrangian/basic/Particle/Particle.C
@@ -320,13 +320,13 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
                 // Pout<< "Ambiguous face crossing, correcting towards cell "
                 //     << "centre and not crossing face" << endl;
 
-                const point& cc = mesh.cellCentres()[celli_];
-                position_ +=
-                    Cloud<ParticleType>::trackingRescueTolerance
-                   *(cc - position_);
+                // const point& cc = mesh.cellCentres()[celli_];
+                // position_ +=
+                //     Cloud<ParticleType>::trackingRescueTolerance
+                //    *(cc - position_);
 
                 // Pout<< "Ambiguous face crossing. " << endl;
-                // correctLambda += Cloud<ParticleType>::minValidTrackFraction;
+                correctLambda += Cloud<ParticleType>::minValidTrackFraction;
 
                 facei_ = -1;
             }

From cbac388de8f30a39602e74278cc4df7340f02ba3 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Fri, 9 Oct 2009 19:22:46 +0100
Subject: [PATCH 13/16] Moving a successful out-of-cell track off of the
 crossed face by an additional minValidTrackFraction so to avoid an ambiguous
 track next time if the particle leaves the new cell.

---
 src/lagrangian/basic/Particle/Particle.C | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/lagrangian/basic/Particle/Particle.C b/src/lagrangian/basic/Particle/Particle.C
index 45c14e536b..7c0530b58a 100644
--- a/src/lagrangian/basic/Particle/Particle.C
+++ b/src/lagrangian/basic/Particle/Particle.C
@@ -326,11 +326,11 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
                 //    *(cc - position_);
 
                 // Pout<< "Ambiguous face crossing. " << endl;
-                correctLambda += Cloud<ParticleType>::minValidTrackFraction;
 
                 facei_ = -1;
             }
 
+            correctLambda += Cloud<ParticleType>::minValidTrackFraction;
             trackFraction = correctLambda;
             position_ += trackFraction*(endPosition - position_);
         }

From bf68452326269e16aed6f7216c444e430a08a33f Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Sun, 11 Oct 2009 00:01:59 +0100
Subject: [PATCH 14/16] Using rescue tolerance for all rescues.

---
 src/lagrangian/basic/Cloud/Cloud.C       |  2 +-
 src/lagrangian/basic/Particle/Particle.C | 10 +++++++---
 2 files changed, 8 insertions(+), 4 deletions(-)

diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
index e0679c2d45..7b069d5221 100644
--- a/src/lagrangian/basic/Cloud/Cloud.C
+++ b/src/lagrangian/basic/Cloud/Cloud.C
@@ -39,7 +39,7 @@ template<class ParticleType>
 const Foam::scalar Foam::Cloud<ParticleType>::minValidTrackFraction = 1e-6;
 
 template<class ParticleType>
-const Foam::scalar Foam::Cloud<ParticleType>::trackingRescueTolerance = 1e-6;
+const Foam::scalar Foam::Cloud<ParticleType>::trackingRescueTolerance = 1e-4;
 
 template<class ParticleType>
 const Foam::scalar Foam::Cloud<ParticleType>::intersectionTolerance = 0;
diff --git a/src/lagrangian/basic/Particle/Particle.C b/src/lagrangian/basic/Particle/Particle.C
index 7c0530b58a..be397937a3 100644
--- a/src/lagrangian/basic/Particle/Particle.C
+++ b/src/lagrangian/basic/Particle/Particle.C
@@ -304,7 +304,9 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
                         //     << endl;
 
                         const point& cc = mesh.cellCentres()[celli_];
-                        position_ += 1e-2*(cc - position_);
+                        position_ +=
+                            Cloud<ParticleType>::trackingRescueTolerance
+                           *(cc - position_);
 
                         facei_ = -1;
                     }
@@ -342,7 +344,8 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
             //     << endl;
 
             const point& cc = mesh.cellCentres()[celli_];
-            position_ += 1e-2*(cc - position_);
+            position_ +=
+                Cloud<ParticleType>::trackingRescueTolerance*(cc - position_);
         }
     }
     else
@@ -394,7 +397,8 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
             //     << checkedCells << endl;
 
             const point& cc = mesh.cellCentres()[celli_];
-            position_ += 1e-2*(cc - position_);
+            position_ +=
+                Cloud<ParticleType>::trackingRescueTolerance*(cc - position_);
         }
         // else
         // {

From 0cbddca2c8dad71a097748ab8bd5d3c04fcbd3ae Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Sun, 11 Oct 2009 00:02:36 +0100
Subject: [PATCH 15/16] Neater handling of non-patch face crossing.  Not
 applying a small tracking increase to patch interactions.

---
 src/lagrangian/basic/Particle/Particle.C | 41 ++++++++++++++++--------
 1 file changed, 27 insertions(+), 14 deletions(-)

diff --git a/src/lagrangian/basic/Particle/Particle.C b/src/lagrangian/basic/Particle/Particle.C
index be397937a3..952f5e357b 100644
--- a/src/lagrangian/basic/Particle/Particle.C
+++ b/src/lagrangian/basic/Particle/Particle.C
@@ -286,6 +286,9 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
         {
             if (!cloud_.internalFace(facei_))
             {
+                // For a patch face, allow a small value of lambda to
+                // ensure patch interactions occur.
+
                 label patchi = patch(facei_);
                 const polyPatch& patch = mesh.boundaryMesh()[patchi];
 
@@ -312,27 +315,37 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
                     }
                 }
             }
-            else if (correctLambda < Cloud<ParticleType>::minValidTrackFraction)
+            else
             {
-                // The particle is not far enough away from the face
-                // to decide if it is valid crossing. Let it move a
-                // little without crossing the face to resolve the
-                // ambiguity.
+                if (correctLambda < Cloud<ParticleType>::minValidTrackFraction)
+                {
+                    // The particle is not far enough away from the face
+                    // to decide if it is valid crossing. Let it move a
+                    // little without crossing the face to resolve the
+                    // ambiguity.
 
-                // Pout<< "Ambiguous face crossing, correcting towards cell "
-                //     << "centre and not crossing face" << endl;
+                    // Pout<< "Ambiguous face crossing, correcting towards cell "
+                    //     << "centre and not crossing face" << endl;
 
-                // const point& cc = mesh.cellCentres()[celli_];
-                // position_ +=
-                //     Cloud<ParticleType>::trackingRescueTolerance
-                //    *(cc - position_);
+                    // const point& cc = mesh.cellCentres()[celli_];
+                    // position_ +=
+                    //     Cloud<ParticleType>::trackingRescueTolerance
+                    //    *(cc - position_);
 
-                // Pout<< "Ambiguous face crossing. " << endl;
+                    // Pout<< "Ambiguous face crossing. " << endl;
 
-                facei_ = -1;
+                    facei_ = -1;
+                }
+
+                // If the face hit was not on a wall, add a small
+                // amount to the track to move it off the face, If it
+                // was not an ambiguous face crossing, this makes sure
+                // the face is not ambiguous next tracking step.  If
+                // it was ambiguous, this should resolve it.
+
+                correctLambda += Cloud<ParticleType>::minValidTrackFraction;
             }
 
-            correctLambda += Cloud<ParticleType>::minValidTrackFraction;
             trackFraction = correctLambda;
             position_ += trackFraction*(endPosition - position_);
         }

From 0e5b1affdb7401902acfbbce2d1cd65b1b49ec91 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Mon, 12 Oct 2009 15:06:48 +0100
Subject: [PATCH 16/16] Removing the 'intrudes into owner/neighbour' data
 members.

---
 src/lagrangian/basic/Cloud/Cloud.C | 48 ------------------------------
 src/lagrangian/basic/Cloud/Cloud.H | 13 --------
 2 files changed, 61 deletions(-)

diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
index 7b069d5221..babfb33298 100644
--- a/src/lagrangian/basic/Cloud/Cloud.C
+++ b/src/lagrangian/basic/Cloud/Cloud.C
@@ -62,9 +62,6 @@ bool Foam::Cloud<ParticleType>::isConcaveCell(const label cellI) const
     const pointField& cCentres = pMesh().cellCentres();
     const pointField& points = pMesh().points();
 
-    PackedBoolList& intrudesIntoOwner = intrudesIntoOwnerPtr_();
-    PackedBoolList& intrudesIntoNeighbour = intrudesIntoNeighbourPtr_();
-
     const cell& cFaces = cells[cellI];
 
     bool concave = false;
@@ -130,14 +127,6 @@ bool Foam::Cloud<ParticleType>::isConcaveCell(const label cellI) const
                     if (d < 0)
                     {
                         // Concave face
-                        if (fOwner[f0I] == cellI)
-                        {
-                            intrudesIntoOwner[f0I] = 1;
-                        }
-                        else
-                        {
-                            intrudesIntoNeighbour[f0I] = 1;
-                        }
 
                         concave = true;
                     }
@@ -159,14 +148,6 @@ bool Foam::Cloud<ParticleType>::isConcaveCell(const label cellI) const
                 if ((fn0 & fn1) > planarCosAngle)
                 {
                     // Planar face
-                    if (fOwner[f0I] == cellI)
-                    {
-                        intrudesIntoOwner[f0I] = 1;
-                    }
-                    else
-                    {
-                        intrudesIntoNeighbour[f0I] = 1;
-                    }
 
                     concave = true;
                 }
@@ -186,9 +167,6 @@ void Foam::Cloud<ParticleType>::calcConcaveCells() const
     concaveCellPtr_.reset(new PackedBoolList(pMesh().nCells()));
     PackedBoolList& concaveCell = concaveCellPtr_();
 
-    intrudesIntoOwnerPtr_.reset(new PackedBoolList(pMesh().nFaces()));
-    intrudesIntoNeighbourPtr_.reset(new PackedBoolList(pMesh().nFaces()));
-
     forAll(cells, cellI)
     {
         // if (isConcaveCell(cellI))
@@ -272,8 +250,6 @@ template<class ParticleType>
 void Foam::Cloud<ParticleType>::clearOut()
 {
     concaveCellPtr_.clear();
-    intrudesIntoOwnerPtr_.clear();
-    intrudesIntoNeighbourPtr_.clear();
     labels_.clearStorage();
 }
 
@@ -293,30 +269,6 @@ const
 }
 
 
-template<class ParticleType>
-const Foam::PackedBoolList& Foam::Cloud<ParticleType>::intrudesIntoOwner()
-const
-{
-    if (!intrudesIntoOwnerPtr_.valid())
-    {
-        calcConcaveCells();
-    }
-    return intrudesIntoOwnerPtr_();
-}
-
-
-template<class ParticleType>
-const Foam::PackedBoolList& Foam::Cloud<ParticleType>::intrudesIntoNeighbour()
-const
-{
-    if (!intrudesIntoNeighbourPtr_.valid())
-    {
-        calcConcaveCells();
-    }
-    return intrudesIntoNeighbourPtr_();
-}
-
-
 template<class ParticleType>
 Foam::label Foam::Cloud<ParticleType>::getNewParticleID() const
 {
diff --git a/src/lagrangian/basic/Cloud/Cloud.H b/src/lagrangian/basic/Cloud/Cloud.H
index b151827d0a..e30204abd6 100644
--- a/src/lagrangian/basic/Cloud/Cloud.H
+++ b/src/lagrangian/basic/Cloud/Cloud.H
@@ -81,12 +81,6 @@ class Cloud
         //- Is the cell concave
         mutable autoPtr<PackedBoolList> concaveCellPtr_;
 
-        //- Does face intrude into owner cell
-        mutable autoPtr<PackedBoolList> intrudesIntoOwnerPtr_;
-
-        //- Does face intrude into neighbour cell
-        mutable autoPtr<PackedBoolList> intrudesIntoNeighbourPtr_;
-
         //- Overall count of particles ever created. Never decreases.
         mutable label particleCount_;
 
@@ -227,13 +221,6 @@ public:
             //- Whether each cell is concave (demand driven data)
             const PackedBoolList& concaveCell() const;
 
-            //- Per face whether it intrudes into the owner cell.(demand driven)
-            const PackedBoolList& intrudesIntoOwner() const;
-
-            //- Per face whether it intrudes into the neighbour cell. (demand
-            //  driven.
-            const PackedBoolList& intrudesIntoNeighbour() const;
-
 
             // Iterators