From d7ff97f49029f83dc3b841adbbd43a1837f3a20f Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Tue, 12 Sep 2017 08:39:58 +0100
Subject: [PATCH] lagrangian: Optimised reduced-dimension constraints

The 4.x tracking enforces reduced dimensionality on the parcels by
moving them to the centre of the mesh at the start of each track,
without considering the topology. This can leave the parcel outside it's
associated tetrahedron.

The barycentric algorithm isn't tolerant to incorrect topology, so
instead of changing position, it was written to track to the mesh
centre. This worked, but effectively doubled the number of tracking
calls. This additional cost has now been removed by absorbing the
constraint displacement into the existing motion track, so that the same
number of tracking steps are performed as before.

Partially resolves bug report https://bugs.openfoam.org/view.php?id=2688
---
 .../parcels/Templates/DSMCParcel/DSMCParcel.C |  8 +++----
 src/lagrangian/basic/particle/particle.C      | 22 ++++++-------------
 src/lagrangian/basic/particle/particle.H      |  7 +++---
 .../KinematicParcel/KinematicParcel.C         |  8 +++----
 4 files changed, 19 insertions(+), 26 deletions(-)

diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C
index c5c4e56234..cdc5ec680c 100644
--- a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C
+++ b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C
@@ -53,17 +53,17 @@ bool Foam::DSMCParcel<ParcelType>::move
 
     while (td.keepParticle && !td.switchProcessor && p.stepFraction() < 1)
     {
-        // Apply correction to position for reduced-D cases
-        p.constrainToMeshCentre();
-
         Utracking = U_;
 
         // Apply correction to velocity to constrain tracking for
         // reduced-D cases
         meshTools::constrainDirection(mesh, mesh.solutionD(), Utracking);
 
+        // Deviation from the mesh centre for reduced-D cases
+        const vector d = p.deviationFromMeshCentre();
+
         const scalar f = 1 - p.stepFraction();
-        p.trackToAndHitFace(f*trackTime*Utracking, f, cloud, td);
+        p.trackToAndHitFace(f*trackTime*Utracking - d, f, cloud, td);
     }
 
     return td.keepParticle;
diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C
index e52ab0ba19..904f4c5723 100644
--- a/src/lagrangian/basic/particle/particle.C
+++ b/src/lagrangian/basic/particle/particle.C
@@ -984,25 +984,17 @@ Foam::scalar Foam::particle::trackToTri
 }
 
 
-void Foam::particle::constrainToMeshCentre()
+Foam::vector Foam::particle::deviationFromMeshCentre() const
 {
-    const Vector<label>& dirs = mesh_.geometricD();
-
-    bool isConstrained = false;
-    forAll(dirs, dirI)
+    if (cmptMin(mesh_.geometricD()) == -1)
     {
-        if (dirs[dirI] == -1)
-        {
-            isConstrained = true;
-            break;
-        }
+        vector pos = position(), posC = pos;
+        meshTools::constrainToMeshCentre(mesh_, posC);
+        return pos - posC;
     }
-
-    if (isConstrained)
+    else
     {
-        vector pos = position();
-        meshTools::constrainToMeshCentre(mesh_, pos);
-        track(pos - position(), 0);
+        return vector::zero;
     }
 }
 
diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H
index 21a8536b05..67c4dc393d 100644
--- a/src/lagrangian/basic/particle/particle.H
+++ b/src/lagrangian/basic/particle/particle.H
@@ -598,9 +598,10 @@ public:
             trackingData& td
         );
 
-        //- Set the constrained components of the particle position to the
-        //  mesh centre.
-        void constrainToMeshCentre();
+        //- Get the displacement from the mesh centre. Used to correct the
+        //  particle position in cases with reduced dimensionality. Returns a
+        //  zero vector for three-dimensional cases.
+        vector deviationFromMeshCentre() const;
 
 
     // Patch data
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index bc5d13c69e..394cf8ce2a 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -276,9 +276,6 @@ bool Foam::KinematicParcel<ParcelType>::move
 
     while (ttd.keepParticle && !ttd.switchProcessor && p.stepFraction() < 1)
     {
-        // Apply correction to position for reduced-D cases
-        p.constrainToMeshCentre();
-
         // Cache the current position, cell and step-fraction
         const point start = p.position();
         const scalar sfrac = p.stepFraction();
@@ -289,6 +286,9 @@ bool Foam::KinematicParcel<ParcelType>::move
         // Cell length scale
         const scalar l = cellLengthScale[p.cell()];
 
+        // Deviation from the mesh centre for reduced-D cases
+        const vector d = p.deviationFromMeshCentre();
+
         // Fraction of the displacement to track in this loop. This is limited
         // to ensure that the both the time and distance tracked is less than
         // maxCo times the total value.
@@ -298,7 +298,7 @@ bool Foam::KinematicParcel<ParcelType>::move
         if (p.active())
         {
             // Track to the next face
-            p.trackToFace(f*s, f);
+            p.trackToFace(f*s - d, f);
         }
         else
         {