From 1f160083b15083f4f66ffd989d8bb0d1b53bb26b Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Thu, 24 Feb 2011 16:18:27 +0000
Subject: [PATCH] ENH: moved findCellFacePt and findTetFacePt methods from
 Cloud to polyMesh

---
 src/OpenFOAM/meshes/polyMesh/polyMesh.C | 115 +++++++++++++++++++++++-
 src/OpenFOAM/meshes/polyMesh/polyMesh.H |  26 +++++-
 2 files changed, 139 insertions(+), 2 deletions(-)

diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
index f83861aae9..da7ffdc266 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -1183,10 +1183,123 @@ void Foam::polyMesh::removeFiles(const fileName& instanceDir) const
     }
 }
 
+
 void Foam::polyMesh::removeFiles() const
 {
     removeFiles(instance());
 }
 
 
+void Foam::polyMesh::findCellFacePt
+(
+    const point& pt,
+    label& cellI,
+    label& tetFaceI,
+    label& tetPtI
+) const
+{
+    cellI = -1;
+    tetFaceI = -1;
+    tetPtI = -1;
+
+    const indexedOctree<treeDataCell>& tree = cellTree();
+
+    // Find nearest cell to the point
+
+    pointIndexHit info = tree.findNearest(pt, sqr(GREAT));
+
+    if (info.hit())
+    {
+        label nearestCellI = tree.shapes().cellLabels()[info.index()];
+
+        // Check the nearest cell to see if the point is inside.
+        findTetFacePt(nearestCellI, pt, tetFaceI, tetPtI);
+
+        if (tetFaceI != -1)
+        {
+            // Point was in the nearest cell
+
+            cellI = nearestCellI;
+
+            return;
+        }
+        else
+        {
+            // Check the other possible cells that the point may be in
+
+            labelList testCells = tree.findIndices(pt);
+
+            forAll(testCells, pCI)
+            {
+                label testCellI = tree.shapes().cellLabels()[testCells[pCI]];
+
+                if (testCellI == nearestCellI)
+                {
+                    // Don't retest the nearest cell
+
+                    continue;
+                }
+
+                // Check the test cell to see if the point is inside.
+                findTetFacePt(testCellI, pt, tetFaceI, tetPtI);
+
+                if (tetFaceI != -1)
+                {
+                    // Point was in the test cell
+
+                    cellI = testCellI;
+
+                    return;
+                }
+            }
+        }
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "void Foam::polyMesh::findCellFacePt"
+            "("
+                "const point&, "
+                "label&, "
+                "label&, "
+                "label&"
+            ") const"
+        )   << "Did not find nearest cell in search tree."
+            << abort(FatalError);
+    }
+}
+
+
+void Foam::polyMesh::findTetFacePt
+(
+    const label cellI,
+    const point& pt,
+    label& tetFaceI,
+    label& tetPtI
+) const
+{
+    const polyMesh& mesh = *this;
+
+    tetFaceI = -1;
+    tetPtI = -1;
+
+    List<tetIndices> cellTets =
+        polyMeshTetDecomposition::cellTetIndices(mesh, cellI);
+
+    forAll(cellTets, tetI)
+    {
+        const tetIndices& cellTetIs = cellTets[tetI];
+
+        if (cellTetIs.tet(mesh).inside(pt))
+        {
+            tetFaceI = cellTetIs.face();
+            tetPtI = cellTetIs.tetPt();
+
+            return;
+        }
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.H b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
index 65cf88dbe5..2be783a25e 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,6 +34,7 @@ SourceFiles
     polyMeshFromShapeMesh.C
     polyMeshIO.C
     polyMeshUpdate.C
+    polyMeshFindCell.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -391,6 +392,7 @@ public:
                 return *this;
             }
 
+
         // Mesh motion
 
             //- Is mesh moving
@@ -508,6 +510,28 @@ public:
 
             //- Remove all files from mesh instance()
             void removeFiles() const;
+
+
+        // Helper functions
+
+            //- Find the cell, tetFaceI and tetPtI for the given position
+            void findCellFacePt
+            (
+                const point& pt,
+                label& cellI,
+                label& tetFaceI,
+                label& tetPtI
+            ) const;
+
+            //- Find the tetFaceI and tetPtI for the given position in
+            //  the supplied cell, tetFaceI and tetPtI = -1 if not found
+            void findTetFacePt
+            (
+                const label cellI,
+                const point& pt,
+                label& tetFaceI,
+                label& tetPtI
+            ) const;
 };