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

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;

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

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

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
         {