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] 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