Lagrangian: Enabled tracking through ACMI patches and minor code improvements

Particle collisions with ACMI patches are now handled. The hit detects
whether the location is within the overlap or the coupled region and
recurses, calling the hit routine appropriate for the region.

The low level tracking methods are now more consistently named. There is
now a distinction between tracking to a face and hitting it. Function
object side effects have been moved out of the base layer and into the
parcels on which they are meaningful.

src/functionObjects/field/nearWallFields/findCellParticle.C

@@ -110,7 +110,7 @@ bool Foam::findCellParticle::move
     while (td.keepParticle && !td.switchProcessor && stepFraction() < 1)
     {
         const scalar f = 1 - stepFraction();
-        trackToFace(f*(end_ - start_), f, td);
+        trackToAndHitFace(f*(end_ - start_), f, td);
     }
 
     if (stepFraction() == 1 || !td.keepParticle)

src/functionObjects/field/streamLine/streamLineParticle.C

@@ -203,7 +203,7 @@ bool Foam::streamLineParticle::move
                 dt = maxDt;
             }
 
-            trackToFace(dt*U, 0, td);
+            trackToAndHitFace(dt*U, 0, td);
 
             if
             (

src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C

@@ -59,7 +59,7 @@ bool Foam::DSMCParcel<ParcelType>::move(TrackData& td, const scalar trackTime)
         meshTools::constrainDirection(mesh, mesh.solutionD(), Utracking);
 
         const scalar f = 1 - p.stepFraction();
-        p.trackToFace(f*trackTime*Utracking, f, td);
+        p.trackToAndHitFace(f*trackTime*Utracking, f, td);
 
         if (p.onBoundaryFace() && td.keepParticle)
         {

src/lagrangian/basic/particle/particle.C

@@ -344,6 +344,7 @@ void Foam::particle::changeFace(const label tetTriI)
     {
         const label newFaceI = mesh_.cells()[celli_][cellFaceI];
         const class face& newFace = mesh_.faces()[newFaceI];
+        const label newOwner = mesh_.faceOwner()[newFaceI];
 
         // Exclude the current face
         if (tetFacei_ == newFaceI)
@@ -351,21 +352,21 @@ void Foam::particle::changeFace(const label tetTriI)
             continue;
         }
 
-        // Loop over the edges, looking for the shared one
+        // Loop over the edges, looking for the shared one. Note that we have to
-        label edgeComp = 0;
+        // match the direction of the edge as well as the end points in order to
+        // avoid false positives when dealing with coincident ACMI faces.
+        const label edgeComp = newOwner == celli_ ? -1 : +1;
         label edgeI = 0;
-        for (; edgeI < newFace.size(); ++ edgeI)
+        for
-        {
+        (
-            edgeComp = edge::compare(sharedEdge, newFace.faceEdge(edgeI));
+            ;
-
+            edgeI < newFace.size()
-            if (edgeComp != 0)
+         && edge::compare(sharedEdge, newFace.faceEdge(edgeI)) != edgeComp;
-            {
+            ++ edgeI
-                break;
+        );
-            }
-        }
 
         // If the face does not contain the edge, then move on to the next face
-        if (edgeComp == 0)
+        if (edgeI >= newFace.size())
         {
             continue;
         }
@@ -643,6 +644,17 @@ Foam::scalar Foam::particle::trackToFace
         {
             // The track has hit a face, so set the current face and return
             facei_ = tetFacei_;
+
+            /*
+            if (!mesh_.isInternalFace(facei_))
+            {
+                const label patchi = mesh_.boundaryMesh().whichPatch(facei_);
+                const polyPatch& patch = mesh_.boundaryMesh()[patchi];
+                Info<< "Tet face is on " << patch.type() << " patch #" << patchi
+                    << endl << endl;
+            }
+            */
+
             return f;
         }
         else
@@ -667,8 +679,8 @@ Foam::scalar Foam::particle::trackToStationaryTri
 
     if (debug)
     {
-        Info<< "Particle " << origId() << endl
+        Info<< "Particle " << origId() << endl << "Tracking from " << x0
-            << "Tracking from " << x0 << " to " << x0 + x1 << endl;
+            << " along " << x1 << " to " << x0 + x1 << endl;
     }
 
     // Get the tet geometry
@@ -778,8 +790,8 @@ Foam::scalar Foam::particle::trackToMovingTri
 
     if (debug)
     {
-        Info<< "Particle " << origId() << endl
+        Info<< "Particle " << origId() << endl << "Tracking from " << x0
-            << "Tracking from " << x0 << " to " << x0 + x1 << endl;
+            << " along " << x1 << " to " << x0 + x1 << endl;
     }
 
     // Get the tet geometry

src/lagrangian/basic/particle/particle.H

@@ -57,6 +57,8 @@ class particle;
 class polyPatch;
 
 class cyclicPolyPatch;
+class cyclicAMIPolyPatch;
+class cyclicACMIPolyPatch;
 class processorPolyPatch;
 class symmetryPlanePolyPatch;
 class symmetryPolyPatch;
@@ -294,10 +296,6 @@ protected:
 
     // Patch interactions
 
-        //- Overridable function to handle the particle hitting a face
-        template<class TrackData>
-        void hitFace(TrackData& td);
-
         //- Overridable function to handle the particle hitting a
         //  patch.  Executed before other patch-hitting functions.
         //  trackFraction is passed in to allow mesh motion to
@@ -343,6 +341,16 @@ protected:
             const vector& direction
         );
 
+        //- Overridable function to handle the particle hitting a
+        //  cyclicACMIPatch
+        template<class TrackData>
+        void hitCyclicACMIPatch
+        (
+            const cyclicACMIPolyPatch&,
+            TrackData& td,
+            const vector& direction
+        );
+
         //- Overridable function to handle the particle hitting a
         //  processorPatch
         template<class TrackData>
@@ -579,12 +587,21 @@ public:
             label& tetTriI
         );
 
-        //- As non-templated particle::trackToFace, but with additional
+        //- Hit the current face. If the current face is internal than this
-        //  boundary handling.
+        //  crosses into the next cell. If it is a boundary face then this will
+        //  interact the particle with the relevant patch.
         template<class TrackData>
-        void trackToFace
+        void hitFace
         (
-            const vector& displacement,
+            const vector& direction,
+            TrackData& td
+        );
+
+        //- Convenience function. Cobines trackToFace and hitFace
+        template<class TrackData>
+        void trackToAndHitFace
+        (
+            const vector& direction,
             const scalar fraction,
             TrackData& td
         );

src/lagrangian/basic/particle/particleTemplates.C

@@ -27,6 +27,7 @@ License
 
 #include "cyclicPolyPatch.H"
 #include "cyclicAMIPolyPatch.H"
+#include "cyclicACMIPolyPatch.H"
 #include "processorPolyPatch.H"
 #include "symmetryPlanePolyPatch.H"
 #include "symmetryPolyPatch.H"
@@ -163,38 +164,30 @@ void Foam::particle::writeFields(const CloudType& c)
 
 
 template<class TrackData>
-void Foam::particle::trackToFace
+void Foam::particle::hitFace
 (
-    const vector& displacement,
+    const vector& direction,
-    const scalar fraction,
     TrackData& td
 )
 {
-    // Track
-    trackToFace(displacement, fraction);
-
-    // If the track is complete, return
     if (!onFace())
     {
         return;
     }
-
+    else if (onInternalFace())
-    // Hit face/patch processing
-    typedef typename TrackData::cloudType::particleType particleType;
-    particleType& p = static_cast<particleType&>(*this);
-    p.hitFace(td);
-    if (onInternalFace())
     {
         changeCell();
     }
-    else
+    else if (onBoundaryFace())
     {
-        label origFacei = facei_;
+        typename TrackData::cloudType::particleType& p =
-        label patchi = mesh_.boundaryMesh().whichPatch(facei_);
+            static_cast<typename TrackData::cloudType::particleType&>(*this);
 
         // No action is taken for tetPti_ for tetFacei_ here. These are handled
         // by the patch interaction call or later during processor transfer.
 
+        const label origFacei = facei_;
+        label patchi = mesh_.boundaryMesh().whichPatch(facei_);
         const tetIndices faceHitTetIs(celli_, tetFacei_, tetPti_);
 
         if
@@ -245,13 +238,22 @@ void Foam::particle::trackToFace
                     static_cast<const cyclicPolyPatch&>(patch), td
                 );
             }
+            else if (isA<cyclicACMIPolyPatch>(patch))
+            {
+                p.hitCyclicACMIPatch
+                (
+                    static_cast<const cyclicACMIPolyPatch&>(patch),
+                    td,
+                    direction
+                );
+            }
             else if (isA<cyclicAMIPolyPatch>(patch))
             {
                 p.hitCyclicAMIPatch
                 (
                     static_cast<const cyclicAMIPolyPatch&>(patch),
                     td,
-                    displacement
+                    direction
                 );
             }
             else if (isA<processorPolyPatch>(patch))
@@ -278,8 +280,16 @@ void Foam::particle::trackToFace
 
 
 template<class TrackData>
-void Foam::particle::hitFace(TrackData&)
+void Foam::particle::trackToAndHitFace
-{}
+(
+    const vector& direction,
+    const scalar fraction,
+    TrackData& td
+)
+{
+    trackToFace(direction, fraction);
+    hitFace(direction, td);
+}
 
 
 template<class TrackData>
@@ -455,6 +465,46 @@ void Foam::particle::hitCyclicAMIPatch
 }
 
 
+template<class TrackData>
+void Foam::particle::hitCyclicACMIPatch
+(
+    const cyclicACMIPolyPatch& cpp,
+    TrackData& td,
+    const vector& direction
+)
+{
+    const label localFacei = cpp.whichFace(facei_);
+
+    // If the mask is within the patch tolerance at either end, then we can
+    // assume an interaction with the appropriate part of the ACMI pair.
+    const scalar mask = cpp.mask()[localFacei];
+    bool couple = mask >= 1 - cpp.tolerance();
+    bool nonOverlap = mask <= cpp.tolerance();
+
+    // If the mask is an intermediate value, then we search for a location on
+    // the other side of the AMI. If we can't find a location, then we assume
+    // that we have hit the non-overlap patch.
+    if (!couple && !nonOverlap)
+    {
+        vector pos = position();
+        couple = cpp.pointFace(localFacei, direction, pos) >= 0;
+        nonOverlap = !couple;
+    }
+
+    if (couple)
+    {
+        hitCyclicAMIPatch(cpp, td, direction);
+    }
+    else
+    {
+        // Move to the face associated with the non-overlap patch and redo the
+        // face interaction.
+        tetFacei_ = facei_ = cpp.nonOverlapPatch().start() + localFacei;
+        hitFace(direction, td);
+    }
+}
+
+
 template<class TrackData>
 void Foam::particle::hitProcessorPatch(const processorPolyPatch&, TrackData&)
 {}

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C

@@ -291,8 +291,8 @@ bool Foam::KinematicParcel<ParcelType>::move
         f = min(f, maxCo*l/max(SMALL*l, mag(s)));
         if (p.active())
         {
-            // Track to the next face
+            // Track to and hit the next face
-            p.trackToFace(f*s, f, td);
+            p.trackToAndHitFace(f*s, f, td);
         }
         else
         {
@@ -304,7 +304,6 @@ bool Foam::KinematicParcel<ParcelType>::move
             p.stepFraction() += f;
         }
 
-
         const scalar dt = (p.stepFraction() - sfrac)*trackTime;
 
         // Avoid problems with extremely small timesteps
@@ -331,6 +330,11 @@ bool Foam::KinematicParcel<ParcelType>::move
 
         p.age() += dt;
 
+        if (p.onFace())
+        {
+            td.cloud().functions().postFace(p, p.face(), td.keepParticle);
+        }
+
         td.cloud().functions().postMove(p, celli, dt, start, td.keepParticle);
     }
 
@@ -338,22 +342,6 @@ bool Foam::KinematicParcel<ParcelType>::move
 }
 
 
-template<class ParcelType>
-template<class TrackData>
-void Foam::KinematicParcel<ParcelType>::hitFace(TrackData& td)
-{
-    typename TrackData::cloudType::parcelType& p =
-        static_cast<typename TrackData::cloudType::parcelType&>(*this);
-
-    td.cloud().functions().postFace(p, p.face(), td.keepParticle);
-}
-
-
-template<class ParcelType>
-void Foam::KinematicParcel<ParcelType>::hitFace(int& td)
-{}
-
-
 template<class ParcelType>
 template<class TrackData>
 bool Foam::KinematicParcel<ParcelType>::hitPatch

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H

@@ -577,14 +577,6 @@ public:
 
         // Patch interactions
 
-            //- Overridable function to handle the particle hitting a face
-            //  without trackData
-            void hitFace(int& td);
-
-            //- Overridable function to handle the particle hitting a face
-            template<class TrackData>
-            void hitFace(TrackData& td);
-
             //- Overridable function to handle the particle hitting a patch
             //  Executed before other patch-hitting functions
             template<class TrackData>

src/lagrangian/molecularDynamics/molecule/molecule/molecule.C

@@ -88,7 +88,7 @@ bool Foam::molecule::move(molecule::trackingData& td, const scalar trackTime)
         while (td.keepParticle && !td.switchProcessor && stepFraction() < 1)
         {
             const scalar f = 1 - stepFraction();
-            trackToFace(f*trackTime*v_, f, td);
+            trackToAndHitFace(f*trackTime*v_, f, td);
         }
     }
     else if (td.part() == 2)

src/lagrangian/solidParticle/solidParticle.C

@@ -58,7 +58,7 @@ bool Foam::solidParticle::move
         const scalar sfrac = stepFraction();
 
         const scalar f = 1 - stepFraction();
-        trackToFace(f*trackTime*U_, f, td);
+        trackToAndHitFace(f*trackTime*U_, f, td);
 
         const scalar dt = (stepFraction() - sfrac)*trackTime;
 

src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C

@@ -144,7 +144,8 @@ bool Foam::trackedParticle::move
             td.maxLevel_[cell()] = max(td.maxLevel_[cell()], level_);
 
             const scalar f = 1 - stepFraction();
-            trackToFace(f*(end_ - start_), f, td);
+            const vector s = end_ - start_;
+            trackToAndHitFace(f*s, f, td);
         }
     }
 

src/sampling/sampledSet/face/faceOnlySet.C

@@ -62,7 +62,7 @@ bool Foam::faceOnlySet::trackToBoundary
     {
         point oldPoint = trackPt;
 
-        singleParticle.trackToFace(end_ - start_, 0, trackData);
+        singleParticle.trackToAndHitFace(end_ - start_, 0, trackData);
 
         trackPt = singleParticle.position();