diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
index 9248cd4abc..ceaf1c44d4 100644
--- a/src/lagrangian/basic/Cloud/Cloud.C
+++ b/src/lagrangian/basic/Cloud/Cloud.C
@@ -30,189 +30,6 @@ 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-6;
-
-template<class ParticleType>
-const Foam::scalar Foam::Cloud<ParticleType>::trackingRescueTolerance = 1e-4;
-
-template<class ParticleType>
-const Foam::scalar Foam::Cloud<ParticleType>::intersectionTolerance = 0;
-
-template<class ParticleType>
-const Foam::scalar Foam::Cloud<ParticleType>::planarCosAngle = (1 - 1e-6);
-
-
-// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
-
-// Calculate using face properties only
-template<class ParticleType>
-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 pointField& cCentres = pMesh().cellCentres();
-    const pointField& points = pMesh().points();
-
-    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];
-
-                const face& f1 = pMesh().faces()[f1I];
-
-                // 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)
-                    {
-                        // Concave face
-
-                        concave = 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
-
-                    concave = true;
-                }
-            }
-        }
-    }
-
-    return concave;
-}
-
-
-template<class ParticleType>
-void Foam::Cloud<ParticleType>::calcConcaveCells() const
-{
-    const cellList& cells = pMesh().cells();
-
-    concaveCellPtr_.reset(new PackedBoolList(pMesh().nCells()));
-    PackedBoolList& concaveCell = concaveCellPtr_();
-
-    forAll(cells, cellI)
-    {
-        // if (isConcaveCell(cellI))
-        // {
-        //     concaveCell[cellI] = 1;
-
-        //     // TODO: extend selection to include any point connected
-        //     // cell if there are conflicts between cells using the two
-        //     // different methods
-        // }
-
-        // Force all cells to be treated exactly
-        concaveCell[cellI] = 1;
-
-        // Force all cells to be treated by planes
-        // concaveCell[cellI] = 0;
-    }
-
-    // {
-    //     // Write cells that are concave to file
-
-    //     DynamicList<label> tmpConcaveCells;
-
-    //     forAll(cells, cellI)
-    //     {
-    //         if (concaveCell[cellI])
-    //         {
-    //             tmpConcaveCells.append(cellI);
-    //         }
-    //     }
-
-    //     Pout<< "Cloud<ParticleType>::calcConcaveCells() :"
-    //         << " overall cells in mesh : " << pMesh().nCells() << nl
-    //         << "    of these concave : " << tmpConcaveCells.size() << endl;
-
-    //     fileName fName = pMesh().time().path()/"concaveCells";
-
-    //     Pout<< "    Writing " << fName.name() << endl;
-
-    //     OFstream file(fName);
-
-    //     file << tmpConcaveCells;
-
-    //     file.flush();
-    // }
-}
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -220,14 +37,12 @@ template<class ParticleType>
 Foam::Cloud<ParticleType>::Cloud
 (
     const polyMesh& pMesh,
-    const IDLList<ParticleType>& particles,
-    const bool concaveCheck
+    const IDLList<ParticleType>& particles
 )
 :
     cloud(pMesh),
     IDLList<ParticleType>(),
     polyMesh_(pMesh),
-    concaveCheck_(concaveCheck),
     particleCount_(0)
 {
     IDLList<ParticleType>::operator=(particles);
@@ -239,43 +54,20 @@ Foam::Cloud<ParticleType>::Cloud
 (
     const polyMesh& pMesh,
     const word& cloudName,
-    const IDLList<ParticleType>& particles,
-    const bool concaveCheck
+    const IDLList<ParticleType>& particles
 )
 :
     cloud(pMesh, cloudName),
     IDLList<ParticleType>(),
     polyMesh_(pMesh),
-    concaveCheck_(concaveCheck),
     particleCount_(0)
 {
     IDLList<ParticleType>::operator=(particles);
 }
 
 
-template<class ParticleType>
-void Foam::Cloud<ParticleType>::clearOut()
-{
-    concaveCellPtr_.clear();
-    labels_.clearStorage();
-}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class ParticleType>
-const Foam::PackedBoolList& Foam::Cloud<ParticleType>::concaveCell()
-const
-{
-    if (!concaveCellPtr_.valid())
-    {
-        calcConcaveCells();
-    }
-
-    return concaveCellPtr_();
-}
-
-
 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 2c2075b8b4..0520c51869 100644
--- a/src/lagrangian/basic/Cloud/Cloud.H
+++ b/src/lagrangian/basic/Cloud/Cloud.H
@@ -39,7 +39,6 @@ SourceFiles
 #include "IDLList.H"
 #include "IOField.H"
 #include "polyMesh.H"
-#include "PackedBoolList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -75,11 +74,6 @@ class Cloud
 
         const polyMesh& polyMesh_;
 
-        const bool concaveCheck_;
-
-        //- Is the cell concave
-        mutable autoPtr<PackedBoolList> concaveCellPtr_;
-
         //- Overall count of particles ever created. Never decreases.
         mutable label particleCount_;
 
@@ -89,12 +83,6 @@ class Cloud
 
     // Private Member Functions
 
-        //- Determine if a particular cell is concave
-        bool isConcaveCell(const label cellI) const;
-
-        //- Analyse mesh for faces not being convex.
-        void calcConcaveCells() const;
-
         //- Initialise cloud on IO constructor
         void initCloud(const bool checkClass);
 
@@ -123,20 +111,6 @@ public:
 
     // Static data
 
-        //- 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;
-
         //- Name of cloud properties dictionary
         static word cloudPropertiesName;
 
@@ -147,8 +121,7 @@ public:
         Cloud
         (
             const polyMesh& mesh,
-            const IDLList<ParticleType>& particles,
-            const bool concaveCheck = true
+            const IDLList<ParticleType>& particles
         );
 
         //- Construct from mesh, cloud name, and a list of particles
@@ -156,8 +129,7 @@ public:
         (
             const polyMesh& mesh,
             const word& cloudName,
-            const IDLList<ParticleType>& particles,
-            const bool concaveCheck = true
+            const IDLList<ParticleType>& particles
         );
 
         //- Construct from mesh by reading from file
@@ -165,8 +137,7 @@ public:
         Cloud
         (
             const polyMesh& mesh,
-            const bool checkClass = true,
-            const bool concaveCheck = true
+            const bool checkClass = true
         );
 
 
@@ -176,17 +147,10 @@ public:
         (
             const polyMesh& pMesh,
             const word& cloudName,
-            const bool checkClass = true,
-            const bool concaveCheck = true
+            const bool checkClass = true
         );
 
 
-    // Destructor
-
-        //- Destroy any demand-driven data
-        void clearOut();
-
-
     // Member Functions
 
         // Access
@@ -226,9 +190,6 @@ public:
                 return IDLList<ParticleType>::size();
             };
 
-            //- 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 f4d7e29e8b..24a65c1685 100644
--- a/src/lagrangian/basic/Cloud/CloudIO.C
+++ b/src/lagrangian/basic/Cloud/CloudIO.C
@@ -135,13 +135,11 @@ template<class ParticleType>
 Foam::Cloud<ParticleType>::Cloud
 (
     const polyMesh& pMesh,
-    const bool checkClass,
-    const bool concaveCheck
+    const bool checkClass
 )
 :
     cloud(pMesh),
     polyMesh_(pMesh),
-    concaveCheck_(concaveCheck),
     particleCount_(0)
 {
     initCloud(checkClass);
@@ -153,13 +151,11 @@ Foam::Cloud<ParticleType>::Cloud
 (
     const polyMesh& pMesh,
     const word& cloudName,
-    const bool checkClass,
-    const bool concaveCheck
+    const bool checkClass
 )
 :
     cloud(pMesh, cloudName),
     polyMesh_(pMesh),
-    concaveCheck_(concaveCheck),
     particleCount_(0)
 {
     initCloud(checkClass);
diff --git a/src/lagrangian/basic/Particle/Particle.C b/src/lagrangian/basic/Particle/Particle.C
index 306dd85b38..b413211499 100644
--- a/src/lagrangian/basic/Particle/Particle.C
+++ b/src/lagrangian/basic/Particle/Particle.C
@@ -37,7 +37,7 @@ License
 template<class ParticleType>
 void Foam::Particle<ParticleType>::findFaces
 (
-    const vector& endPosition,
+    const vector& position,
     DynamicList<label>& faceList
 ) const
 {
@@ -46,11 +46,10 @@ void Foam::Particle<ParticleType>::findFaces
     const vector& C = mesh.cellCentres()[celli_];
 
     faceList.clear();
-
     forAll(faces, i)
     {
         label facei = faces[i];
-        scalar lam = lambda(C, endPosition, facei);
+        scalar lam = lambda(C, position, facei);
 
         if ((lam > 0) && (lam < 1.0))
         {
@@ -63,7 +62,7 @@ void Foam::Particle<ParticleType>::findFaces
 template<class ParticleType>
 void Foam::Particle<ParticleType>::findFaces
 (
-    const vector& endPosition,
+    const vector& position,
     const label celli,
     const scalar stepFraction,
     DynamicList<label>& faceList
@@ -77,7 +76,7 @@ void Foam::Particle<ParticleType>::findFaces
     forAll(faces, i)
     {
         label facei = faces[i];
-        scalar lam = lambda(C, endPosition, facei, stepFraction);
+        scalar lam = lambda(C, position, facei, stepFraction);
 
         if ((lam > 0) && (lam < 1.0))
         {
@@ -87,490 +86,6 @@ void Foam::Particle<ParticleType>::findFaces
 }
 
 
-template<class ParticleType>
-bool Foam::Particle<ParticleType>::insideCellExact
-(
-    const vector& testPt,
-    const label celli,
-    bool beingOnAFaceMeansOutside
-) const
-{
-    const polyMesh& mesh = cloud_.pMesh();
-    const labelList& faces = mesh.cells()[celli];
-    const vector& C = mesh.cellCentres()[celli];
-
-    label nFaceCrossings = 0;
-
-    // The vector from the cell centre to the end point
-    vector delta = testPt - C;
-
-    forAll (faces, i)
-    {
-        label facei = faces[i];
-
-        pointHit inter = mesh.faces()[facei].intersection
-        (
-            C,
-            delta,
-            mesh.faceCentres()[facei],
-            mesh.points(),
-            intersection::HALF_RAY,
-            Cloud<ParticleType>::intersectionTolerance
-        );
-
-        if (inter.hit())
-        {
-            if (beingOnAFaceMeansOutside)
-            {
-                if (inter.distance() <= 1.0)
-                {
-                    // This face was actually crossed.
-
-                    nFaceCrossings++;
-                }
-            }
-            else
-            {
-                if (inter.distance() < 1.0)
-                {
-                    // This face was actually crossed.
-
-                    nFaceCrossings++;
-                }
-            }
-        }
-    }
-
-    // if (nFaceCrossings > 1)
-    // {
-    //     Pout<< "In cell " << celli_ << " there were " << nFaceCrossings
-    //         << " face crossings detected tracking from concave cell "
-    //         << " centre to endPosition"
-    //         << 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>::trackToFaceExact
-(
-    scalar& trackFraction,
-    const vector& endPosition,
-    TrackData& td
-)
-{
-    facei_ = -1;
-
-    const polyMesh& mesh = cloud_.pMesh();
-    const labelList& faces = mesh.cells()[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.
-
-    if (insideCellExact(endPosition, celli_, false))
-    {
-        // 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 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
-    //    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.
-
-    const vector deltaPosition = endPosition - position_;
-
-    if (insideCellExact(position_, celli_, false))
-    {
-        // The particle started inside the cell and finished outside
-        // of it, find which face to cross
-
-        scalar tmpLambda = GREAT;
-        scalar correctLambda = GREAT;
-
-        forAll(faces, i)
-        {
-            label facei = faces[i];
-
-            // 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
-                (
-                    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)
-        {
-            if (cloud_.boundaryFace(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
-
-                        const point& cc = mesh.cellCentres()[celli_];
-                        position_ +=
-                            Cloud<ParticleType>::trackingRescueTolerance
-                           *(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 a valid crossing. Let
-                    // it move a little without crossing the face to
-                    // resolve the ambiguity.
-
-                    facei_ = -1;
-                }
-
-                // If the face hit was not on a patch, 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;
-            }
-
-            trackFraction = correctLambda;
-            position_ += trackFraction*(endPosition - position_);
-        }
-        else
-        {
-            // The particle started inside of the cell and finished
-            // outside of it, but did not find a face to cross.
-            // Applying a rescuing correction.
-
-            const point& cc = mesh.cellCentres()[celli_];
-            position_ +=
-                Cloud<ParticleType>::trackingRescueTolerance*(cc - position_);
-        }
-    }
-    else
-    {
-        // The particle started outside of the cell.  Find which cell
-        // it should be in.
-
-        const labelList& cPts = mesh.cellPoints(celli_);
-
-        DynamicList<label> checkedCells;
-
-        bool found = false;
-
-        forAll(cPts, cPtI)
-        {
-            label ptI = cPts[cPtI];
-
-            const labelList& pCs = mesh.pointCells(ptI);
-
-            forAll(pCs, pCI)
-            {
-                label cellI = pCs[pCI];
-
-                if (findIndex(checkedCells, cellI) == -1)
-                {
-                    checkedCells.append(cellI);
-
-                    if (insideCellExact(position_, cellI, false))
-                    {
-                        found = true;
-
-                        celli_ = cellI;
-
-                        break;
-                    }
-                }
-            }
-
-            if (found)
-            {
-                break;
-            }
-        }
-
-        if (!found)
-        {
-            // Didn't find a new cell after searching point connected
-            // cells.  Applying a rescuing correction.
-
-            const point& cc = mesh.cellCentres()[celli_];
-            position_ +=
-                Cloud<ParticleType>::trackingRescueTolerance*(cc - position_);
-        }
-    }
-
-    if (facei_ > -1)
-    {
-        faceAction(trackFraction, endPosition, td);
-    }
-}
-
-
-template<class ParticleType>
-template<class TrackData>
-void Foam::Particle<ParticleType>::trackToFacePlanes
-(
-    scalar& trackFraction,
-    const vector& endPosition,
-    TrackData& td
-)
-{
-    facei_ = -1;
-
-    DynamicList<label>& faces = cloud_.labels_;
-    findFaces(endPosition, faces);
-
-    if (faces.empty())
-    {
-        // 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];
-            }
-        }
-    }
-
-    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 && 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 (except by a small value to move
-        // it off the face if it is an internal face), only change
-        // cell.
-
-        trackFraction = Cloud<ParticleType>::minValidTrackFraction;
-        position_ += trackFraction*(endPosition - position_);
-    }
-    else
-    {
-        if (lambdaMin <= 1.0)
-        {
-            trackFraction = lambdaMin;
-            position_ += trackFraction*(endPosition - position_);
-        }
-        else
-        {
-            // 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 < SMALL)
-    {
-        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::faceAction")
-                << "face-cell 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>
 void Foam::Particle<ParticleType>::prepareForParallelTransfer
@@ -701,6 +216,7 @@ Foam::label Foam::Particle<ParticleType>::track
 }
 
 
+
 template<class ParticleType>
 Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
 {
@@ -708,7 +224,6 @@ 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
@@ -717,30 +232,170 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     TrackData& td
 )
 {
+    const polyMesh& mesh = cloud_.polyMesh_;
+
+    DynamicList<label>& faces = cloud_.labels_;
+    findFaces(endPosition, faces);
+
+    facei_ = -1;
     scalar trackFraction = 0.0;
 
-    if (cloud_.concaveCheck_)
+    if (faces.empty())  // inside cell
     {
-        if (cloud_.concaveCell()[celli_])
+        trackFraction = 1.0;
+        position_ = endPosition;
+    }
+    else // hit face
+    {
+        scalar lambdaMin = GREAT;
+
+        if (faces.size() == 1)
         {
-            // Use a more careful tracking algorithm if the cell is concave
-            trackToFaceExact(trackFraction, endPosition, td);
+            lambdaMin = lambda(position_, endPosition, faces[0], stepFraction_);
+            facei_ = faces[0];
         }
         else
         {
-            // Use the original tracking algorithm if the cell is convex
-            trackToFacePlanes(trackFraction, endPosition, td);
+            // 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
+            forAll(faces, i)
+            {
+                scalar lam =
+                    lambda(position_, endPosition, faces[i], stepFraction_);
+
+                if (lam < lambdaMin)
+                {
+                    lambdaMin = lam;
+                    facei_ = faces[i];
+                }
+            }
+        }
+
+        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_];
+            }
+            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);
+                }
+            }
         }
     }
-    else
+
+    // 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 < SMALL)
     {
-        trackToFacePlanes(trackFraction, endPosition, td);
+        position_ += 1.0e-3*(mesh.cellCentres()[celli_] - position_);
     }
 
     return trackFraction;
 }
 
-
 template<class ParticleType>
 Foam::scalar Foam::Particle<ParticleType>::trackToFace
 (
@@ -751,7 +406,6 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     return trackToFace(endPosition, dummyTd);
 }
 
-
 template<class ParticleType>
 void Foam::Particle<ParticleType>::transformPosition(const tensor& T)
 {
diff --git a/src/lagrangian/basic/Particle/Particle.H b/src/lagrangian/basic/Particle/Particle.H
index bc0385e3b8..077a6d4b8e 100644
--- a/src/lagrangian/basic/Particle/Particle.H
+++ b/src/lagrangian/basic/Particle/Particle.H
@@ -163,60 +163,22 @@ protected:
             const label facei
         ) const;
 
-        //- Find the faces between endPosition and cell centre
+        //- Find the faces between position and cell centre
         void findFaces
         (
-            const vector& endPosition,
+            const vector& position,
             DynamicList<label>& faceList
         ) const;
 
-        //- Find the faces between endPosition and cell centre
+        //- Find the faces between position and cell centre
         void findFaces
         (
-            const vector& endPosition,
+            const vector& position,
             const label celli,
             const scalar stepFraction,
             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 beingOnAFaceMeansOutside
-        ) const;
-
-        //- Track to cell face using face decomposition, used for
-        //  concave cells.
-        template<class TrackData>
-        void trackToFaceExact
-        (
-            scalar& trackFraction,
-            const vector& endPosition,
-            TrackData& td
-        );
-
-        //- Track to cell face using the infinite planes of the faces.
-        //  The cell *must* be convex.
-        template<class TrackData>
-        void trackToFacePlanes
-        (
-            scalar& trackFraction,
-            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
 
diff --git a/src/lagrangian/basic/Particle/ParticleI.H b/src/lagrangian/basic/Particle/ParticleI.H
index 58bbc827b8..d3098325eb 100644
--- a/src/lagrangian/basic/Particle/ParticleI.H
+++ b/src/lagrangian/basic/Particle/ParticleI.H
@@ -47,7 +47,7 @@ inline Foam::scalar Foam::Particle<ParticleType>::lambda
         vector Cf = mesh.faceCentres()[facei];
 
         // patch interaction
-        if (cloud_.boundaryFace(facei))
+        if (!cloud_.internalFace(facei))
         {
             label patchi = patch(facei);
             const polyPatch& patch = mesh.boundaryMesh()[patchi];
@@ -195,7 +195,7 @@ inline Foam::scalar Foam::Particle<ParticleType>::lambda
     vector Cf = mesh.faceCentres()[facei];
 
     // patch interaction
-    if (cloud_.boundaryFace(facei))
+    if (!cloud_.internalFace(facei))
     {
         label patchi = patch(facei);
         const polyPatch& patch = mesh.boundaryMesh()[patchi];
@@ -207,11 +207,10 @@ inline Foam::scalar Foam::Particle<ParticleType>::lambda
             scalar CCf = mag((C - Cf) & Sf);
             // check if distance between cell centre and face centre
             // is larger than wallImpactDistance
-
             const ParticleType& p = static_cast<const ParticleType&>(*this);
             if (CCf > p.wallImpactDistance(Sf))
             {
-                Cf -= p.wallImpactDistance(Sf)*Sf;
+                Cf -=p.wallImpactDistance(Sf)*Sf;
             }
         }
     }