sampledSet: Further improvements in robustness for midPoint and midPointAndFace

src/sampling/sampledSet/face/faceOnlySet.C

@@ -56,7 +56,6 @@ bool Foam::faceOnlySet::trackToBoundary
 {
     particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
 
-    // Alias
     const point& trackPt = singleParticle.position();
 
     while(true)
@@ -76,12 +75,12 @@ bool Foam::faceOnlySet::trackToBoundary
 
         if (mag(trackPt - end_) < smallDist)
         {
-            // end reached
+            // End reached
             return false;
         }
         else if (singleParticle.onBoundary())
         {
-            // Boundary reached.
+            // Boundary reached
             return true;
         }
     }

src/sampling/sampledSet/sampledSet/sampledSet.C

@@ -47,28 +47,46 @@ Foam::label Foam::sampledSet::getBoundaryCell(const label faceI) const
 }
 
 
+Foam::label Foam::sampledSet::getNeighbourCell(const label faceI) const
+{
+    if (faceI >= mesh().nInternalFaces())
+    {
+        return mesh().faceOwner()[faceI];
+    }
+    else
+    {
+        return mesh().faceNeighbour()[faceI];
+    }
+}
+
+
 Foam::label Foam::sampledSet::pointInCell
 (
     const point& p,
     const label samplei
 ) const
 {
-    const label cellio = mesh().faceOwner()[faces_[samplei]];
+    // Collect the face owner and neighbour cells of the sample into an array
-    const label cellin = mesh().faceNeighbour()[faces_[samplei]];
+    // for convenience
-    const label celljo = mesh().faceOwner()[faces_[samplei+1]];
+    label cells[4] =
-    const label celljn = mesh().faceNeighbour()[faces_[samplei+1]];
+    {
-
+        mesh().faceOwner()[faces_[samplei]],
-    // If mid-point is in the cell including the sampled faces
+        getNeighbourCell(faces_[samplei]),
-    // include in list otherwise ignore
+        mesh().faceOwner()[faces_[samplei+1]],
+        getNeighbourCell(faces_[samplei+1])
+    };
 
+    // Find the sampled cell by checking the owners and neighbours of the
+    // sampled faces
     label cellm =
-        (cellio == celljo || cellio == celljn) ? cellio
+        (cells[0] == cells[2] || cells[0] == cells[3]) ? cells[0]
-      : (cellin == celljo || cellin == celljn) ? cellin
+      : (cells[1] == cells[2] || cells[1] == cells[3]) ? cells[1]
       : -1;
 
     if (cellm != -1)
     {
-        // Check the mid-point is in the cell otherwise ignore
+        // If found the sampled cell check the point is in the cell
+        // otherwise ignore
         if (!mesh().pointInCell(p, cellm, searchEngine_.decompMode()))
         {
            cellm = -1;
@@ -81,7 +99,19 @@ Foam::label Foam::sampledSet::pointInCell
             }
         }
     }
-    else if (debug)
+    else
+    {
+        // If the sample does not pass through a single cell check if the point
+        // is in any of the owners or neighbours otherwise ignore
+        for (label i=0; i<4; i++)
+        {
+            if (mesh().pointInCell(p, cells[i], searchEngine_.decompMode()))
+            {
+                return cells[i];
+            }
+        }
+
+        if (debug)
         {
             WarningInFunction
                 << "Could not find cell for mid-point" << nl
@@ -90,12 +120,13 @@ Foam::label Foam::sampledSet::pointInCell
                 << "  face[samplei]: " << faces_[samplei]
                 << "  pts[samplei+1]: " << operator[](samplei+1)
                 << "  face[samplei+1]: " << faces_[samplei+1]
-            << "  cellio: " << cellio
+                << "  cellio: " << cells[0]
-            << "  cellin: " << cellin
+                << "  cellin: " << cells[1]
-            << "  celljo: " << celljo
+                << "  celljo: " << cells[2]
-            << "  celljn: " << celljn
+                << "  celljn: " << cells[3]
                 << endl;
         }
+    }
 
     return cellm;
 }

src/sampling/sampledSet/sampledSet/sampledSet.H

@@ -92,6 +92,9 @@ protected:
         //- Returns cell next to boundary face
         label getBoundaryCell(const label) const;
 
+        //- Returns the neigbour cell or the owner if face in on the boundary
+        label getNeighbourCell(const label) const;
+
         //- Return the cell in which the point on the sample line
         //  resides if found otherwise return -1
         label pointInCell(const point& p, const label samplei) const;