diff --git a/applications/utilities/mesh/advanced/autoRefineMesh/autoRefineMesh.C b/applications/utilities/mesh/advanced/autoRefineMesh/autoRefineMesh.C
index db7585039a..2f99404f9a 100644
--- a/applications/utilities/mesh/advanced/autoRefineMesh/autoRefineMesh.C
+++ b/applications/utilities/mesh/advanced/autoRefineMesh/autoRefineMesh.C
@@ -715,7 +715,7 @@ int main(int argc, char *argv[])
     triSurfaceSearch querySurf(surf);
 
     // Search engine on mesh. No face decomposition since mesh unwarped.
-    meshSearch queryMesh(mesh, false);
+    meshSearch queryMesh(mesh, polyMesh::FACEPLANES);
 
     // Check all 'outside' points
     forAll(outsidePts, outsideI)
diff --git a/applications/utilities/mesh/advanced/selectCells/selectCells.C b/applications/utilities/mesh/advanced/selectCells/selectCells.C
index 424826284d..ae7525f127 100644
--- a/applications/utilities/mesh/advanced/selectCells/selectCells.C
+++ b/applications/utilities/mesh/advanced/selectCells/selectCells.C
@@ -381,7 +381,7 @@ int main(int argc, char *argv[])
     (void)edgeCalc.minLen(Info);
 
     // Search engine on mesh. Face decomposition since faces might be warped.
-    meshSearch queryMesh(mesh, true);
+    meshSearch queryMesh(mesh, polyMesh::FACEDIAGTETS);
 
     // Check all 'outside' points
     forAll(outsidePts, outsideI)
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
index f8ce438b1f..ad420d096f 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
@@ -971,7 +971,7 @@ Foam::backgroundMeshDecomposition::distribute
         {
             // Map cellVertices, cellVertexIndices and cellVertexTypes
 
-            meshSearch cellSearch(mesh_);
+            meshSearch cellSearch(mesh_, polyMesh::FACEPLANES);
 
             const labelList& reverseCellMap = map().reverseCellMap();
 
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 5d4961e93e..3c7720e92a 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -1405,7 +1405,7 @@ bool Foam::conformalVoronoiMesh::distributeBackground()
             zeroGradientFvPatchScalarField::typeName
         );
 
-        meshSearch cellSearch(bMesh);
+        meshSearch cellSearch(bMesh, polyMesh::FACEPLANES);
 
         List<DynamicList<Foam::point> > cellVertices(bMesh.nCells());
         List<DynamicList<label> > cellVertexIndices(bMesh.nCells());
diff --git a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
index b34676fac5..2d75cc36ef 100644
--- a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
+++ b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
@@ -2194,7 +2194,7 @@ int main(int argc, char *argv[])
 
             label regionI = -1;
 
-            label cellI = mesh.findCell(insidePoint);
+            label cellI = mesh.findCell(insidePoint, polyMesh::FACEDIAGTETS);
 
             Info<< nl << "Found point " << insidePoint << " in cell " << cellI
                 << endl;
diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
index 50c5f9b592..1a92c9bc39 100644
--- a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
+++ b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
@@ -528,7 +528,11 @@ int main(int argc, char *argv[])
                         << "Cell number should be between 0 and "
                         << mesh.nCells()-1 << nl
                         << "On this mesh the particle should be in cell "
-                        << mesh.findCell(iter().position())
+                        <<  mesh.findCell
+                            (
+                                iter().position(),
+                                polyMesh::FACEDIAGTETS
+                            )
                         << exit(FatalError);
                 }
 
diff --git a/applications/utilities/preProcessing/mapFields/mapLagrangian.C b/applications/utilities/preProcessing/mapFields/mapLagrangian.C
index ee0e2aa4d3..f73029529e 100644
--- a/applications/utilities/preProcessing/mapFields/mapLagrangian.C
+++ b/applications/utilities/preProcessing/mapFields/mapLagrangian.C
@@ -56,7 +56,11 @@ static label findCell(const Cloud<passiveParticle>& cloud, const point& pt)
         // See if particle on face by finding nearest face and shifting
         // particle.
 
-        meshSearch meshSearcher(mesh, false);
+        meshSearch meshSearcher
+        (
+            mesh,
+            polyMesh::FACEPLANES    // no decomposition needed
+        );
 
         label faceI = meshSearcher.findNearestBoundaryFace(pt);
 
diff --git a/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C b/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C
index 9dab19abe4..1bb57a44d7 100644
--- a/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C
+++ b/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C
@@ -25,7 +25,7 @@ License
 
 #include "treeDataCell.H"
 #include "indexedOctree.H"
-#include "primitiveMesh.H"
+#include "polyMesh.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -83,13 +83,15 @@ void Foam::treeDataCell::update()
 Foam::treeDataCell::treeDataCell
 (
     const bool cacheBb,
-    const primitiveMesh& mesh,
-    const labelUList& cellLabels
+    const polyMesh& mesh,
+    const labelUList& cellLabels,
+    const polyMesh::cellRepresentation decompMode
 )
 :
     mesh_(mesh),
     cellLabels_(cellLabels),
-    cacheBb_(cacheBb)
+    cacheBb_(cacheBb),
+    decompMode_(decompMode)
 {
     update();
 }
@@ -98,13 +100,15 @@ Foam::treeDataCell::treeDataCell
 Foam::treeDataCell::treeDataCell
 (
     const bool cacheBb,
-    const primitiveMesh& mesh,
-    const Xfer<labelList>& cellLabels
+    const polyMesh& mesh,
+    const Xfer<labelList>& cellLabels,
+    const polyMesh::cellRepresentation decompMode
 )
 :
     mesh_(mesh),
     cellLabels_(cellLabels),
-    cacheBb_(cacheBb)
+    cacheBb_(cacheBb),
+    decompMode_(decompMode)
 {
     update();
 }
@@ -113,12 +117,14 @@ Foam::treeDataCell::treeDataCell
 Foam::treeDataCell::treeDataCell
 (
     const bool cacheBb,
-    const primitiveMesh& mesh
+    const polyMesh& mesh,
+    const polyMesh::cellRepresentation decompMode
 )
 :
     mesh_(mesh),
     cellLabels_(identity(mesh_.nCells())),
-    cacheBb_(cacheBb)
+    cacheBb_(cacheBb),
+    decompMode_(decompMode)
 {
     update();
 }
@@ -162,7 +168,7 @@ bool Foam::treeDataCell::contains
     const point& sample
 ) const
 {
-    return mesh_.pointInCell(sample, cellLabels_[index]);
+    return mesh_.pointInCell(sample, cellLabels_[index], decompMode_);
 }
 
 
diff --git a/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H b/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H
index 62bf4b4a4f..476266faf9 100644
--- a/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H
+++ b/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H
@@ -36,6 +36,7 @@ SourceFiles
 #ifndef treeDataCell_H
 #define treeDataCell_H
 
+#include "polyMesh.H"
 #include "treeBoundBoxList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -44,7 +45,6 @@ namespace Foam
 {
 
 // Forward declaration of classes
-class primitiveMesh;
 template<class Type> class indexedOctree;
 
 /*---------------------------------------------------------------------------*\
@@ -55,7 +55,7 @@ class treeDataCell
 {
     // Private data
 
-        const primitiveMesh& mesh_;
+        const polyMesh& mesh_;
 
         //- Subset of cells to work on
         const labelList cellLabels_;
@@ -63,6 +63,9 @@ class treeDataCell
         //- Whether to precalculate and store cell bounding box
         const bool cacheBb_;
 
+        //- How to decide if point is inside cell
+        const polyMesh::cellRepresentation decompMode_;
+
         //- cell bounding boxes (valid only if cacheBb_)
         treeBoundBoxList bbs_;
 
@@ -87,20 +90,27 @@ public:
         treeDataCell
         (
             const bool cacheBb,
-            const primitiveMesh&,
-            const labelUList&
+            const polyMesh&,
+            const labelUList&,
+            const polyMesh::cellRepresentation decompMode
         );
 
         //- Construct from mesh and subset of cells, transferring contents
         treeDataCell
         (
             const bool cacheBb,
-            const primitiveMesh&,
-            const Xfer<labelList>&
+            const polyMesh&,
+            const Xfer<labelList>&,
+            const polyMesh::cellRepresentation decompMode
         );
 
         //- Construct from mesh. Uses all cells in mesh.
-        treeDataCell(const bool cacheBb, const primitiveMesh&);
+        treeDataCell
+        (
+            const bool cacheBb,
+            const polyMesh&,
+            const polyMesh::cellRepresentation decompMode
+        );
 
 
     // Member Functions
@@ -112,11 +122,15 @@ public:
                 return cellLabels_;
             }
 
-            inline const primitiveMesh& mesh() const
+            inline const polyMesh& mesh() const
             {
                 return mesh_;
             }
 
+            inline polyMesh::cellRepresentation decompMode() const
+            {
+                return decompMode_;
+            }
 
             inline label size() const
             {
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
index 7a2f5afa96..cbc223c583 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
@@ -32,7 +32,6 @@ License
 #include "globalMeshData.H"
 #include "processorPolyPatch.H"
 #include "OSspecific.H"
-#include "demandDrivenData.H"
 #include "polyMeshTetDecomposition.H"
 #include "indexedOctree.H"
 #include "treeDataCell.H"
@@ -209,6 +208,7 @@ Foam::polyMesh::polyMesh(const IOobject& io)
     geometricD_(Vector<label>::zero),
     solutionD_(Vector<label>::zero),
     tetBasePtIsPtr_(NULL),
+    cellTreePtr_(NULL),
     pointZones_
     (
         IOobject
@@ -401,6 +401,7 @@ Foam::polyMesh::polyMesh
     geometricD_(Vector<label>::zero),
     solutionD_(Vector<label>::zero),
     tetBasePtIsPtr_(NULL),
+    cellTreePtr_(NULL),
     pointZones_
     (
         IOobject
@@ -558,6 +559,7 @@ Foam::polyMesh::polyMesh
     geometricD_(Vector<label>::zero),
     solutionD_(Vector<label>::zero),
     tetBasePtIsPtr_(NULL),
+    cellTreePtr_(NULL),
     pointZones_
     (
         IOobject
@@ -859,7 +861,7 @@ Foam::label Foam::polyMesh::nSolutionD() const
 
 const Foam::labelList& Foam::polyMesh::tetBasePtIs() const
 {
-    if (!tetBasePtIsPtr_)
+    if (tetBasePtIsPtr_.empty())
     {
         if (debug)
         {
@@ -869,13 +871,51 @@ const Foam::labelList& Foam::polyMesh::tetBasePtIs() const
                 << endl;
         }
 
-        tetBasePtIsPtr_ = new labelList
+        tetBasePtIsPtr_.reset
         (
-            polyMeshTetDecomposition::findFaceBasePts(*this)
+            new labelList
+            (
+                polyMeshTetDecomposition::findFaceBasePts(*this)
+            )
         );
     }
 
-    return *tetBasePtIsPtr_;
+    return tetBasePtIsPtr_();
+}
+
+
+const Foam::indexedOctree<Foam::treeDataCell>&
+Foam::polyMesh::cellTree() const
+{
+    if (cellTreePtr_.empty())
+    {
+        treeBoundBox overallBb(points());
+
+        Random rndGen(261782);
+
+        overallBb = overallBb.extend(rndGen, 1E-4);
+        overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+        overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+        cellTreePtr_.reset
+        (
+            new indexedOctree<treeDataCell>
+            (
+                treeDataCell
+                (
+                    false,          // not cache bb
+                    *this,
+                    FACEDIAGTETS    // use tetDecomposition for any inside test
+                ),
+                overallBb,
+                8,              // maxLevel
+                10,             // leafsize
+                5.0             // duplicity
+            )
+        );
+    }
+
+    return cellTreePtr_();
 }
 
 
@@ -911,7 +951,7 @@ void Foam::polyMesh::addPatches
     // recalculation. Problem: should really be done in removeBoundary but
     // there is some info in parallelData which might be interesting inbetween
     // removeBoundary and addPatches.
-    deleteDemandDrivenData(globalMeshDataPtr_);
+    globalMeshDataPtr_.clear();
 
     if (validBoundary)
     {
@@ -1033,7 +1073,7 @@ const Foam::labelList& Foam::polyMesh::faceNeighbour() const
 // Return old mesh motion points
 const Foam::pointField& Foam::polyMesh::oldPoints() const
 {
-    if (!oldPointsPtr_)
+    if (oldPointsPtr_.empty())
     {
         if (debug)
         {
@@ -1042,11 +1082,11 @@ const Foam::pointField& Foam::polyMesh::oldPoints() const
                 << endl;
         }
 
-        oldPointsPtr_ = new pointField(points_);
+        oldPointsPtr_.reset(new pointField(points_));
         curMotionTimeIndex_ = time().timeIndex();
     }
 
-    return *oldPointsPtr_;
+    return oldPointsPtr_();
 }
 
 
@@ -1068,8 +1108,8 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
     if (curMotionTimeIndex_ != time().timeIndex())
     {
         // Mesh motion in the new time step
-        deleteDemandDrivenData(oldPointsPtr_);
-        oldPointsPtr_ = new pointField(points_);
+        oldPointsPtr_.clear();
+        oldPointsPtr_.reset(new pointField(points_));
         curMotionTimeIndex_ = time().timeIndex();
     }
 
@@ -1099,9 +1139,9 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
     );
 
     // Adjust parallel shared points
-    if (globalMeshDataPtr_)
+    if (globalMeshDataPtr_.valid())
     {
-        globalMeshDataPtr_->movePoints(points_);
+        globalMeshDataPtr_().movePoints(points_);
     }
 
     // Force recalculation of all geometric data with new points
@@ -1141,14 +1181,14 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
 void Foam::polyMesh::resetMotion() const
 {
     curMotionTimeIndex_ = 0;
-    deleteDemandDrivenData(oldPointsPtr_);
+    oldPointsPtr_.clear();
 }
 
 
 // Return parallel info
 const Foam::globalMeshData& Foam::polyMesh::globalData() const
 {
-    if (!globalMeshDataPtr_)
+    if (globalMeshDataPtr_.empty())
     {
         if (debug)
         {
@@ -1158,10 +1198,10 @@ const Foam::globalMeshData& Foam::polyMesh::globalData() const
                 << endl;
         }
         // Construct globalMeshData using processorPatch information only.
-        globalMeshDataPtr_ = new globalMeshData(*this);
+        globalMeshDataPtr_.reset(new globalMeshData(*this));
     }
 
-    return *globalMeshDataPtr_;
+    return globalMeshDataPtr_();
 }
 
 
@@ -1308,4 +1348,111 @@ void Foam::polyMesh::findTetFacePt
 }
 
 
+bool Foam::polyMesh::pointInCell
+(
+    const point& p,
+    label cellI,
+    const cellRepresentation decompMode
+) const
+{
+    switch (decompMode)
+    {
+        case FACEPLANES:
+        {
+            return primitiveMesh::pointInCell(p, cellI);
+        }
+        break;
+
+        case FACECENTRETETS:
+        {
+            const point& cc = cellCentres()[cellI];
+            const cell& cFaces = cells()[cellI];
+            forAll(cFaces, cFaceI)
+            {
+                label faceI = cFaces[cFaceI];
+                const face& f = faces_[faceI];
+                const point& fc = faceCentres()[faceI];
+                bool isOwn = (owner_[faceI] == cellI);
+
+                forAll(f, fp)
+                {
+                    label pointI;
+                    label nextPointI;
+
+                    if (isOwn)
+                    {
+                        pointI = f[fp];
+                        nextPointI = f.nextLabel(fp);
+                    }
+                    else
+                    {
+                        pointI = f.nextLabel(fp);
+                        nextPointI = f[fp];
+                    }
+
+                    if
+                    (
+                        tetPointRef
+                        (
+                            points()[nextPointI],
+                            points()[pointI],
+                            fc,
+                            cc
+                        ).inside(p)
+                    )
+                    {
+                        return true;
+                    }
+                }
+            }
+            return false;
+        }
+        break;
+
+        case FACEDIAGTETS:
+        {
+            label tetFaceI, tetPtI;
+            findTetFacePt(cellI, p, tetFaceI, tetPtI);
+
+            return tetFaceI != -1;
+        }
+        break;
+    }
+    return false;
+}
+
+
+Foam::label Foam::polyMesh::findCell
+(
+    const point& location,
+    const cellRepresentation decompMode
+) const
+{
+    if (nCells() == 0)
+    {
+        return -1;
+    }
+
+    // Find the nearest cell centre to this location
+    label cellI = findNearestCell(location);
+
+    // If point is in the nearest cell return
+    if (pointInCell(location, cellI, decompMode))
+    {
+        return cellI;
+    }
+    else // point is not in the nearest cell so search all cells
+    {
+        for (label cellI = 0; cellI < nCells(); cellI++)
+        {
+            if (pointInCell(location, cellI, decompMode))
+            {
+                return cellI;
+            }
+        }
+        return -1;
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.H b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
index 4e97886ac8..ed67ef6a32 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
@@ -34,7 +34,6 @@ SourceFiles
     polyMeshFromShapeMesh.C
     polyMeshIO.C
     polyMeshUpdate.C
-    polyMeshFindCell.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -61,9 +60,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class globalMeshData;
 class mapPolyMesh;
 class polyMeshTetDecomposition;
+class treeDataCell;
+template<class Type> class indexedOctree;
 
 /*---------------------------------------------------------------------------*\
                           Class polyMesh Declaration
@@ -91,6 +93,15 @@ public:
             TOPO_PATCH_CHANGE
         };
 
+        //- Enumeration defining the representation of the cell for
+        //  inside checking
+        enum cellRepresentation
+        {
+            FACEPLANES,     // cell bound by planes of faces
+            FACECENTRETETS, // tet decomposition using facectr and cellctr
+            FACEDIAGTETS    // tet decomposition using face diagonal and cellctr
+        };
+
 
 private:
 
@@ -130,7 +141,10 @@ private:
             mutable Vector<label> solutionD_;
 
             //- Base point for face decomposition into tets
-            mutable labelList* tetBasePtIsPtr_;
+            mutable autoPtr<labelList> tetBasePtIsPtr_;
+
+            //- Search tree to allow spatial cell searching
+            mutable autoPtr<indexedOctree<treeDataCell> > cellTreePtr_;
 
 
         // Zoning information
@@ -146,7 +160,7 @@ private:
 
 
         //- Parallel info
-        mutable globalMeshData* globalMeshDataPtr_;
+        mutable autoPtr<globalMeshData> globalMeshDataPtr_;
 
 
         // Mesh motion related data
@@ -161,7 +175,7 @@ private:
             mutable label curMotionTimeIndex_;
 
             //- Old points (for the last mesh motion)
-            mutable pointField* oldPointsPtr_;
+            mutable autoPtr<pointField> oldPointsPtr_;
 
 
     // Private Member Functions
@@ -365,6 +379,9 @@ public:
             //- Return the tetBasePtIs
             const labelList& tetBasePtIs() const;
 
+            //- Return the cell search tree
+            const indexedOctree<treeDataCell>& cellTree() const;
+
             //- Return point zone mesh
             const pointZoneMesh& pointZones() const
             {
@@ -505,6 +522,9 @@ public:
             //- Clear primitive data (points, faces and cells)
             void clearPrimitives();
 
+            //- Clear cell tree data
+            void clearCellTree();
+
             //- Remove all files from mesh instance
             void removeFiles(const fileName& instanceDir) const;
 
@@ -532,6 +552,21 @@ public:
                 label& tetFaceI,
                 label& tetPtI
             ) const;
+
+            //- Is the point in the cell
+            bool pointInCell
+            (
+                const point&,
+                label cellI,
+                const cellRepresentation
+            ) const;
+
+            //- Find cell enclosing this location (-1 if not in mesh)
+            label findCell
+            (
+                const point&,
+                const cellRepresentation
+            ) const;
 };
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C b/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
index 941ebca0cf..8960e932ab 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
@@ -26,9 +26,10 @@ License
 #include "polyMesh.H"
 #include "primitiveMesh.H"
 #include "globalMeshData.H"
-#include "demandDrivenData.H"
 #include "pointMesh.H"
 #include "Time.H"
+#include "indexedOctree.H"
+#include "treeDataCell.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -72,7 +73,9 @@ void Foam::polyMesh::clearGeom()
     solutionD_ = Vector<label>::zero;
 
     // Remove the stored tet base points
-    deleteDemandDrivenData(tetBasePtIsPtr_);
+    tetBasePtIsPtr_.clear();
+    // Remove the cell tree
+    cellTreePtr_.clear();
 
     pointMesh::Delete(*this);
 }
@@ -91,7 +94,7 @@ void Foam::polyMesh::clearAddressing()
 
     // parallelData depends on the processorPatch ordering so force
     // recalculation
-    deleteDemandDrivenData(globalMeshDataPtr_);
+    globalMeshDataPtr_.clear();
 
     // Reset valid directions
     geometricD_ = Vector<label>::zero;
@@ -121,4 +124,10 @@ void Foam::polyMesh::clearOut()
 }
 
 
+void Foam::polyMesh::clearCellTree()
+{
+    cellTreePtr_.clear();
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C
index 6b37e60e5c..bf8c028ccc 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C
@@ -30,6 +30,9 @@ Description
 #include "Time.H"
 #include "primitiveMesh.H"
 #include "DynamicList.H"
+#include "indexedOctree.H"
+#include "treeDataCell.H"
+#include "globalMeshData.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -504,6 +507,7 @@ Foam::polyMesh::polyMesh
     geometricD_(Vector<label>::zero),
     solutionD_(Vector<label>::zero),
     tetBasePtIsPtr_(NULL),
+    cellTreePtr_(NULL),
     pointZones_
     (
         IOobject
@@ -548,6 +552,7 @@ Foam::polyMesh::polyMesh
     ),
     globalMeshDataPtr_(NULL),
     moving_(false),
+    changing_(false),
     curMotionTimeIndex_(time().timeIndex()),
     oldPointsPtr_(NULL)
 {
@@ -778,6 +783,7 @@ Foam::polyMesh::polyMesh
     geometricD_(Vector<label>::zero),
     solutionD_(Vector<label>::zero),
     tetBasePtIsPtr_(NULL),
+    cellTreePtr_(NULL),
     pointZones_
     (
         IOobject
@@ -822,6 +828,7 @@ Foam::polyMesh::polyMesh
     ),
     globalMeshDataPtr_(NULL),
     moving_(false),
+    changing_(false),
     curMotionTimeIndex_(time().timeIndex()),
     oldPointsPtr_(NULL)
 {
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C b/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C
index 4fa4ce6279..bba7367e8b 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C
@@ -45,7 +45,7 @@ void Foam::polyMesh::updateMesh(const mapPolyMesh& mpm)
     cellZones_.clearAddressing();
 
     // Update parallel data
-    if (globalMeshDataPtr_)
+    if (globalMeshDataPtr_.valid())
     {
         globalMeshDataPtr_->updateMesh();
     }
@@ -53,12 +53,12 @@ void Foam::polyMesh::updateMesh(const mapPolyMesh& mpm)
     setInstance(time().timeName());
 
     // Map the old motion points if present
-    if (oldPointsPtr_)
+    if (oldPointsPtr_.valid())
     {
         // Make a copy of the original points
-        pointField oldMotionPoints = *oldPointsPtr_;
+        pointField oldMotionPoints = oldPointsPtr_();
 
-        pointField& newMotionPoints = *oldPointsPtr_;
+        pointField& newMotionPoints = oldPointsPtr_();
 
         // Resize the list to new size
         newMotionPoints.setSize(points_.size());
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C
index 07108b67c4..6d7516d7ea 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C
@@ -25,8 +25,6 @@ License
 
 #include "primitiveMesh.H"
 #include "demandDrivenData.H"
-#include "indexedOctree.H"
-#include "treeDataCell.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -63,8 +61,6 @@ Foam::primitiveMesh::primitiveMesh()
     ppPtr_(NULL),
     cpPtr_(NULL),
 
-    cellTreePtr_(NULL),
-
     labels_(0),
 
     cellCentresPtr_(NULL),
@@ -107,8 +103,6 @@ Foam::primitiveMesh::primitiveMesh
     ppPtr_(NULL),
     cpPtr_(NULL),
 
-    cellTreePtr_(NULL),
-
     labels_(0),
 
     cellCentresPtr_(NULL),
@@ -351,36 +345,4 @@ const Foam::cellShapeList& Foam::primitiveMesh::cellShapes() const
 }
 
 
-const Foam::indexedOctree<Foam::treeDataCell>&
-Foam::primitiveMesh::cellTree() const
-{
-    if (!cellTreePtr_)
-    {
-        treeBoundBox overallBb(points());
-
-        Random rndGen(261782);
-
-        overallBb = overallBb.extend(rndGen, 1E-4);
-        overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-        overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
-        cellTreePtr_ =
-            new indexedOctree<treeDataCell>
-            (
-                treeDataCell
-                (
-                    false,      // not cache bb
-                    *this
-                ),
-                overallBb,
-                8,              // maxLevel
-                10,             // leafsize
-                5.0             // duplicity
-            );
-    }
-
-    return *cellTreePtr_;
-}
-
-
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
index 03157bc337..8c647d8454 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
@@ -69,10 +69,6 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
-class treeDataCell;
-template<class Type> class indexedOctree;
-
 /*---------------------------------------------------------------------------*\
                       Class primitiveMesh Declaration
 \*---------------------------------------------------------------------------*/
@@ -155,9 +151,6 @@ class primitiveMesh
             //- Cell-points
             mutable labelListList* cpPtr_;
 
-            //- Search tree to allow spatial tet searching
-            mutable indexedOctree<treeDataCell>* cellTreePtr_;
-
 
         // On-the-fly edge addresing storage
 
@@ -485,10 +478,6 @@ public:
                 const labelListList& cellPoints() const;
 
 
-            //- Build (if necessary) and return the cell search tree
-            const indexedOctree<treeDataCell>& cellTree() const;
-
-
             // Geometric data (raw!)
 
                 const vectorField& cellCentres() const;
@@ -821,9 +810,6 @@ public:
             //- Clear topological data
             void clearAddressing();
 
-            //- Clear cell tree data
-            void clearCellTree();
-
             //- Clear all geometry and addressing unnecessary for CFD
             void clearOut();
 };
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshClear.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshClear.C
index d3ef679c7e..438176cd39 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshClear.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshClear.C
@@ -24,11 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "primitiveMesh.H"
-#include "indexedOctree.H"
-#include "treeDataCell.H"
 #include "demandDrivenData.H"
-#include "indexedOctree.H"
-#include "treeDataCell.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -102,11 +98,6 @@ void Foam::primitiveMesh::printAllocated() const
         Pout<< "    Cell-point" << endl;
     }
 
-    if (cellTreePtr_)
-    {
-        Pout<< "    Cell-tree" << endl;
-    }
-
     // Geometry
     if (cellCentresPtr_)
     {
@@ -173,14 +164,6 @@ void Foam::primitiveMesh::clearAddressing()
     deleteDemandDrivenData(pePtr_);
     deleteDemandDrivenData(ppPtr_);
     deleteDemandDrivenData(cpPtr_);
-
-    deleteDemandDrivenData(cellTreePtr_);
-}
-
-
-void Foam::primitiveMesh::clearCellTree()
-{
-    deleteDemandDrivenData(cellTreePtr_);
 }
 
 
diff --git a/src/engine/ignition/ignitionSite.C b/src/engine/ignition/ignitionSite.C
index 3ecdd63271..758d84449d 100644
--- a/src/engine/ignition/ignitionSite.C
+++ b/src/engine/ignition/ignitionSite.C
@@ -36,7 +36,7 @@ void Foam::ignitionSite::findIgnitionCells(const fvMesh& mesh)
     const volVectorField& centres = mesh.C();
     const scalarField& vols = mesh.V();
 
-    label ignCell = mesh.findCell(location_);
+    label ignCell = mesh.findCell(location_, polyMesh::FACEDIAGTETS);
     if (ignCell == -1)
     {
         return;
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.C
index cd7394a173..c59b203618 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.C
+++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.C
@@ -107,7 +107,11 @@ void Foam::basicSource::setCellSet()
 
             forAll(points_, i)
             {
-                label cellI = mesh_.findCell(points_[i]);
+                label cellI = mesh_.findCell
+                (
+                    points_[i],
+                    polyMesh::FACEDIAGTETS
+                );
                 if (cellI >= 0)
                 {
                     selectedCells.insert(cellI);
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/timeActivatedExplicitSource/TimeActivatedExplicitSource.C b/src/finiteVolume/cfdTools/general/fieldSources/timeActivatedExplicitSource/TimeActivatedExplicitSource.C
index e9123bee34..cce9bc97ae 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/timeActivatedExplicitSource/TimeActivatedExplicitSource.C
+++ b/src/finiteVolume/cfdTools/general/fieldSources/timeActivatedExplicitSource/TimeActivatedExplicitSource.C
@@ -234,7 +234,11 @@ void Foam::TimeActivatedExplicitSource<Type>::setCellSet()
 
             forAll(points_, i)
             {
-                label cellI = mesh_.findCell(points_[i]);
+                label cellI = mesh_.findCell
+                (
+                    points_[i],
+                    polyMesh::FACEDIAGTETS
+                );
                 if (cellI >= 0)
                 {
                     selectedCells.insert(cellI);
diff --git a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
index 45c8d2edba..eb6048f939 100644
--- a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
+++ b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
@@ -76,7 +76,11 @@ void Foam::setRefCell
         else if (dict.found(refPointName))
         {
             point refPointi(dict.lookup(refPointName));
-            refCelli = field.mesh().findCell(refPointi);
+            refCelli = field.mesh().findCell
+            (
+                refPointi,
+                polyMesh::FACEDIAGTETS
+            );
             label hasRef = (refCelli >= 0 ? 1 : 0);
             label sumHasRef = returnReduce<label>(hasRef, sumOp<label>());
             if (sumHasRef != 1)
diff --git a/src/lagrangian/basic/InteractionLists/InteractionLists.C b/src/lagrangian/basic/InteractionLists/InteractionLists.C
index 572c66522e..408bfc0966 100644
--- a/src/lagrangian/basic/InteractionLists/InteractionLists.C
+++ b/src/lagrangian/basic/InteractionLists/InteractionLists.C
@@ -163,7 +163,13 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
 
     indexedOctree<treeDataCell> coupledPatchRangeTree
     (
-        treeDataCell(true, mesh_, coupledPatchRangeCells),
+        treeDataCell
+        (
+            true,                   // cache cell bb
+            mesh_,
+            coupledPatchRangeCells, // subset of mesh
+            polyMesh::FACEDIAGTETS  // consistent with tracking
+        ),
         procBbRndExt,
         8,              // maxLevel,
         10,             // leafSize,
@@ -382,7 +388,7 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
 
     indexedOctree<treeDataCell> allCellsTree
     (
-        treeDataCell(true, mesh_),
+        treeDataCell(true, mesh_, polyMesh::FACEDIAGTETS),
         procBbRndExt,
         8,              // maxLevel,
         10,             // leafSize,
diff --git a/src/lagrangian/dieselSpray/spray/findInjectorCell.H b/src/lagrangian/dieselSpray/spray/findInjectorCell.H
index c50d276b56..bf1ca82cec 100644
--- a/src/lagrangian/dieselSpray/spray/findInjectorCell.H
+++ b/src/lagrangian/dieselSpray/spray/findInjectorCell.H
@@ -9,7 +9,8 @@ if (injectorCell >= 0)
     foundCell = mesh_.pointInCell
     (
         injectionPosition,
-        injectorCell
+        injectorCell,
+        polyMesh::FACEDIAGTETS
     );
 }
 
@@ -39,7 +40,8 @@ if (!foundCell)
         foundCell = mesh_.pointInCell
         (
             injectionPosition,
-            injectorCell
+            injectorCell,
+            polyMesh::FACEDIAGTETS
         );
     }
     reduce(foundCell, orOp<bool>());
@@ -61,7 +63,8 @@ if (!foundCell)
             foundCell = mesh_.pointInCell
             (
                 injectionPosition,
-                injectorCell
+                injectorCell,
+                polyMesh::FACEDIAGTETS
             );
         }
         reduce(foundCell, orOp<bool>());
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
index 28d211397e..256ff8d166 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
@@ -205,7 +205,15 @@ bool Foam::InjectionModel<CloudType>::findCellAtPosition
         {
             position += SMALL*(cellCentres[cellI] - position);
 
-            if (this->owner().mesh().pointInCell(position, cellI))
+            if
+            (
+                this->owner().mesh().pointInCell
+                (
+                    position,
+                    cellI,
+                    polyMesh::FACEDIAGTETS
+                )
+            )
             {
                 procI = Pstream::myProcNo();
             }
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
index ba8d0f4200..653e34c1ca 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
@@ -86,7 +86,7 @@ Foam::labelList Foam::refinementParameters::findCells(const polyMesh& mesh)
     {
         const point& keepPoint = keepPoints_[i];
 
-        label localCellI = mesh.findCell(keepPoint);
+        label localCellI = mesh.findCell(keepPoint, polyMesh::FACEDIAGTETS);
 
         label globalCellI = -1;
 
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index 059e2a9d6d..71349d061e 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -1825,7 +1825,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMeshRegions
 
     label regionI = -1;
 
-    label cellI = mesh_.findCell(keepPoint);
+    label cellI = mesh_.findCell(keepPoint, polyMesh::FACEDIAGTETS);
 
     if (cellI != -1)
     {
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index a9c6badb3b..175ce9ca53 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -1248,7 +1248,7 @@ void Foam::meshRefinement::findCellZoneInsideWalk
         // Find the region containing the insidePoint
         label keepRegionI = -1;
 
-        label cellI = mesh_.findCell(insidePoint);
+        label cellI = mesh_.findCell(insidePoint, polyMesh::FACEDIAGTETS);
 
         if (cellI != -1)
         {
@@ -1418,7 +1418,7 @@ void Foam::meshRefinement::findCellZoneTopo
     // Find the region containing the keepPoint
     label keepRegionI = -1;
 
-    label cellI = mesh_.findCell(keepPoint);
+    label cellI = mesh_.findCell(keepPoint, polyMesh::FACEDIAGTETS);
 
     if (cellI != -1)
     {
@@ -1959,7 +1959,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
     // Find the region containing the keepPoint
     label keepRegionI = -1;
 
-    label cellI = mesh_.findCell(keepPoint);
+    label cellI = mesh_.findCell(keepPoint, polyMesh::FACEDIAGTETS);
 
     if (cellI != -1)
     {
diff --git a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
index 716c161084..b78495e443 100644
--- a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
+++ b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
@@ -174,7 +174,7 @@ void Foam::mappedPatchBase::findSamples
             }
 
             // Octree based search engine
-            meshSearch meshSearchEngine(mesh, false);
+            meshSearch meshSearchEngine(mesh, polyMesh::FACEDIAGTETS);
 
             forAll(samples, sampleI)
             {
@@ -291,7 +291,7 @@ void Foam::mappedPatchBase::findSamples
             }
 
             // Octree based search engine
-            meshSearch meshSearchEngine(mesh, false);
+            meshSearch meshSearchEngine(mesh, polyMesh::FACEDIAGTETS);
 
             forAll(samples, sampleI)
             {
diff --git a/src/meshTools/meshSearch/meshSearch.C b/src/meshTools/meshSearch/meshSearch.C
index 115fc3272a..713b6ddbcc 100644
--- a/src/meshTools/meshSearch/meshSearch.C
+++ b/src/meshTools/meshSearch/meshSearch.C
@@ -302,7 +302,7 @@ Foam::label Foam::meshSearch::findCellLinear(const point& location) const
 
     while ((!cellFound) && (n < mesh_.nCells()))
     {
-        if (pointInCell(location, n))
+        if (mesh_.pointInCell(location, n, cellDecompMode_))
         {
             cellFound = true;
             cellI = n;
@@ -338,7 +338,7 @@ Foam::label Foam::meshSearch::findCellWalk
         )   << "illegal seedCell:" << seedCellI << exit(FatalError);
     }
 
-    if (pointInCell(location, seedCellI))
+    if (mesh_.pointInCell(location, seedCellI, cellDecompMode_))
     {
         return seedCellI;
     }
@@ -368,7 +368,7 @@ Foam::label Foam::meshSearch::findCellWalk
                 }
 
                 // Check if this is the correct cell
-                if (pointInCell(location, cellI))
+                if (mesh_.pointInCell(location, cellI, cellDecompMode_))
                 {
                     return cellI;
                 }
@@ -500,10 +500,14 @@ Foam::vector Foam::meshSearch::offset
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
-Foam::meshSearch::meshSearch(const polyMesh& mesh, const bool faceDecomp)
+Foam::meshSearch::meshSearch
+(
+    const polyMesh& mesh,
+    const polyMesh::cellRepresentation cellDecompMode
+)
 :
     mesh_(mesh),
-    faceDecomp_(faceDecomp)
+    cellDecompMode_(cellDecompMode)
 {}
 
 
@@ -512,11 +516,11 @@ Foam::meshSearch::meshSearch
 (
     const polyMesh& mesh,
     const treeBoundBox& bb,
-    const bool faceDecomp
+    const polyMesh::cellRepresentation cellDecompMode
 )
 :
     mesh_(mesh),
-    faceDecomp_(faceDecomp)
+    cellDecompMode_(cellDecompMode)
 {
     overallBbPtr_.reset(new treeBoundBox(bb));
 }
@@ -615,8 +619,9 @@ const
             (
                 treeDataCell
                 (
-                    false,  // not cache bb
-                    mesh_
+                    false,          // not cache bb
+                    mesh_,
+                    cellDecompMode_ // cell decomposition mode for inside tests
                 ),
                 overallBbPtr_(),
                 8,              // maxLevel
@@ -630,90 +635,90 @@ const
 }
 
 
-// Is the point in the cell
-// Works by checking if there is a face inbetween the point and the cell
-// centre.
-// Check for internal uses proper face decomposition or just average normal.
-bool Foam::meshSearch::pointInCell(const point& p, label cellI) const
-{
-    if (faceDecomp_)
-    {
-        const point& ctr = mesh_.cellCentres()[cellI];
-
-        vector dir(p - ctr);
-        scalar magDir = mag(dir);
-
-        // Check if any faces are hit by ray from cell centre to p.
-        // If none -> p is in cell.
-        const labelList& cFaces = mesh_.cells()[cellI];
-
-        // Make sure half_ray does not pick up any faces on the wrong
-        // side of the ray.
-        scalar oldTol = intersection::setPlanarTol(0.0);
-
-        forAll(cFaces, i)
-        {
-            label faceI = cFaces[i];
-
-            pointHit inter = mesh_.faces()[faceI].ray
-            (
-                ctr,
-                dir,
-                mesh_.points(),
-                intersection::HALF_RAY,
-                intersection::VECTOR
-            );
-
-            if (inter.hit())
-            {
-                scalar dist = inter.distance();
-
-                if (dist < magDir)
-                {
-                    // Valid hit. Hit face so point is not in cell.
-                    intersection::setPlanarTol(oldTol);
-
-                    return false;
-                }
-            }
-        }
-
-        intersection::setPlanarTol(oldTol);
-
-        // No face inbetween point and cell centre so point is inside.
-        return true;
-    }
-    else
-    {
-        const labelList& f = mesh_.cells()[cellI];
-        const labelList& owner = mesh_.faceOwner();
-        const vectorField& cf = mesh_.faceCentres();
-        const vectorField& Sf = mesh_.faceAreas();
-
-        forAll(f, facei)
-        {
-            label nFace = f[facei];
-            vector proj = p - cf[nFace];
-            vector normal = Sf[nFace];
-            if (owner[nFace] == cellI)
-            {
-                if ((normal & proj) > 0)
-                {
-                    return false;
-                }
-            }
-            else
-            {
-                if ((normal & proj) < 0)
-                {
-                    return false;
-                }
-            }
-        }
-
-        return true;
-    }
-}
+//// Is the point in the cell
+//// Works by checking if there is a face inbetween the point and the cell
+//// centre.
+//// Check for internal uses proper face decomposition or just average normal.
+//bool Foam::meshSearch::pointInCell(const point& p, label cellI) const
+//{
+//    if (faceDecomp_)
+//    {
+//        const point& ctr = mesh_.cellCentres()[cellI];
+//
+//        vector dir(p - ctr);
+//        scalar magDir = mag(dir);
+//
+//        // Check if any faces are hit by ray from cell centre to p.
+//        // If none -> p is in cell.
+//        const labelList& cFaces = mesh_.cells()[cellI];
+//
+//        // Make sure half_ray does not pick up any faces on the wrong
+//        // side of the ray.
+//        scalar oldTol = intersection::setPlanarTol(0.0);
+//
+//        forAll(cFaces, i)
+//        {
+//            label faceI = cFaces[i];
+//
+//            pointHit inter = mesh_.faces()[faceI].ray
+//            (
+//                ctr,
+//                dir,
+//                mesh_.points(),
+//                intersection::HALF_RAY,
+//                intersection::VECTOR
+//            );
+//
+//            if (inter.hit())
+//            {
+//                scalar dist = inter.distance();
+//
+//                if (dist < magDir)
+//                {
+//                    // Valid hit. Hit face so point is not in cell.
+//                    intersection::setPlanarTol(oldTol);
+//
+//                    return false;
+//                }
+//            }
+//        }
+//
+//        intersection::setPlanarTol(oldTol);
+//
+//        // No face inbetween point and cell centre so point is inside.
+//        return true;
+//    }
+//    else
+//    {
+//        const labelList& f = mesh_.cells()[cellI];
+//        const labelList& owner = mesh_.faceOwner();
+//        const vectorField& cf = mesh_.faceCentres();
+//        const vectorField& Sf = mesh_.faceAreas();
+//
+//        forAll(f, facei)
+//        {
+//            label nFace = f[facei];
+//            vector proj = p - cf[nFace];
+//            vector normal = Sf[nFace];
+//            if (owner[nFace] == cellI)
+//            {
+//                if ((normal & proj) > 0)
+//                {
+//                    return false;
+//                }
+//            }
+//            else
+//            {
+//                if ((normal & proj) < 0)
+//                {
+//                    return false;
+//                }
+//            }
+//        }
+//
+//        return true;
+//    }
+//}
 
 
 Foam::label Foam::meshSearch::findNearestCell
diff --git a/src/meshTools/meshSearch/meshSearch.H b/src/meshTools/meshSearch/meshSearch.H
index bf39bca8ff..637e3f3e13 100644
--- a/src/meshTools/meshSearch/meshSearch.H
+++ b/src/meshTools/meshSearch/meshSearch.H
@@ -38,6 +38,7 @@ SourceFiles
 
 #include "pointIndexHit.H"
 #include "pointField.H"
+#include "polyMesh.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,7 +46,6 @@ namespace Foam
 {
 
 // Forward declaration of classes
-class polyMesh;
 class treeDataCell;
 class treeDataFace;
 template<class Type> class indexedOctree;
@@ -62,8 +62,8 @@ class meshSearch
         //- Reference to mesh
         const polyMesh& mesh_;
 
-        //- Whether to use face decomposition for all geometric tests
-        const bool faceDecomp_;
+        //- Whether to use cell decomposition for all geometric tests
+        const polyMesh::cellRepresentation cellDecompMode_;
 
         //- data bounding box
         mutable autoPtr<treeBoundBox> overallBbPtr_;
@@ -168,7 +168,11 @@ public:
 
         //- Construct from components. Constructs bb slightly bigger than
         //  mesh points bb.
-        meshSearch(const polyMesh& mesh, const bool faceDecomp = true);
+        meshSearch
+        (
+            const polyMesh& mesh,
+            const polyMesh::cellRepresentation = polyMesh::FACEDIAGTETS
+        );
 
         //- Construct with a custom bounding box. Any mesh element outside
         //  bb will not be found. Up to user to make sure bb
@@ -177,7 +181,7 @@ public:
         (
             const polyMesh& mesh,
             const treeBoundBox& bb,
-            const bool faceDecomp = true
+            const polyMesh::cellRepresentation = polyMesh::FACEDIAGTETS
         );
 
     //- Destructor
@@ -193,6 +197,11 @@ public:
                 return mesh_;
             }
 
+            polyMesh::cellRepresentation decompMode() const
+            {
+                return cellDecompMode_;
+            }
+
             //- Get (demand driven) reference to octree holding all
             //  boundary faces
             const indexedOctree<treeDataFace>& boundaryTree() const;
@@ -203,10 +212,6 @@ public:
 
         // Queries
 
-            //- test for point in cell. Does not handle cells with center
-            //  outside cell.
-            bool pointInCell(const point& p, const label celli) const;
-
             //- Find nearest cell in terms of cell centre.
             //  Options:
             //  - use octree
@@ -227,7 +232,7 @@ public:
                 const bool useTreeSearch = true
             ) const;
 
-            //- Find cell containing (using pointInCell) location.
+            //- Find cell containing location.
             //  If seed provided walks and falls back to linear/tree search.
             //  (so handles holes correctly)s
             //  Returns -1 if not in domain.
diff --git a/src/meshTools/octree/octreeDataCell.C b/src/meshTools/octree/octreeDataCell.C
index 4dcc838a0e..9ddf44f9ae 100644
--- a/src/meshTools/octree/octreeDataCell.C
+++ b/src/meshTools/octree/octreeDataCell.C
@@ -116,7 +116,12 @@ bool Foam::octreeDataCell::contains
     const point& sample
 ) const
 {
-    return mesh_.pointInCell(sample, cellLabels_[index]);
+    return mesh_.pointInCell
+    (
+        sample,
+        cellLabels_[index],
+        polyMesh::FACEDIAGTETS
+    );
 }
 
 
diff --git a/src/meshTools/sets/cellSources/nearestToCell/nearestToCell.C b/src/meshTools/sets/cellSources/nearestToCell/nearestToCell.C
index 7c64d11a7e..bf7bb5d7e5 100644
--- a/src/meshTools/sets/cellSources/nearestToCell/nearestToCell.C
+++ b/src/meshTools/sets/cellSources/nearestToCell/nearestToCell.C
@@ -25,8 +25,6 @@ License
 
 #include "nearestToCell.H"
 #include "polyMesh.H"
-#include "meshSearch.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
diff --git a/src/meshTools/sets/cellSources/regionToCell/regionToCell.C b/src/meshTools/sets/cellSources/regionToCell/regionToCell.C
index 77948ca5a5..375d9d1894 100644
--- a/src/meshTools/sets/cellSources/regionToCell/regionToCell.C
+++ b/src/meshTools/sets/cellSources/regionToCell/regionToCell.C
@@ -60,7 +60,7 @@ Foam::topoSetSource::addToUsageTable Foam::regionToCell::usage_
 
 void Foam::regionToCell::combine(topoSet& set, const bool add) const
 {
-    label cellI = mesh_.findCell(insidePoint_);
+    label cellI = mesh_.findCell(insidePoint_, polyMesh::FACEDIAGTETS);
 
     // Load the subset of cells
     boolList blockedFace(mesh_.nFaces(), false);
diff --git a/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C b/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C
index 819c6dc476..16c031f1be 100644
--- a/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C
+++ b/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C
@@ -166,7 +166,7 @@ void Foam::surfaceToCell::combine(topoSet& set, const bool add) const
 
         // Construct search engine on mesh
 
-        meshSearch queryMesh(mesh_, true);
+        meshSearch queryMesh(mesh_, polyMesh::FACEDIAGTETS);
 
 
         // Check all 'outside' points
diff --git a/src/meshTools/sets/pointSources/nearestToPoint/nearestToPoint.C b/src/meshTools/sets/pointSources/nearestToPoint/nearestToPoint.C
index 8af43613da..6ce211033f 100644
--- a/src/meshTools/sets/pointSources/nearestToPoint/nearestToPoint.C
+++ b/src/meshTools/sets/pointSources/nearestToPoint/nearestToPoint.C
@@ -25,8 +25,6 @@ License
 
 #include "nearestToPoint.H"
 #include "polyMesh.H"
-#include "meshSearch.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
diff --git a/src/meshTools/sets/pointSources/surfaceToPoint/surfaceToPoint.C b/src/meshTools/sets/pointSources/surfaceToPoint/surfaceToPoint.C
index 0a7fafd89b..4aaf9eee17 100644
--- a/src/meshTools/sets/pointSources/surfaceToPoint/surfaceToPoint.C
+++ b/src/meshTools/sets/pointSources/surfaceToPoint/surfaceToPoint.C
@@ -25,7 +25,6 @@ License
 
 #include "surfaceToPoint.H"
 #include "polyMesh.H"
-#include "meshSearch.H"
 #include "triSurfaceSearch.H"
 #include "cpuTime.H"
 
diff --git a/src/meshTools/surfaceSets/surfaceSets.C b/src/meshTools/surfaceSets/surfaceSets.C
index 47f553517a..a9689e0d2a 100644
--- a/src/meshTools/surfaceSets/surfaceSets.C
+++ b/src/meshTools/surfaceSets/surfaceSets.C
@@ -235,7 +235,7 @@ void Foam::surfaceSets::getSurfaceSets
 )
 {
     // Construct search engine on mesh
-    meshSearch queryMesh(mesh, true);
+    meshSearch queryMesh(mesh, polyMesh::FACEDIAGTETS);
 
     // Cut faces with surface and classify cells
     cellClassification cellType
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLine.C b/src/postProcessing/functionObjects/field/streamLine/streamLine.C
index 9e8896d7cc..67fe11c29a 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLine.C
+++ b/src/postProcessing/functionObjects/field/streamLine/streamLine.C
@@ -363,7 +363,7 @@ void Foam::streamLine::read(const dictionary& dict)
 
         const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
 
-        meshSearchPtr_.reset(new meshSearch(mesh, false));
+        meshSearchPtr_.reset(new meshSearch(mesh, polyMesh::FACEDIAGTETS));
 
         const dictionary& coeffsDict = dict.subDict(seedSet_ + "Coeffs");
         sampledSetPtr_ = sampledSet::New
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C b/src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C
index a1566c28b3..c93db3f316 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C
+++ b/src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C
@@ -102,7 +102,7 @@ void Foam::meshToMesh::calcAddressing()
 
     indexedOctree<treeDataCell> oc
     (
-        treeDataCell(false, fromMesh_),
+        treeDataCell(false, fromMesh_, polyMesh::FACEDIAGTETS),
         shiftedBb,      // overall bounding box
         8,              // maxLevel
         10,             // leafsize
@@ -267,7 +267,7 @@ void Foam::meshToMesh::cellAddresses
         cellAddressing_[toI] = -1;
 
         // Check point is actually in the nearest cell
-        if (fromMesh.pointInCell(p, curCell))
+        if (fromMesh.pointInCell(p, curCell, polyMesh::FACEDIAGTETS))
         {
             cellAddressing_[toI] = curCell;
         }
@@ -292,7 +292,15 @@ void Foam::meshToMesh::cellAddresses
                 {
                     // search through all the neighbours.
                     // If point is in neighbour reset current cell
-                    if (fromMesh.pointInCell(p, neighbours[nI]))
+                    if
+                    (
+                        fromMesh.pointInCell
+                        (
+                            p,
+                            neighbours[nI],
+                            polyMesh::FACEDIAGTETS
+                        )
+                    )
                     {
                         cellAddressing_[toI] = neighbours[nI];
                         found = true;
@@ -316,7 +324,15 @@ void Foam::meshToMesh::cellAddresses
                         {
                             // search through all the neighbours.
                             // If point is in neighbour reset current cell
-                            if (fromMesh.pointInCell(p, nn[nI]))
+                            if
+                            (
+                                fromMesh.pointInCell
+                                (
+                                    p,
+                                    nn[nI],
+                                    polyMesh::FACEDIAGTETS
+                                )
+                            )
                             {
                                 cellAddressing_[toI] = nn[nI];
                                 found = true;
diff --git a/src/sampling/probes/patchProbes.C b/src/sampling/probes/patchProbes.C
index 49b0bc413f..e1aaadc923 100644
--- a/src/sampling/probes/patchProbes.C
+++ b/src/sampling/probes/patchProbes.C
@@ -28,7 +28,6 @@ License
 #include "IOmanip.H"
 // For 'nearInfo' helper class only
 #include "mappedPatchBase.H"
-//#include "meshSearch.H"
 #include "treeBoundBox.H"
 #include "treeDataFace.H"
 
diff --git a/src/sampling/probes/probes.C b/src/sampling/probes/probes.C
index b0a073e1e4..8bca017d4b 100644
--- a/src/sampling/probes/probes.C
+++ b/src/sampling/probes/probes.C
@@ -45,7 +45,7 @@ void Foam::probes::findElements(const fvMesh& mesh)
     {
         const vector& location = operator[](probeI);
 
-        elementList_[probeI] = mesh.findCell(location);
+        elementList_[probeI] = mesh.findCell(location, polyMesh::FACEDIAGTETS);
 
         if (debug && elementList_[probeI] != -1)
         {
diff --git a/src/sampling/sampledSet/sampledSet/sampledSet.C b/src/sampling/sampledSet/sampledSet/sampledSet.C
index 36704bbbef..05256451b1 100644
--- a/src/sampling/sampledSet/sampledSet/sampledSet.C
+++ b/src/sampling/sampledSet/sampledSet/sampledSet.C
@@ -68,7 +68,7 @@ Foam::label Foam::sampledSet::getCell
     {
         label cellI = getBoundaryCell(faceI);
 
-        if (!mesh().pointInCell(sample, cellI))
+        if (!mesh().pointInCell(sample, cellI, searchEngine_.decompMode()))
         {
             FatalErrorIn
             (
@@ -85,7 +85,7 @@ Foam::label Foam::sampledSet::getCell
 
         label cellI = mesh().faceOwner()[faceI];
 
-        if (mesh().pointInCell(sample, cellI))
+        if (mesh().pointInCell(sample, cellI, searchEngine_.decompMode()))
         {
             return cellI;
         }
@@ -93,7 +93,7 @@ Foam::label Foam::sampledSet::getCell
         {
             cellI = mesh().faceNeighbour()[faceI];
 
-            if (mesh().pointInCell(sample, cellI))
+            if (mesh().pointInCell(sample, cellI, searchEngine_.decompMode()))
             {
                 return cellI;
             }
@@ -261,12 +261,17 @@ bool Foam::sampledSet::getTrackingPoint
     if (bFaceI == -1)
     {
         // No boundary intersection. Try and find cell samplePt is in
-        trackCellI = mesh().findCell(samplePt);
+        trackCellI = mesh().findCell(samplePt, searchEngine_.decompMode());
 
         if
         (
             (trackCellI == -1)
-         || !mesh().pointInCell(samplePt, trackCellI)
+        || !mesh().pointInCell
+            (
+                samplePt,
+                trackCellI,
+                searchEngine_.decompMode()
+            )
         )
         {
             // Line samplePt - end_ does not intersect domain at all.
@@ -316,7 +321,7 @@ bool Foam::sampledSet::getTrackingPoint
             // samplePt inside or marginally outside.
             trackPt = samplePt;
             trackFaceI = -1;
-            trackCellI = mesh().findCell(trackPt);
+            trackCellI = mesh().findCell(trackPt, searchEngine_.decompMode());
 
             isGoodSample = true;
         }
diff --git a/src/sampling/sampledSet/sampledSets/sampledSets.C b/src/sampling/sampledSet/sampledSets/sampledSets.C
index f2bbced7a5..6882649782 100644
--- a/src/sampling/sampledSet/sampledSets/sampledSets.C
+++ b/src/sampling/sampledSet/sampledSets/sampledSets.C
@@ -138,7 +138,7 @@ Foam::sampledSets::sampledSets
     mesh_(refCast<const fvMesh>(obr)),
     loadFromFiles_(loadFromFiles),
     outputPath_(fileName::null),
-    searchEngine_(mesh_, true),
+    searchEngine_(mesh_, polyMesh::FACEDIAGTETS),
     interpolationScheme_(word::null),
     writeFormat_(word::null)
 {
diff --git a/src/sampling/sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.C b/src/sampling/sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.C
index 3704f950e0..e88c87ebc7 100644
--- a/src/sampling/sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.C
+++ b/src/sampling/sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.C
@@ -53,7 +53,11 @@ void Foam::triSurfaceMeshPointSet::calcSamples
 {
     forAll(sampleCoords_, sampleI)
     {
-        label cellI = searchEngine().findCell(sampleCoords_[sampleI]);
+        label cellI = searchEngine().findCell
+        (
+            sampleCoords_[sampleI],
+            polyMesh::FACEDIAGTETS
+        );
 
         if (cellI != -1)
         {
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
index a232a3eb83..4c050ff485 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
@@ -241,7 +241,7 @@ bool Foam::sampledTriSurfaceMesh::update()
     const pointField& fc = surface_.faceCentres();
 
     // Mesh search engine, no triangulation of faces.
-    meshSearch meshSearcher(mesh(), false);
+    meshSearch meshSearcher(mesh(), polyMesh::FACEPLANES);
 
 
     List<nearInfo> nearest(fc.size());
@@ -435,7 +435,15 @@ bool Foam::sampledTriSurfaceMesh::update()
                 sampleElements_[pointI] = cellI;
 
                 // Check if point inside cell
-                if (mesh().pointInCell(pt, sampleElements_[pointI]))
+                if
+                (
+                    mesh().pointInCell
+                    (
+                        pt,
+                        sampleElements_[pointI],
+                        meshSearcher.decompMode()
+                    )
+                )
                 {
                     samplePoints_[pointI] = pt;
                 }