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)
 {