From f3ef79b8295d5f2033287598b98fdf530dbcab2b Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Wed, 5 Jul 2017 09:26:56 +0100
Subject: [PATCH] ParticleCollector: Prevented missing and duplicate
 collections

The particle collector was collecting some particles twice due to a
tolerance extending the tracked path. This has been removed. The new
tracking algorithm does not generate the same sorts of spurious
tolerance-scale motions that the old one did, so this extension of the
tracking path is unnecessary.

Some particles were also not being collected at all as they were hitting
a diagonal of the collection polygon and registering as not having hit
either of the adjacent triangles. The hit criteria has been rewritten. A
hit now occurs when the normals of the triangles created by joining the
intersection point with the polygon edges are all in the same direction
as the overall polygon normal. This calculation is not affected by the
polygon's diagonals.

The issue was raised by, and resolved with support from, Karl Meredith
at FM Global.

This resolves bug-report https://bugs.openfoam.org/view.php?id=2595
---
 .../ParticleCollector/ParticleCollector.C     | 46 +++++++++----------
 1 file changed, 23 insertions(+), 23 deletions(-)

diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
index 3f278e8c5..a0a5fd744 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
@@ -290,9 +290,6 @@ void Foam::ParticleCollector<CloudType>::collectParcelPolygon
     const point& p2
 ) const
 {
-    label dummyNearType = -1;
-    label dummyNearLabel = -1;
-
     forAll(faces_, facei)
     {
         const label facePoint0 = faces_[facei][0];
@@ -311,24 +308,31 @@ void Foam::ParticleCollector<CloudType>::collectParcelPolygon
         // Intersection point
         const point pIntersect = p1 + (d1/(d1 - d2))*(p2 - p1);
 
-        const List<face>& tris = faceTris_[facei];
-
-        // Identify if point is within poly bounds
-        forAll(tris, triI)
+        // Identify if point is within the bounds of the face. Create triangles
+        // between the intersection point and each edge of the face. If all the
+        // triangle normals point in the same direction as the face normal, then
+        // the particle is within the face. Note that testing for pointHits on
+        // the face's decomposed triangles does not work due to ambiguity along
+        // the diagonals.
+        const face& f = faces_[facei];
+        const vector n = f.normal(points_);
+        bool inside = true;
+        for (label i = 0; i < f.size(); ++ i)
         {
-            const face& tri = tris[triI];
-            triPointRef t
-            (
-                points_[tri[0]],
-                points_[tri[1]],
-                points_[tri[2]]
-            );
-
-            if (t.classify(pIntersect, dummyNearType, dummyNearLabel))
+            const label j = f.fcIndex(i);
+            const triPointRef t(pIntersect, points_[f[i]], points_[f[j]]);
+            if ((n & t.normal()) < 0)
             {
-                hitFaceIDs_.append(facei);
+                inside = false;
+                break;
             }
         }
+
+        // Add to the list of hits
+        if (inside)
+        {
+            hitFaceIDs_.append(facei);
+        }
     }
 }
 
@@ -654,28 +658,24 @@ void Foam::ParticleCollector<CloudType>::postMove
         return;
     }
 
-    // Slightly extend end position to avoid falling within tracking tolerances
-    const point position1 = position0 + 1.0001*(p.position() - position0);
-
     hitFaceIDs_.clear();
 
     switch (mode_)
     {
         case mtPolygon:
         {
-            collectParcelPolygon(position0, position1);
+            collectParcelPolygon(position0, p.position());
             break;
         }
         case mtConcentricCircle:
         {
-            collectParcelConcentricCircles(position0, position1);
+            collectParcelConcentricCircles(position0, p.position());
             break;
         }
         default:
         {}
     }
 
-
     forAll(hitFaceIDs_, i)
     {
         label facei = hitFaceIDs_[i];