Merge branch 'concaveTrackingByCell' into particleInteractions

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());
             }

src/lagrangian/basic/Cloud/Cloud.C

@@ -31,6 +31,181 @@ 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;
+        // }
+
+        concaveCell[cellI] = 1;
+    }
+
+    {
+        // Write cells that are a problem 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  * * * * * * * * * * * * * * //
 
@@ -38,12 +213,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 +232,43 @@ 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()
+{
+    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
 {

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,11 @@ 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_;
 
@@ -84,6 +90,12 @@ 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);
 
@@ -104,13 +116,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 +148,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 +157,8 @@ public:
         Cloud
         (
             const polyMesh& mesh,
-            const bool checkClass = true
+            const bool checkClass = true,
+            const bool concaveCheck = true
         );
 
 
@@ -136,10 +168,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 +218,9 @@ public:
                 return IDLList<ParticleType>::size();
             };
 
+            //- Whether each cell is concave (demand driven data)
+            const PackedBoolList& concaveCell() const;
+
 
             // Iterators
 

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,11 +80,13 @@ 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);

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
 {
@@ -47,10 +47,11 @@ 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);
+        scalar lam = lambda(C, endPosition, facei);
 
         if ((lam > 0) && (lam < 1.0))
         {
@@ -63,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
@@ -77,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))
         {
@@ -87,6 +88,565 @@ 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())
+        {
+            // Pout<< "insideCellExact cell " << celli
+            //     << " face " << facei << " "
+            //     << inter.distance() << endl;
+
+            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>::trackToFaceConcave
+(
+    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.
+
+    // 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.
+
+        position_ = endPosition;
+        trackFraction = 1.0;
+        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 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_;
+    vector deltaTrack =
+        mag(mesh.cellCentres()[celli_] - position_)
+       *deltaPosition/(mag(deltaPosition) + VSMALL);
+
+    // Pout<< "Inside test:" << endl;
+
+    if (insideCellExact(position_, celli_, false))
+    {
+        // 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
+
+        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();
+
+                // 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)
+        {
+            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];
+
+                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_ +=
+                            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 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;
+
+                    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;
+            }
+
+            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_ +=
+                Cloud<ParticleType>::trackingRescueTolerance*(cc - position_);
+        }
+    }
+    else
+    {
+        // Pout<< "The particle started outside of the cell" << endl;
+
+        // Find which cell the particle 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)
+        {
+            // Pout<< "Didn't find a new cell after searching "
+            //     << checkedCells << endl;
+
+            const point& cc = mesh.cellCentres()[celli_];
+            position_ +=
+                Cloud<ParticleType>::trackingRescueTolerance*(cc - position_);
+        }
+        // else
+        // {
+        //     Pout<< "Found new cell " << celli_
+        //         << " by searching " << checkedCells
+        //         << endl;
+        // }
+    }
+
+    // Pout<< facei_ << " " << celli_ << endl;
+
+    if (facei_ > -1)
+    {
+        faceAction(trackFraction, endPosition, td);
+    }
+
+    // Pout<< facei_ << " " << celli_ << endl;
+}
+
+
+template<class ParticleType>
+template<class TrackData>
+void Foam::Particle<ParticleType>::trackToFaceConvex
+(
+    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];
+            }
+        }
+    }
+
+    // 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>::trackingRescueTolerance;
+        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 < 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>
 void Foam::Particle<ParticleType>::prepareForParallelTransfer
@@ -217,7 +777,6 @@ Foam::label Foam::Particle<ParticleType>::track
 }
 
 
-
 template<class ParticleType>
 Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
 {
@@ -225,6 +784,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
@@ -233,172 +793,31 @@ 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 (faces.empty())  // inside cell
+    if (cloud_.concaveCheck_)
     {
-        trackFraction = 1.0;
+        if (cloud_.concaveCell()[celli_])
-        position_ = endPosition;
-    }
-    else // hit face
-    {
-        scalar lambdaMin = GREAT;
-
-        if (faces.size() == 1)
         {
-            lambdaMin = lambda(position_, endPosition, faces[0], stepFraction_);
+            // Use a more careful tracking algorithm if the cell is concave
-            facei_ = faces[0];
+            trackToFaceConcave(trackFraction, endPosition, td);
         }
         else
         {
-            // If the particle has to cross more than one cell to reach the
+            // Use a more careful tracking algorithm if the cell is concave
-            // endPosition, we check which way to go.
+            trackToFaceConvex(trackFraction, endPosition, td);
-            // 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)
     {
-        position_ += 1.0e-3*(mesh.cellCentres()[celli_] - position_);
+        // Use the original tracking algorithm if the cell is convex
+        trackToFaceConvex(trackFraction, endPosition, td);
     }
 
     return trackFraction;
 }
 
+
 template<class ParticleType>
 Foam::scalar Foam::Particle<ParticleType>::trackToFace
 (
@@ -409,6 +828,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     return trackToFace(endPosition, dummyTd);
 }
 
+
 template<class ParticleType>
 void Foam::Particle<ParticleType>::transformPosition(const tensor& T)
 {

src/lagrangian/basic/Particle/Particle.H

@@ -157,22 +157,60 @@ 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
         ) 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 trackToFaceConcave
+        (
+            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 trackToFaceConvex
+        (
+            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
 

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)/thermoCloud/thermoCloud.C
-$(BASECLOUDS)/reactingCloud/reactingCloud.C
+// $(BASECLOUDS)/reactingCloud/reactingCloud.C
-$(BASECLOUDS)/reactingMultiphaseCloud/reactingMultiphaseCloud.C
+// $(BASECLOUDS)/reactingMultiphaseCloud/reactingMultiphaseCloud.C
 
 
 /* kinematic parcel sub-models */
@@ -33,27 +33,27 @@ $(INTERACTINGKINEMATICPARCEL)/makeBasicInteractingKinematicParcelSubmodels.C
 
 
 /* thermo parcel sub-models */
-THERMOPARCEL=$(DERIVEDPARCELS)/basicThermoParcel
+// THERMOPARCEL=$(DERIVEDPARCELS)/basicThermoParcel
-$(THERMOPARCEL)/basicThermoParcel.C
+// $(THERMOPARCEL)/basicThermoParcel.C
-$(THERMOPARCEL)/defineBasicThermoParcel.C
+// $(THERMOPARCEL)/defineBasicThermoParcel.C
-$(THERMOPARCEL)/makeBasicThermoParcelSubmodels.C
+// $(THERMOPARCEL)/makeBasicThermoParcelSubmodels.C
 
 
 /* reacting parcel sub-models */
-REACTINGPARCEL=$(DERIVEDPARCELS)/BasicReactingParcel
+// REACTINGPARCEL=$(DERIVEDPARCELS)/BasicReactingParcel
-$(REACTINGPARCEL)/defineBasicReactingParcel.C
+// $(REACTINGPARCEL)/defineBasicReactingParcel.C
-$(REACTINGPARCEL)/makeBasicReactingParcelSubmodels.C
+// $(REACTINGPARCEL)/makeBasicReactingParcelSubmodels.C
 
 
 /* reacting multiphase parcel sub-models */
-REACTINGMPPARCEL=$(DERIVEDPARCELS)/BasicReactingMultiphaseParcel
+// REACTINGMPPARCEL=$(DERIVEDPARCELS)/BasicReactingMultiphaseParcel
-$(REACTINGMPPARCEL)/defineBasicReactingMultiphaseParcel.C
+// $(REACTINGMPPARCEL)/defineBasicReactingMultiphaseParcel.C
-$(REACTINGMPPARCEL)/makeBasicReactingMultiphaseParcelSubmodels.C
+// $(REACTINGMPPARCEL)/makeBasicReactingMultiphaseParcelSubmodels.C
 
 
 /* bolt-on models */
-submodels/addOns/radiation/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C
+// submodels/addOns/radiation/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C
-submodels/addOns/radiation/scatter/cloudScatter/cloudScatter.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/phaseProperties.C
-phaseProperties/phaseProperties/phasePropertiesIO.C
+// phaseProperties/phaseProperties/phasePropertiesIO.C
-phaseProperties/phasePropertiesList/phasePropertiesList.C
+// phaseProperties/phasePropertiesList/phasePropertiesList.C
 
 
 LIB = $(FOAM_LIBBIN)/liblagrangianIntermediate