From e2d12e73067ddc5557796f5340f5581c7393e86f Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@Germany>
Date: Mon, 28 Jul 2008 10:30:50 +0200
Subject: [PATCH] cuttingPlane : handle cutting planes that coincide with a
 mesh face

---
 src/sampling/cuttingPlane/cuttingPlane.C | 31 ++++++++++++++++++++----
 src/sampling/cuttingPlane/cuttingPlane.H |  4 +++
 2 files changed, 30 insertions(+), 5 deletions(-)

diff --git a/src/sampling/cuttingPlane/cuttingPlane.C b/src/sampling/cuttingPlane/cuttingPlane.C
index f546c31dbe..472d30e782 100644
--- a/src/sampling/cuttingPlane/cuttingPlane.C
+++ b/src/sampling/cuttingPlane/cuttingPlane.C
@@ -29,6 +29,15 @@ License
 #include "linePointRef.H"
 #include "meshTools.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+// set values for what is close to zero and what is considered to
+// be positive (and not just rounding noise)
+//! @cond localScope
+const Foam::scalar zeroish  = Foam::SMALL;
+const Foam::scalar positive = Foam::SMALL * 1E3;
+//! @endcond localScope
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 // Find cut cells
@@ -71,8 +80,8 @@ void Foam::cuttingPlane::calcCutCells
 
             if
             (
-                (dotProducts[e[0]] < 0 && dotProducts[e[1]] > 0)
-             || (dotProducts[e[1]] < 0 && dotProducts[e[0]] > 0)
+                (dotProducts[e[0]] < zeroish && dotProducts[e[1]] > positive)
+             || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive)
             )
             {
                 nCutEdges++;
@@ -116,8 +125,8 @@ Foam::labelList Foam::cuttingPlane::intersectEdges
 
         if
         (
-            (dotProducts[e[0]] < 0 && dotProducts[e[1]] > 0)
-         || (dotProducts[e[1]] < 0 && dotProducts[e[0]] > 0)
+            (dotProducts[e[0]] < zeroish && dotProducts[e[1]] > positive)
+         || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive)
         )
         {
             // Edge is cut.
@@ -126,7 +135,19 @@ Foam::labelList Foam::cuttingPlane::intersectEdges
 
             scalar alpha = lineIntersect(linePointRef(p0, p1));
 
-            dynCuttingPoints.append((1-alpha)*p0 + alpha*p1);
+            if (alpha < zeroish)
+            {
+                dynCuttingPoints.append(p0);
+            }
+            else if (alpha > 1.0)
+            {
+                dynCuttingPoints.append(p1);
+            }
+            else
+            {
+                dynCuttingPoints.append((1-alpha)*p0 + alpha*p1);
+            }
+
             edgePoint[edgeI] = dynCuttingPoints.size() - 1;
         }
     }
diff --git a/src/sampling/cuttingPlane/cuttingPlane.H b/src/sampling/cuttingPlane/cuttingPlane.H
index 32dec6908d..3e33be9da7 100644
--- a/src/sampling/cuttingPlane/cuttingPlane.H
+++ b/src/sampling/cuttingPlane/cuttingPlane.H
@@ -30,6 +30,10 @@ Description
 
     No attempt at resolving degenerate cases.
 
+Note
+    When the cutting plane coincides with a mesh face, the cell edge on the
+    positive side of the plane is taken.
+
 SourceFiles
     cuttingPlane.C