From 0db70bd8bdf61b5ec30e99b4a4da62a51a06c78d Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Tue, 19 Nov 2013 15:34:11 +0000
Subject: [PATCH] ENH: primitive/poly-Mesh - updated find cell routines

---
 src/OpenFOAM/meshes/polyMesh/polyMesh.C       | 64 ++++++++++++++-----
 .../primitiveMesh/primitiveMeshFindCell.C     | 12 ++--
 2 files changed, 55 insertions(+), 21 deletions(-)

diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
index 33c445f416..1f27017968 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
@@ -1351,8 +1351,10 @@ bool Foam::polyMesh::pointInCell
 
         case FACECENTRETETS:
         {
-            const point& cc = cellCentres()[cellI];
+            // only test that point is on inside of plane defined by cell face
+            // triangles
             const cell& cFaces = cells()[cellI];
+
             forAll(cFaces, cFaceI)
             {
                 label faceI = cFaces[cFaceI];
@@ -1376,31 +1378,61 @@ bool Foam::polyMesh::pointInCell
                         nextPointI = f[fp];
                     }
 
-                    if
-                    (
-                        tetPointRef
-                        (
-                            points()[nextPointI],
-                            points()[pointI],
-                            fc,
-                            cc
-                        ).inside(p)
-                    )
+                    const point& p0 = points()[pointI];
+                    const point& p1 = points()[nextPointI];
+                    const point& p2 = fc;
+
+                    vector twoFaceArea = (p1 - p0)^(p2 - p0);
+                    point centre = (p0 + p1 + p2)/3.0;
+                    vector proj = p - centre;
+
+                    if ((twoFaceArea & proj) > 0)
                     {
-                        return true;
+                        return false;
                     }
                 }
             }
-            return false;
+            return true;
         }
         break;
 
         case FACEDIAGTETS:
         {
-            label tetFaceI, tetPtI;
-            findTetFacePt(cellI, p, tetFaceI, tetPtI);
+            // only test that point is on inside of plane defined by cell face
+            // triangles
+            const cell& cFaces = cells()[cellI];
 
-            return tetFaceI != -1;
+            forAll(cFaces, cFaceI)
+            {
+                label faceI = cFaces[cFaceI];
+                const face& f = faces_[faceI];
+
+                for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
+                {
+                    // Get tetIndices of face triangle
+                    tetIndices faceTetIs
+                    (
+                        polyMeshTetDecomposition::triangleTetIndices
+                        (
+                            *this,
+                            faceI,
+                            cellI,
+                            tetPtI
+                        )
+                    );
+
+                    triPointRef faceTri = faceTetIs.faceTri(*this);
+
+                    vector proj = p - faceTri.centre();
+
+                    if ((faceTri.normal() & proj) > 0)
+                    {
+                        return false;
+                    }
+                }
+            }
+
+            return true;
         }
         break;
     }
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C
index bfa67c5574..edb261c601 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,8 +66,6 @@ bool Foam::primitiveMesh::pointInCell(const point& p, label celli) const
     const vectorField& cf = faceCentres();
     const vectorField& Sf = faceAreas();
 
-    bool inCell = true;
-
     forAll(f, facei)
     {
         label nFace = f[facei];
@@ -77,10 +75,14 @@ bool Foam::primitiveMesh::pointInCell(const point& p, label celli) const
         {
             normal = -normal;
         }
-        inCell = inCell && ((normal & proj) <= 0);
+
+        if ((normal & proj) > 0)
+        {
+            return false;
+        }
     }
 
-    return inCell;
+    return true;
 }