From e879da1341cbd34cb7aef89b385397b2a1d0727c Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Tue, 18 Jul 2017 09:19:32 +0100
Subject: [PATCH] lagrangian: Fixed infinite loops

Tracking through an inverted region of the mesh happens in a reversed
direction relative to a non-inverted region. Usually, this allows the
tracking to propagate normally, regardless of the sign of the space.
However, in rare cases, it is possible for a straight trajectory to form
a closed loop through both positive and negative regions. This causes
the tracking to loop indefinitely.

To fix this, the displacement through inverted regions has been
artifically increased by a small amount (1% at the moment). This has the
effect that the change in track fraction over the negative part of the
loop no longer exactly cancels the change over the positive part, and
the track therefore terminates.
---
 src/lagrangian/basic/particle/particle.C | 19 +++++++++++++++----
 src/lagrangian/basic/particle/particle.H | 10 ++++++++++
 2 files changed, 25 insertions(+), 4 deletions(-)

diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C
index 80d3683a00..2cef6b8f8b 100644
--- a/src/lagrangian/basic/particle/particle.C
+++ b/src/lagrangian/basic/particle/particle.C
@@ -30,6 +30,8 @@ License
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
+const Foam::scalar Foam::particle::negativeSpaceDisplacementFactor = 1.01;
+
 Foam::label Foam::particle::particleCount_ = 0;
 
 namespace Foam
@@ -685,8 +687,11 @@ Foam::scalar Foam::particle::trackToStationaryTri
             << "Start local coordinates = " << y0 << endl;
     }
 
+    // Get the factor by which the displacement is increased
+    const scalar f = detA >= 0 ? 1 : negativeSpaceDisplacementFactor;
+
     // Calculate the local tracking displacement
-    barycentric Tx1(x1 & T);
+    barycentric Tx1(f*x1 & T);
 
     if (debug)
     {
@@ -783,16 +788,22 @@ Foam::scalar Foam::particle::trackToMovingTri
     FixedList<barycentricTensor, 3> T;
     movingTetReverseTransform(fraction, centre, detA, T);
 
+    // Get the factor by which the displacement is increased
+    const scalar f = detA[0] >= 0 ? 1 : negativeSpaceDisplacementFactor;
+
+    // Get the relative global position
     const vector x0Rel = x0 - centre[0];
-    const vector x1Rel = x1 - centre[1];
+    const vector x1Rel = f*x1 - centre[1];
 
     // Form the determinant and hit equations
     cubicEqn detAEqn(sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
     const barycentric yC(1, 0, 0, 0);
-    const barycentric hitEqnA = ((x1Rel & T[2]) + detA[3]*yC)*sqr(detA[0]);
+    const barycentric hitEqnA =
+        ((x1Rel & T[2])                  + detA[3]*yC)*sqr(detA[0]);
     const barycentric hitEqnB =
         ((x1Rel & T[1]) + (x0Rel & T[2]) + detA[2]*yC)*detA[0];
-    const barycentric hitEqnC = (x1Rel & T[0]) + (x0Rel & T[1]) + detA[1]*yC;
+    const barycentric hitEqnC =
+        ((x1Rel & T[0]) + (x0Rel & T[1]) + detA[1]*yC);
     const barycentric hitEqnD = y0;
     FixedList<cubicEqn, 4> hitEqn;
     forAll(hitEqn, i)
diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H
index 7f6d486f79..cdb3b78f3f 100644
--- a/src/lagrangian/basic/particle/particle.H
+++ b/src/lagrangian/basic/particle/particle.H
@@ -91,6 +91,16 @@ class particle
         //- Size in bytes of the fields
         static const std::size_t sizeofFields_;
 
+        //- The factor by which the displacement is increased when passing
+        //  through negative space. This should be slightly bigger than one.
+        //  This is needed as a straight trajectory can form a closed loop
+        //  through regions of overlapping positive and negative space, leading
+        //  to a track which never ends. By increasing the rate of displacement
+        //  through negative regions, the change in track fraction over this
+        //  part of the loop no longer exactly cancels the change over the
+        //  positive part, and the track therefore terminates.
+        static const scalar negativeSpaceDisplacementFactor;
+
 
 public: