diff --git a/applications/solvers/combustion/XiFoam/bEqn.H b/applications/solvers/combustion/XiFoam/bEqn.H
index e4a2f1ac39..df37bf34c5 100644
--- a/applications/solvers/combustion/XiFoam/bEqn.H
+++ b/applications/solvers/combustion/XiFoam/bEqn.H
@@ -63,6 +63,7 @@ if (ign.ignited())
 
     // Solve for b
     // ~~~~~~~~~~~
+    bEqn.relax();
     bEqn.solve();
 
     Info<< "min(b) = " << min(b).value() << endl;
@@ -201,7 +202,7 @@ if (ign.ignited())
         volScalarField XiEq
         (
             scalar(1.001)
-          + (scalar(1) + (2*XiShapeCoef)*(scalar(0.5) - b))
+          + (scalar(1) + (2*XiShapeCoef)*(scalar(0.5) - min(max(b, 0.0), 1.0)))
             *(XiEqStar - scalar(1.001))
         );
 
diff --git a/applications/solvers/stressAnalysis/solidDisplacementFoam/readMechanicalProperties.H b/applications/solvers/stressAnalysis/solidDisplacementFoam/readMechanicalProperties.H
index c8e3ccba56..762245abf4 100644
--- a/applications/solvers/stressAnalysis/solidDisplacementFoam/readMechanicalProperties.H
+++ b/applications/solvers/stressAnalysis/solidDisplacementFoam/readMechanicalProperties.H
@@ -13,70 +13,176 @@
     );
 
     const dictionary& rhoDict(mechanicalProperties.subDict("rho"));
-    word rhoType(rhoDict.lookup("rho"));
+    word rhoType(rhoDict.lookup("type"));
 
-    volScalarField rho
+    autoPtr<volScalarField> rhoPtr;
+
+    IOobject rhoIO
     (
-        IOobject
-        (
-            "rho",
-            runTime.timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        ),
+        "rho",
+        runTime.timeName(0),
         mesh,
-        dimensionedScalar("zero", dimMass/dimVolume, 0.0)
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
     );
 
-    if (rhoType == "rhoInf")
+    if (rhoType == "uniform")
     {
-        rho = rhoDict.lookup("rhoInf");
+        scalar rhoValue(readScalar(rhoDict.lookup("value")));
+
+        rhoPtr.reset
+        (
+            new volScalarField
+            (
+                rhoIO,
+                mesh,
+                dimensionedScalar
+                (
+                    "rho",
+                    dimMass/dimVolume,
+                    rhoValue
+                ),
+                zeroGradientFvPatchField<scalar>::typeName
+            )
+        );
+    }
+    else if (rhoType == "field")
+    {
+        rhoIO.readOpt() = IOobject::MUST_READ;
+
+        rhoPtr.reset
+        (
+            new volScalarField
+            (
+                rhoIO,
+                mesh
+            )
+        );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "readMechanicalProperties.H"
+        )   << "Valid type entries are uniform or field for rho"
+            << abort(FatalError);
     }
 
-    volScalarField rhoE
-    (
-        IOobject
-        (
-            "E",
-            runTime.timeName(0),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("0", dimMass/dimLength/sqr(dimTime), 0.0)
-    );
+    volScalarField& rho = rhoPtr();
 
     const dictionary& EDict(mechanicalProperties.subDict("E"));
-    word EType(EDict.lookup("E"));
-    if (EType == "EInf")
+    word EType(EDict.lookup("type"));
+
+    autoPtr<volScalarField> EPtr;
+
+    IOobject EIO
+    (
+        "E",
+        runTime.timeName(0),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (EType == "uniform")
     {
-        rhoE = EDict.lookup("EInf");
+        scalar rhoEValue(readScalar(EDict.lookup("value")));
+
+        EPtr.reset
+        (
+            new volScalarField
+            (
+                EIO,
+                mesh,
+                dimensionedScalar
+                (
+                    "Erho",
+                    dimMass/dimLength/sqr(dimTime),
+                    rhoEValue
+                ),
+                zeroGradientFvPatchField<scalar>::typeName
+            )
+        );
+    }
+    else if (EType == "field")
+    {
+        EIO.readOpt() = IOobject::MUST_READ;
+
+        EPtr.reset
+        (
+            new volScalarField
+            (
+                EIO,
+                mesh
+            )
+        );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "readMechanicalProperties.H"
+        )   << "Valid type entries are uniform or field for E"
+            << abort(FatalError);
     }
 
+    volScalarField& rhoE = EPtr();
 
-    volScalarField nu
+    autoPtr<volScalarField> nuPtr;
+
+    IOobject nuIO
     (
-        IOobject
-        (
-            "nu",
-            runTime.timeName(0),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        ),
+        "nu",
+        runTime.timeName(0),
         mesh,
-        dimensionedScalar("0", dimless, 0.0)
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
     );
 
     const dictionary& nuDict(mechanicalProperties.subDict("nu"));
-    word nuType(nuDict.lookup("nu"));
+    word nuType(nuDict.lookup("type"));
 
-    if (nuType == "nuInf")
+    if (nuType == "uniform")
     {
-        nu = nuDict.lookup("nuInf");
+        scalar nuValue(readScalar(nuDict.lookup("value")));
+        nuPtr.reset
+        (
+            new volScalarField
+            (
+                nuIO,
+                mesh,
+                dimensionedScalar
+                (
+                    "nu",
+                    dimless,
+                    nuValue
+                ),
+                zeroGradientFvPatchField<scalar>::typeName
+            )
+        );
     }
+    else if(nuType == "field")
+    {
+        nuIO.readOpt() = IOobject::MUST_READ;
+        nuPtr.reset
+        (
+            new volScalarField
+            (
+                nuIO,
+                mesh
+            )
+        );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "readMechanicalProperties.H"
+        )   << "Valid type entries are uniform or field for nu"
+            << abort(FatalError);
+    }
+
+    volScalarField& nu = nuPtr();
 
     Info<< "Normalising E : E/rho\n" << endl;
     volScalarField E = rhoE/rho;
diff --git a/applications/solvers/stressAnalysis/solidDisplacementFoam/readThermalProperties.H b/applications/solvers/stressAnalysis/solidDisplacementFoam/readThermalProperties.H
index 4d7bb43fc9..54ac7d6db9 100644
--- a/applications/solvers/stressAnalysis/solidDisplacementFoam/readThermalProperties.H
+++ b/applications/solvers/stressAnalysis/solidDisplacementFoam/readThermalProperties.H
@@ -46,70 +46,180 @@ volScalarField DT
 
 if (thermalStress)
 {
-    volScalarField C
+
+    autoPtr<volScalarField> CPtr;
+
+    IOobject CIO
     (
-        IOobject
-        (
-            "C",
-            runTime.timeName(0),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        ),
+        "C",
+        runTime.timeName(0),
         mesh,
-        dimensionedScalar("0", dimensionSet(0, 2, -2 , -1, 0), 0.0)
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
     );
 
     const dictionary& CDict(thermalProperties.subDict("C"));
-    word CType(CDict.lookup("C"));
-    if (CType == "CInf")
+    word CType(CDict.lookup("type"));
+    if (CType == "uniform")
     {
-        C = CDict.lookup("CInf");
+        scalar CValue(readScalar(CDict.lookup("value")));
+
+        CPtr.reset
+        (
+            new volScalarField
+            (
+                CIO,
+                mesh,
+                dimensionedScalar
+                (
+                    "C",
+                    dimensionSet(0, 2, -2 , -1, 0),
+                    CValue
+                ),
+                zeroGradientFvPatchField<scalar>::typeName
+            )
+        );
+
+    }
+    else if(CType == "field")
+    {
+        CIO.readOpt() = IOobject::MUST_READ;
+
+        CPtr.reset
+        (
+            new volScalarField
+            (
+                CIO,
+                mesh
+            )
+        );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+             "readThermalProperties.H"
+        )   << "Valid type entries are uniform or field for C"
+            << abort(FatalError);
     }
 
+    volScalarField& C = CPtr();
 
-    volScalarField rhoK
+    autoPtr<volScalarField> rhoKPtr;
+
+    IOobject rhoKIO
     (
-        IOobject
-        (
-            "k",
-            runTime.timeName(0),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        ),
+        "k",
+        runTime.timeName(0),
         mesh,
-        dimensionedScalar("0", dimensionSet(1, 1, -3 , -1, 0), 0.0)
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
     );
 
     const dictionary& kDict(thermalProperties.subDict("k"));
-    word kType(kDict.lookup("k"));
-    if (kType == "kInf")
+    word kType(kDict.lookup("type"));
+    if (kType == "uniform")
     {
-        rhoK = kDict.lookup("kInf");
+        scalar rhoKValue(readScalar(kDict.lookup("value")));
+
+        rhoKPtr.reset
+        (
+            new volScalarField
+            (
+                rhoKIO,
+                mesh,
+                dimensionedScalar
+                (
+                    "rhoK",
+                    dimensionSet(1, 1, -3 , -1, 0),
+                    rhoKValue
+                ),
+                zeroGradientFvPatchField<scalar>::typeName
+            )
+        );
+
+    }
+    else if (kType == "field")
+    {
+        rhoKIO.readOpt() = IOobject::MUST_READ;
+
+        rhoKPtr.reset
+        (
+            new volScalarField
+            (
+                rhoKIO,
+                mesh
+            )
+        );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+             "readThermalProperties.H"
+        )   << "Valid type entries are uniform or field for K"
+            << abort(FatalError);
     }
 
-    volScalarField alpha
+    volScalarField& rhoK = rhoKPtr();
+
+    autoPtr<volScalarField> alphaPtr;
+
+    IOobject alphaIO
     (
-        IOobject
-        (
-            "alpha",
-            runTime.timeName(0),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        ),
+        "alpha",
+        runTime.timeName(0),
         mesh,
-        dimensionedScalar("0", dimensionSet(0, 0, 0 , -1, 0), 0.0)
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
     );
 
-    const dictionary& alphaDict(thermalProperties.subDict("alpha"));
-    word alphaType(alphaDict.lookup("alpha"));
 
-    if (alphaType == "alphaInf")
+    const dictionary& alphaDict(thermalProperties.subDict("alpha"));
+    word alphaType(alphaDict.lookup("type"));
+
+    if (alphaType == "uniform")
     {
-        alpha = alphaDict.lookup("alphaInf");
+        scalar alphaValue(readScalar(alphaDict.lookup("value")));
+        alphaPtr.reset
+        (
+            new volScalarField
+            (
+                alphaIO,
+                mesh,
+                dimensionedScalar
+                (
+                    "alpha",
+                    inv(dimTemperature),
+                    alphaValue
+                ),
+                zeroGradientFvPatchField<scalar>::typeName
+            )
+        );
     }
+    else if (alphaType == "field")
+    {
+        alphaIO.readOpt() = IOobject::MUST_READ;
+
+        alphaPtr.reset
+        (
+            new volScalarField
+            (
+                alphaIO,
+                mesh
+            )
+        );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+             "readThermalProperties.H"
+        )   << "Valid type entries are uniform or field for alpha"
+            << abort(FatalError);
+    }
+
+    volScalarField& alpha = alphaPtr();
 
     Info<< "Normalising k : k/rho\n" << endl;
     volScalarField k = rhoK/rho;
diff --git a/applications/solvers/stressAnalysis/solidEquilibriumDisplacementFoam/tractionDisplacementCorrection/tractionDisplacementCorrectionFvPatchVectorField.C b/applications/solvers/stressAnalysis/solidEquilibriumDisplacementFoam/tractionDisplacementCorrection/tractionDisplacementCorrectionFvPatchVectorField.C
index 8aff8879e7..47a2005992 100644
--- a/applications/solvers/stressAnalysis/solidEquilibriumDisplacementFoam/tractionDisplacementCorrection/tractionDisplacementCorrectionFvPatchVectorField.C
+++ b/applications/solvers/stressAnalysis/solidEquilibriumDisplacementFoam/tractionDisplacementCorrection/tractionDisplacementCorrectionFvPatchVectorField.C
@@ -150,13 +150,18 @@ void tractionDisplacementCorrectionFvPatchVectorField::updateCoeffs()
         "mechanicalProperties"
     );
 
-    dimensionedScalar rho(mechanicalProperties.lookup("rho"));
-    dimensionedScalar rhoE(mechanicalProperties.lookup("E"));
-    dimensionedScalar nu(mechanicalProperties.lookup("nu"));
+    const fvPatchField<scalar>& rho =
+        patch().lookupPatchField<volScalarField, scalar>("rho");
 
-    dimensionedScalar E = rhoE/rho;
-    dimensionedScalar mu = E/(2.0*(1.0 + nu));
-    dimensionedScalar lambda = nu*E/((1.0 + nu)*(1.0 - 2.0*nu));
+    const fvPatchField<scalar>& rhoE =
+        patch().lookupPatchField<volScalarField, scalar>("E");
+
+    const fvPatchField<scalar>& nu =
+        patch().lookupPatchField<volScalarField, scalar>("nu");
+
+    scalarField E = rhoE/rho;
+    scalarField mu = E/(2.0*(1.0 + nu));
+    scalarField lambda = nu*E/((1.0 + nu)*(1.0 - 2.0*nu));
 
     Switch planeStress(mechanicalProperties.lookup("planeStress"));
 
@@ -175,8 +180,8 @@ void tractionDisplacementCorrectionFvPatchVectorField::updateCoeffs()
 
     gradient() =
     (
-        (traction_ + pressure_*n)/rho.value() - (n & (sigmaD + sigmaExp))
-    )/(2.0*mu + lambda).value();
+        (traction_ + pressure_*n)/rho - (n & (sigmaD + sigmaExp))
+    )/(2.0*mu + lambda);
 
     fixedGradientFvPatchVectorField::updateCoeffs();
 }
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/matrices/lduMatrix/solvers/GAMG/interfaces/cyclicGAMGInterface/cyclicGAMGInterface.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/cyclicGAMGInterface/cyclicGAMGInterface.C
index e4d56022f0..39038f5f29 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/cyclicGAMGInterface/cyclicGAMGInterface.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/cyclicGAMGInterface/cyclicGAMGInterface.C
@@ -25,7 +25,8 @@ License
 
 #include "cyclicGAMGInterface.H"
 #include "addToRunTimeSelectionTable.H"
-#include "Map.H"
+#include "labelPair.H"
+#include "HashTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -62,180 +63,65 @@ Foam::cyclicGAMGInterface::cyclicGAMGInterface
     ),
     fineCyclicInterface_(refCast<const cyclicLduInterface>(fineInterface))
 {
-    // Make a lookup table of entries for owner/neighbour
-    Map<SLList<label> > neighboursTable
+    // From coarse face to coarse cell
+    DynamicList<label> dynFaceCells(localRestrictAddressing.size());
+    // From fine face to coarse face
+    DynamicList<label> dynFaceRestrictAddressing
     (
         localRestrictAddressing.size()
     );
 
-    // Table of face-sets to be agglomerated
-    Map<SLList<SLList<label> > > faceFaceTable
+    // From coarse cell pair to coarse face
+    HashTable<label, labelPair, labelPair::Hash<> > cellsToCoarseFace
     (
-        localRestrictAddressing.size()
+        2*localRestrictAddressing.size()
     );
 
-    label nCoarseFaces = 0;
-
     forAll(localRestrictAddressing, ffi)
     {
-        label curMaster = -1;
-        label curSlave = -1;
+        labelPair cellPair;
 
         // Do switching on master/slave indexes based on the owner/neighbour of
         // the processor index such that both sides get the same answer.
         if (owner())
         {
             // Master side
-            curMaster = localRestrictAddressing[ffi];
-            curSlave = neighbourRestrictAddressing[ffi];
+            cellPair = labelPair
+            (
+                localRestrictAddressing[ffi],
+                neighbourRestrictAddressing[ffi]
+            );
         }
         else
         {
             // Slave side
-            curMaster = neighbourRestrictAddressing[ffi];
-            curSlave = localRestrictAddressing[ffi];
+            cellPair = labelPair
+            (
+                neighbourRestrictAddressing[ffi],
+                localRestrictAddressing[ffi]
+            );
         }
 
-        // Look for the master cell.  If it has already got a face,
-        // add the coefficient to the face.  If not, create a new face.
-        if (neighboursTable.found(curMaster))
+        HashTable<label, labelPair, labelPair::Hash<> >::const_iterator fnd =
+            cellsToCoarseFace.find(cellPair);
+
+        if (fnd == cellsToCoarseFace.end())
         {
-            // Check all current neighbours to see if the current slave already
-            // exists and if so, add the fine face to the agglomeration.
-
-            SLList<label>& curNbrs = neighboursTable.find(curMaster)();
-
-            SLList<SLList<label> >& curFaceFaces =
-                faceFaceTable.find(curMaster)();
-
-            bool nbrFound = false;
-
-            SLList<label>::iterator nbrsIter = curNbrs.begin();
-
-            SLList<SLList<label> >::iterator faceFacesIter =
-                curFaceFaces.begin();
-
-            for
-            (
-                ;
-                nbrsIter != curNbrs.end(), faceFacesIter != curFaceFaces.end();
-                ++nbrsIter, ++faceFacesIter
-            )
-            {
-                if (nbrsIter() == curSlave)
-                {
-                    nbrFound = true;
-                    faceFacesIter().append(ffi);
-                    break;
-                }
-            }
-
-            if (!nbrFound)
-            {
-                curNbrs.append(curSlave);
-                curFaceFaces.append(ffi);
-
-                // New coarse face created
-                nCoarseFaces++;
-            }
+            // New coarse face
+            label coarseI = dynFaceCells.size();
+            dynFaceRestrictAddressing.append(coarseI);
+            dynFaceCells.append(localRestrictAddressing[ffi]);
+            cellsToCoarseFace.insert(cellPair, coarseI);
         }
         else
         {
-            // This master has got no neighbours yet.  Add a neighbour
-            // and a coefficient, thus creating a new face
-            neighboursTable.insert(curMaster, SLList<label>(curSlave));
-            faceFaceTable.insert(curMaster, SLList<SLList<label> >(ffi));
-
-            // New coarse face created
-            nCoarseFaces++;
-        }
-    } // end for all fine faces
-
-
-
-    faceCells_.setSize(nCoarseFaces, -1);
-    faceRestrictAddressing_.setSize(localRestrictAddressing.size());
-
-    labelList contents = neighboursTable.toc();
-
-    // Reset face counter for re-use
-    nCoarseFaces = 0;
-
-    if (owner())
-    {
-        // On master side, the owner addressing is stored in table of contents
-        forAll(contents, masterI)
-        {
-            SLList<label>& curNbrs = neighboursTable.find(contents[masterI])();
-
-            SLList<SLList<label> >& curFaceFaces =
-                faceFaceTable.find(contents[masterI])();
-
-            SLList<label>::iterator nbrsIter = curNbrs.begin();
-            SLList<SLList<label> >::iterator faceFacesIter =
-                curFaceFaces.begin();
-
-            for
-            (
-                ;
-                nbrsIter != curNbrs.end(), faceFacesIter != curFaceFaces.end();
-                ++nbrsIter, ++faceFacesIter
-            )
-            {
-                faceCells_[nCoarseFaces] = contents[masterI];
-
-                for
-                (
-                    SLList<label>::iterator facesIter = faceFacesIter().begin();
-                    facesIter != faceFacesIter().end();
-                    ++facesIter
-                )
-                {
-                    faceRestrictAddressing_[facesIter()] = nCoarseFaces;
-                }
-
-                nCoarseFaces++;
-            }
+            // Already have coarse face
+            dynFaceRestrictAddressing.append(fnd());
         }
     }
-    else
-    {
-        // On slave side, the owner addressing is stored in linked lists
-        forAll(contents, masterI)
-        {
-            SLList<label>& curNbrs = neighboursTable.find(contents[masterI])();
 
-            SLList<SLList<label> >& curFaceFaces =
-                faceFaceTable.find(contents[masterI])();
-
-            SLList<label>::iterator nbrsIter = curNbrs.begin();
-
-            SLList<SLList<label> >::iterator faceFacesIter =
-                curFaceFaces.begin();
-
-            for
-            (
-                ;
-                nbrsIter != curNbrs.end(), faceFacesIter != curFaceFaces.end();
-                ++nbrsIter, ++faceFacesIter
-            )
-            {
-                faceCells_[nCoarseFaces] = nbrsIter();
-
-                for
-                (
-                    SLList<label>::iterator facesIter = faceFacesIter().begin();
-                    facesIter != faceFacesIter().end();
-                    ++facesIter
-                )
-                {
-                    faceRestrictAddressing_[facesIter()] = nCoarseFaces;
-                }
-
-                nCoarseFaces++;
-            }
-        }
-    }
+    faceCells_.transfer(dynFaceCells);
+    faceRestrictAddressing_.transfer(dynFaceRestrictAddressing);
 }
 
 
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/processorGAMGInterface/processorGAMGInterface.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/processorGAMGInterface/processorGAMGInterface.C
index 256956faff..e5e640f6f8 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/processorGAMGInterface/processorGAMGInterface.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/processorGAMGInterface/processorGAMGInterface.C
@@ -25,7 +25,8 @@ License
 
 #include "processorGAMGInterface.H"
 #include "addToRunTimeSelectionTable.H"
-#include "Map.H"
+#include "HashTable.H"
+#include "labelPair.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -62,171 +63,65 @@ Foam::processorGAMGInterface::processorGAMGInterface
     ),
     fineProcInterface_(refCast<const processorLduInterface>(fineInterface))
 {
-    // Make a lookup table of entries for owner/neighbour
-    Map<SLList<label> > neighboursTable
+    // From coarse face to coarse cell
+    DynamicList<label> dynFaceCells(localRestrictAddressing.size());
+    // From fine face to coarse face
+    DynamicList<label> dynFaceRestrictAddressing
     (
         localRestrictAddressing.size()
     );
 
-    // Table of face-sets to be agglomerated
-    Map<SLList<SLList<label> > > faceFaceTable
+    // From coarse cell pair to coarse face
+    HashTable<label, labelPair, labelPair::Hash<> > cellsToCoarseFace
     (
-        localRestrictAddressing.size()
+        2*localRestrictAddressing.size()
     );
 
-    label nCoarseFaces = 0;
-
     forAll(localRestrictAddressing, ffi)
     {
-        label curMaster = -1;
-        label curSlave = -1;
+        labelPair cellPair;
 
         // Do switching on master/slave indexes based on the owner/neighbour of
         // the processor index such that both sides get the same answer.
         if (myProcNo() < neighbProcNo())
         {
             // Master side
-            curMaster = localRestrictAddressing[ffi];
-            curSlave = neighbourRestrictAddressing[ffi];
+            cellPair = labelPair
+            (
+                localRestrictAddressing[ffi],
+                neighbourRestrictAddressing[ffi]
+            );
         }
         else
         {
             // Slave side
-            curMaster = neighbourRestrictAddressing[ffi];
-            curSlave = localRestrictAddressing[ffi];
+            cellPair = labelPair
+            (
+                neighbourRestrictAddressing[ffi],
+                localRestrictAddressing[ffi]
+            );
         }
 
-        // Look for the master cell.  If it has already got a face,
-        // add the coefficient to the face.  If not, create a new face.
-        if (neighboursTable.found(curMaster))
+        HashTable<label, labelPair, labelPair::Hash<> >::const_iterator fnd =
+            cellsToCoarseFace.find(cellPair);
+
+        if (fnd == cellsToCoarseFace.end())
         {
-            // Check all current neighbours to see if the current slave already
-            // exists and if so, add the fine face to the agglomeration.
-
-            SLList<label>& curNbrs = neighboursTable.find(curMaster)();
-
-            SLList<SLList<label> >& curFaceFaces =
-                faceFaceTable.find(curMaster)();
-
-            bool nbrFound = false;
-
-            SLList<label>::iterator nbrsIter = curNbrs.begin();
-
-            SLList<SLList<label> >::iterator faceFacesIter =
-                curFaceFaces.begin();
-
-            for
-            (
-                ;
-                nbrsIter != curNbrs.end(), faceFacesIter != curFaceFaces.end();
-                ++nbrsIter, ++faceFacesIter
-            )
-            {
-                if (nbrsIter() == curSlave)
-                {
-                    nbrFound = true;
-                    faceFacesIter().append(ffi);
-                    break;
-                }
-            }
-
-            if (!nbrFound)
-            {
-                curNbrs.append(curSlave);
-                curFaceFaces.append(ffi);
-
-                // New coarse face created
-                nCoarseFaces++;
-            }
+            // New coarse face
+            label coarseI = dynFaceCells.size();
+            dynFaceRestrictAddressing.append(coarseI);
+            dynFaceCells.append(localRestrictAddressing[ffi]);
+            cellsToCoarseFace.insert(cellPair, coarseI);
         }
         else
         {
-            // This master has got no neighbours yet.  Add a neighbour
-            // and a coefficient, thus creating a new face
-            neighboursTable.insert(curMaster, SLList<label>(curSlave));
-            faceFaceTable.insert(curMaster, SLList<SLList<label> >(ffi));
-
-            // New coarse face created
-            nCoarseFaces++;
-        }
-    } // end for all fine faces
-
-
-
-    faceCells_.setSize(nCoarseFaces, -1);
-    faceRestrictAddressing_.setSize(localRestrictAddressing.size());
-
-    labelList contents = neighboursTable.toc();
-
-    // Reset face counter for re-use
-    nCoarseFaces = 0;
-
-    if (myProcNo() < neighbProcNo())
-    {
-        // On master side, the owner addressing is stored in table of contents
-        forAll(contents, masterI)
-        {
-            SLList<label>& curNbrs = neighboursTable.find(contents[masterI])();
-
-            SLList<SLList<label> >& curFaceFaces =
-                faceFaceTable.find(contents[masterI])();
-
-            SLList<label>::iterator nbrsIter = curNbrs.begin();
-
-            SLList<SLList<label> >::iterator faceFacesIter =
-                curFaceFaces.begin();
-
-            for
-            (
-                ;
-                nbrsIter != curNbrs.end(), faceFacesIter != curFaceFaces.end();
-                ++nbrsIter, ++faceFacesIter
-            )
-            {
-                faceCells_[nCoarseFaces] = contents[masterI];
-
-                forAllConstIter(SLList<label>, faceFacesIter(), facesIter)
-                {
-                    faceRestrictAddressing_[facesIter()] = nCoarseFaces;
-                }
-
-                nCoarseFaces++;
-            }
+            // Already have coarse face
+            dynFaceRestrictAddressing.append(fnd());
         }
     }
-    else
-    {
-        // On slave side, the owner addressing is stored in linked lists
-        forAll(contents, masterI)
-        {
-            SLList<label>& curNbrs = neighboursTable.find(contents[masterI])();
 
-            SLList<SLList<label> >& curFaceFaces =
-                faceFaceTable.find(contents[masterI])();
-
-            SLList<label>::iterator nbrsIter = curNbrs.begin();
-
-            SLList<SLList<label> >::iterator faceFacesIter =
-                curFaceFaces.begin();
-
-            for
-            (
-                ;
-                nbrsIter != curNbrs.end(), faceFacesIter != curFaceFaces.end();
-                ++nbrsIter, ++faceFacesIter
-            )
-            {
-                faceCells_[nCoarseFaces] = nbrsIter();
-
-                forAllConstIter(SLList<label>, faceFacesIter(), facesIter)
-                {
-                    faceRestrictAddressing_[facesIter()] = nCoarseFaces;
-                }
-
-                nCoarseFaces++;
-            }
-        }
-    }
+    faceCells_.transfer(dynFaceCells);
+    faceRestrictAddressing_.transfer(dynFaceRestrictAddressing);
 }
 
 
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/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
index 622748be1e..63f533e090 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
@@ -680,13 +680,17 @@ void Foam::polyTopoChange::getFaceOrder
     // First unassigned face
     label newFaceI = 0;
 
+    labelList nbr;
+    labelList order;
+
     forAll(cellMap_, cellI)
     {
         label startOfCell = cellFaceOffsets[cellI];
         label nFaces = cellFaceOffsets[cellI+1] - startOfCell;
 
         // Neighbouring cells
-        SortableList<label> nbr(nFaces);
+        //SortableList<label> nbr(nFaces);
+        nbr.setSize(nFaces);
 
         for (label i = 0; i < nFaces; i++)
         {
@@ -725,14 +729,24 @@ void Foam::polyTopoChange::getFaceOrder
             }
         }
 
-        nbr.sort();
+        //nbr.sort();
+        order.setSize(nFaces);
+        sortedOrder(nbr, order);
 
-        forAll(nbr, i)
+        //forAll(nbr, i)
+        //{
+        //    if (nbr[i] != -1)
+        //    {
+        //        oldToNew[cellFaces[startOfCell + nbr.indices()[i]]] =
+        //            newFaceI++;
+        //    }
+        //}
+        forAll(order, i)
         {
-            if (nbr[i] != -1)
+            label index = order[i];
+            if (nbr[index] != -1)
             {
-                oldToNew[cellFaces[startOfCell + nbr.indices()[i]]] =
-                    newFaceI++;
+                oldToNew[cellFaces[startOfCell + index]] = newFaceI++;
             }
         }
     }
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/postProcessing/functionObjects/utilities/codedFunctionObject/codedFunctionObject.C b/src/postProcessing/functionObjects/utilities/codedFunctionObject/codedFunctionObject.C
index 45a915456a..7745c7e308 100644
--- a/src/postProcessing/functionObjects/utilities/codedFunctionObject/codedFunctionObject.C
+++ b/src/postProcessing/functionObjects/utilities/codedFunctionObject/codedFunctionObject.C
@@ -129,7 +129,8 @@ Foam::codedFunctionObject::codedFunctionObject
 :
     functionObject(name),
     codedBase(),
-    time_(time)
+    time_(time),
+    dict_(dict)
 {
     if (readNow)
     {
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;
                 }
diff --git a/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H b/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H
index 9154126ada..bd8bdb4573 100644
--- a/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H
+++ b/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H
@@ -135,6 +135,9 @@ public:
 
     // Member Functions
 
+        //- Limit the temperature to be in the range Tlow_ to Thigh_
+        inline scalar limit(const scalar T) const;
+
         // Fundamental properties
 
             //- Heat capacity at constant pressure [J/(kmol K)]
diff --git a/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H b/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H
index f031aa63f6..31bd392b5d 100644
--- a/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H
@@ -89,6 +89,16 @@ Foam::eConstThermo<EquationOfState>::New(const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class EquationOfState>
+inline Foam::scalar Foam::eConstThermo<EquationOfState>::limit
+(
+    const scalar T
+) const
+{
+    return T;
+}
+
+
 template<class EquationOfState>
 inline Foam::scalar Foam::eConstThermo<EquationOfState>::cp
 (
diff --git a/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H b/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H
index 014de5eef2..9a8352a064 100644
--- a/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H
+++ b/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H
@@ -133,6 +133,9 @@ public:
 
     // Member Functions
 
+        //- Limit the temperature to be in the range Tlow_ to Thigh_
+        inline scalar limit(const scalar T) const;
+
         // Fundamental properties
 
             //- Heat capacity at constant pressure [J/(kmol K)]
diff --git a/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H b/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H
index 76ec41b620..0a4a35358c 100644
--- a/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H
@@ -89,6 +89,16 @@ Foam::hConstThermo<equationOfState>::New(const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class EquationOfState>
+inline Foam::scalar Foam::hConstThermo<EquationOfState>::limit
+(
+    const scalar T
+) const
+{
+    return T;
+}
+
+
 template<class equationOfState>
 inline Foam::scalar Foam::hConstThermo<equationOfState>::cp
 (
diff --git a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H
index 7a4dcb3867..31954a4a56 100644
--- a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H
+++ b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H
@@ -151,6 +151,9 @@ public:
 
     // Member Functions
 
+        //- Limit the temperature to be in the range Tlow_ to Thigh_
+        inline scalar limit(const scalar T) const;
+
         // Fundamental properties
 
             //- Heat capacity at constant pressure [J/(kmol K)]
diff --git a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
index 620019e6d1..a164392c61 100644
--- a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
@@ -82,6 +82,16 @@ inline Foam::hPolynomialThermo<EquationOfState, PolySize>::hPolynomialThermo
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class EquationOfState, int PolySize>
+inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::limit
+(
+    const scalar T
+) const
+{
+    return T;
+}
+
+
 template<class EquationOfState, int PolySize>
 inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::cp
 (
diff --git a/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H b/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H
index 9c414afaac..b552fad094 100644
--- a/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H
+++ b/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H
@@ -119,9 +119,6 @@ private:
         //- Check that input data is valid
         void checkInputData() const;
 
-        //- Check given temperature is within the range of the fitted coeffs
-        inline void checkT(const scalar T) const;
-
         //- Return the coefficients corresponding to the given temperature
         inline const coeffArray& coeffs(const scalar T) const;
 
@@ -153,6 +150,9 @@ public:
 
     // Member Functions
 
+        //- Limit the temperature to be in the range Tlow_ to Thigh_
+        inline scalar limit(const scalar T) const;
+
         // Fundamental properties
 
             //- Heat capacity at constant pressure [J/(kmol K)]
diff --git a/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H b/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
index 05f0c0c0cc..2135d53e03 100644
--- a/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
@@ -52,22 +52,6 @@ inline Foam::janafThermo<EquationOfState>::janafThermo
 }
 
 
-template<class EquationOfState>
-inline void Foam::janafThermo<EquationOfState>::checkT(const scalar T) const
-{
-    if (T < Tlow_ || T > Thigh_)
-    {
-        FatalErrorIn
-        (
-            "janafThermo<EquationOfState>::checkT(const scalar T) const"
-        )   << "attempt to use janafThermo<EquationOfState>"
-               " out of temperature range "
-            << Tlow_ << " -> " << Thigh_ << ";  T = " << T
-            << abort(FatalError);
-    }
-}
-
-
 template<class EquationOfState>
 inline const typename Foam::janafThermo<EquationOfState>::coeffArray&
 Foam::janafThermo<EquationOfState>::coeffs
@@ -75,8 +59,6 @@ Foam::janafThermo<EquationOfState>::coeffs
     const scalar T
 ) const
 {
-    checkT(T);
-
     if (T < Tcommon_)
     {
         return lowCpCoeffs_;
@@ -112,6 +94,31 @@ inline Foam::janafThermo<EquationOfState>::janafThermo
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class EquationOfState>
+inline Foam::scalar Foam::janafThermo<EquationOfState>::limit
+(
+    const scalar T
+) const
+{
+    if (T < Tlow_ || T > Thigh_)
+    {
+        WarningIn
+        (
+            "janafThermo<EquationOfState>::limit(const scalar T) const"
+        )   << "attempt to use janafThermo<EquationOfState>"
+               " out of temperature range "
+            << Tlow_ << " -> " << Thigh_ << ";  T = " << T
+            << endl;
+
+        return min(max(T, Tlow_), Thigh_);
+    }
+    else
+    {
+        return T;
+    }
+}
+
+
 template<class EquationOfState>
 inline Foam::scalar Foam::janafThermo<EquationOfState>::cp
 (
@@ -196,7 +203,19 @@ inline void Foam::janafThermo<EquationOfState>::operator+=
 
     Tlow_ = max(Tlow_, jt.Tlow_);
     Thigh_ = min(Thigh_, jt.Thigh_);
-    Tcommon_ = molr1*Tcommon_ + molr2*jt.Tcommon_;
+
+    if (notEqual(Tcommon_, jt.Tcommon_))
+    {
+        FatalErrorIn
+        (
+            "janafThermo<EquationOfState>::operator+="
+            "(const janafThermo<EquationOfState>& jt) const"
+        )   << "Tcommon " << Tcommon_ << " for "
+            << (this->name().size() ? this->name() : "others")
+            << " != " << jt.Tcommon_ << " for "
+            << (jt.name().size() ? jt.name() : "others")
+            << exit(FatalError);
+    }
 
     for
     (
@@ -231,7 +250,19 @@ inline void Foam::janafThermo<EquationOfState>::operator-=
 
     Tlow_ = max(Tlow_, jt.Tlow_);
     Thigh_ = min(Thigh_, jt.Thigh_);
-    Tcommon_ = molr1*Tcommon_ - molr2*jt.Tcommon_;
+
+    if (notEqual(Tcommon_, jt.Tcommon_))
+    {
+        FatalErrorIn
+        (
+            "janafThermo<EquationOfState>::operator-="
+            "(const janafThermo<EquationOfState>& jt) const"
+        )   << "Tcommon " << Tcommon_ << " for "
+            << (this->name().size() ? this->name() : "others")
+            << " != " << jt.Tcommon_ << " for "
+            << (jt.name().size() ? jt.name() : "others")
+            << exit(FatalError);
+    }
 
     for
     (
@@ -285,12 +316,26 @@ inline Foam::janafThermo<EquationOfState> Foam::operator+
           + molr2*jt2.lowCpCoeffs_[coefLabel];
     }
 
+    if (notEqual(jt1.Tcommon_, jt2.Tcommon_))
+    {
+        FatalErrorIn
+        (
+            "operator+"
+            "(const janafThermo<EquationOfState>& jt1,"
+            " const janafThermo<EquationOfState>& jt2)"
+        )   << "Tcommon " << jt1.Tcommon_ << " for "
+            << (jt1.name().size() ? jt1.name() : "others")
+            << " != " << jt2.Tcommon_ << " for "
+            << (jt2.name().size() ? jt2.name() : "others")
+            << exit(FatalError);
+    }
+
     return janafThermo<EquationOfState>
     (
         eofs,
         max(jt1.Tlow_, jt2.Tlow_),
         min(jt1.Thigh_, jt2.Thigh_),
-        molr1*jt1.Tcommon_ + molr2*jt2.Tcommon_,
+        jt1.Tcommon_,
         highCpCoeffs,
         lowCpCoeffs
     );
@@ -329,12 +374,26 @@ inline Foam::janafThermo<EquationOfState> Foam::operator-
           - molr2*jt2.lowCpCoeffs_[coefLabel];
     }
 
+    if (notEqual(jt1.Tcommon_, jt2.Tcommon_))
+    {
+        FatalErrorIn
+        (
+            "operator-"
+            "(const janafThermo<EquationOfState>& jt1,"
+            " const janafThermo<EquationOfState>& jt2)"
+        )   << "Tcommon " << jt1.Tcommon_ << " for "
+            << (jt1.name().size() ? jt1.name() : "others")
+            << " != " << jt2.Tcommon_ << " for "
+            << (jt2.name().size() ? jt2.name() : "others")
+            << exit(FatalError);
+    }
+
     return janafThermo<EquationOfState>
     (
         eofs,
         max(jt1.Tlow_, jt2.Tlow_),
         min(jt1.Thigh_, jt2.Thigh_),
-        molr1*jt1.Tcommon_ - molr2*jt2.Tcommon_,
+        jt1.Tcommon_,
         highCpCoeffs,
         lowCpCoeffs
     );
diff --git a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.H b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.H
index c6f0f7a316..93268b9e6b 100644
--- a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.H
+++ b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.H
@@ -103,14 +103,15 @@ class specieThermo
 
     // Private Member Functions
 
-        //- return the temperature corresponding to the value of the
+        //- Return the temperature corresponding to the value of the
         //  thermodynamic property f, given the function f = F(T) and dF(T)/dT
         inline scalar T
         (
             scalar f,
             scalar T0,
             scalar (specieThermo::*F)(const scalar) const,
-            scalar (specieThermo::*dFdT)(const scalar) const
+            scalar (specieThermo::*dFdT)(const scalar) const,
+            scalar (specieThermo::*limit)(const scalar) const
         ) const;
 
 
diff --git a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H
index f8632376c2..b2d2edd5c1 100644
--- a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H
@@ -43,7 +43,8 @@ inline Foam::scalar Foam::specieThermo<Thermo>::T
     scalar f,
     scalar T0,
     scalar (specieThermo<Thermo>::*F)(const scalar) const,
-    scalar (specieThermo<Thermo>::*dFdT)(const scalar) const
+    scalar (specieThermo<Thermo>::*dFdT)(const scalar) const,
+    scalar (specieThermo<Thermo>::*limit)(const scalar) const
 ) const
 {
     scalar Test = T0;
@@ -54,7 +55,8 @@ inline Foam::scalar Foam::specieThermo<Thermo>::T
     do
     {
         Test = Tnew;
-        Tnew = Test - ((this->*F)(Test) - f)/(this->*dFdT)(Test);
+        Tnew =
+            (this->*limit)(Test - ((this->*F)(Test) - f)/(this->*dFdT)(Test));
 
         if (iter++ > maxIter_)
         {
@@ -276,7 +278,14 @@ inline Foam::scalar Foam::specieThermo<Thermo>::TH
     const scalar T0
 ) const
 {
-    return T(h, T0, &specieThermo<Thermo>::H, &specieThermo<Thermo>::Cp);
+    return T
+    (
+        h,
+        T0,
+        &specieThermo<Thermo>::H,
+        &specieThermo<Thermo>::Cp,
+        &specieThermo<Thermo>::limit
+    );
 }
 
 
@@ -287,7 +296,14 @@ inline Foam::scalar Foam::specieThermo<Thermo>::THs
     const scalar T0
 ) const
 {
-    return T(hs, T0, &specieThermo<Thermo>::Hs, &specieThermo<Thermo>::Cp);
+    return T
+    (
+        hs,
+        T0,
+        &specieThermo<Thermo>::Hs,
+        &specieThermo<Thermo>::Cp,
+        &specieThermo<Thermo>::limit
+    );
 }
 
 
@@ -298,7 +314,14 @@ inline Foam::scalar Foam::specieThermo<Thermo>::TE
     const scalar T0
 ) const
 {
-    return T(e, T0, &specieThermo<Thermo>::E, &specieThermo<Thermo>::Cv);
+    return T
+    (
+        e,
+        T0,
+        &specieThermo<Thermo>::E,
+        &specieThermo<Thermo>::Cv,
+        &specieThermo<Thermo>::limit
+    );
 }
 
 
diff --git a/src/turbulenceModels/derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchField.C b/src/turbulenceModels/derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchField.C
index 0d7ead98d1..07cc0c9481 100644
--- a/src/turbulenceModels/derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchField.C
+++ b/src/turbulenceModels/derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchField.C
@@ -70,17 +70,10 @@ Foam::porousBafflePressureFvPatchField<Type>::porousBafflePressureFvPatchField
     I_(readScalar(dict.lookup("I"))),
     length_(readScalar(dict.lookup("length")))
 {
-    if (dict.found("value"))
-    {
-        fvPatchField<Type>::operator=
-        (
-            Field<Type>("value", dict, p.size())
-        );
-    }
-    else
-    {
-        this->evaluate(Pstream::blocking);
-    }
+    fvPatchField<Type>::operator=
+    (
+        Field<Type>("value", dict, p.size())
+    );
 }
 
 
diff --git a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
index 5458ad5e64..eaf09f03c1 100644
--- a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
+++ b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
@@ -426,6 +426,7 @@ void LienCubicKELowRe::correct()
 
     volScalarField G
     (
+        "RASModel::G",
         Cmu_*fMu*sqr(k_)/epsilon_*S2 - (nonlinearStress_ && gradU_)
     );
 
diff --git a/src/turbulenceModels/incompressible/RAS/Make/files b/src/turbulenceModels/incompressible/RAS/Make/files
index f18b124e06..e28a685903 100644
--- a/src/turbulenceModels/incompressible/RAS/Make/files
+++ b/src/turbulenceModels/incompressible/RAS/Make/files
@@ -16,6 +16,7 @@ LienCubicKELowRe/LienCubicKELowRe.C
 NonlinearKEShih/NonlinearKEShih.C
 LienLeschzinerLowRe/LienLeschzinerLowRe.C
 LamBremhorstKE/LamBremhorstKE.C
+kkLOmega/kkLOmega.C
 
 /* Wall functions */
 wallFunctions = derivedFvPatchFields/wallFunctions
diff --git a/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C
new file mode 100644
index 0000000000..ad864bb3c7
--- /dev/null
+++ b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C
@@ -0,0 +1,779 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "kkLOmega.H"
+#include "addToRunTimeSelectionTable.H"
+
+#include "backwardsCompatibilityWallFunctions.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(kkLOmega, 0);
+addToRunTimeSelectionTable(RASModel, kkLOmega, dictionary);
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+tmp<volScalarField> kkLOmega::fv(const volScalarField& Ret) const
+{
+    return(1.0 - exp(-sqrt(Ret)/Av_));
+}
+
+
+tmp<volScalarField> kkLOmega::fINT() const
+{
+    return
+    (
+        min
+        (
+            kl_/(Cint_*(kl_ + kt_ + kMin_)),
+            dimensionedScalar("1.0", dimless, 1.0)
+        )
+    );
+}
+
+
+tmp<volScalarField> kkLOmega::fSS(const volScalarField& omega) const
+{
+    return(exp(-sqr(Css_*nu()*omega/(kt_ + kMin_))));
+}
+
+
+tmp<volScalarField> kkLOmega::Cmu(const volScalarField& S) const
+{
+    return(1.0/(A0_ + As_*(S/(omega_ + omegaMin_))));
+}
+
+
+tmp<volScalarField> kkLOmega::BetaTS(const volScalarField& Rew) const
+{
+    return(scalar(1.0) - exp(-sqr(max(Rew - CtsCrit_, 0.0))/Ats_));
+}
+
+
+tmp<volScalarField> kkLOmega::fTaul
+(
+    const volScalarField& lambdaEff,
+    const volScalarField& ktL
+) const
+{
+    return
+    (
+        scalar(1.0)
+      - exp
+        (
+            -CtauL_*ktL
+          /
+            (
+                sqr
+                (
+                    lambdaEff*omega_
+                  + dimensionedScalar
+                    (
+                        "ROTVSMALL",
+                        dimLength*inv(dimTime),
+                        ROOTVSMALL
+                    )
+                )
+            )
+        )
+    );
+}
+
+
+tmp<volScalarField> kkLOmega::alphaT
+(
+    const volScalarField& lambdaEff,
+    const volScalarField& fv,
+    const volScalarField& ktS
+) const
+{
+    return(fv*CmuStd_*sqrt(ktS)*lambdaEff);
+}
+
+
+tmp<volScalarField> kkLOmega::fOmega
+(
+    const volScalarField& lambdaEff,
+    const volScalarField& lambdaT
+) const
+{
+    return
+    (
+        scalar(1.0)
+      - exp
+        (
+            -0.41
+           * pow4
+            (
+                lambdaEff
+              / (
+                    lambdaT
+                  + dimensionedScalar
+                    (
+                        "ROTVSMALL",
+                        lambdaT.dimensions(),
+                        ROOTVSMALL
+                    )
+                )
+            )
+        )
+    );
+}
+
+
+tmp<volScalarField> kkLOmega::gammaBP(const volScalarField& omega) const
+{
+    return
+    (
+        max
+        (
+            kt_/nu()
+          /
+            (
+                omega
+              + dimensionedScalar("ROTVSMALL", omega.dimensions(), ROOTVSMALL)
+            )
+          -
+            CbpCrit_
+        ,
+            0.0
+        )
+    );
+}
+
+
+tmp<volScalarField> kkLOmega::gammaNAT
+(
+    const volScalarField& ReOmega,
+    const volScalarField& fNatCrit
+) const
+{
+    return
+    (
+        max
+        (
+            ReOmega
+          - CnatCrit_
+          / (
+                fNatCrit + dimensionedScalar("ROTVSMALL", dimless, ROOTVSMALL)
+            )
+          ,
+            0.0
+        )
+    );
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+kkLOmega::kkLOmega
+(
+    const volVectorField& U,
+    const surfaceScalarField& phi,
+    transportModel& transport,
+    const word& turbulenceModelName,
+    const word& modelName
+)
+:
+    RASModel(modelName, U, phi, transport, turbulenceModelName),
+
+    A0_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "A0",
+            coeffDict_,
+            4.04
+        )
+    ),
+    As_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "As",
+            coeffDict_,
+            2.12
+        )
+    ),
+    Av_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Av",
+            coeffDict_,
+            6.75
+        )
+    ),
+    Abp_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Abp",
+            coeffDict_,
+            0.6
+        )
+    ),
+    Anat_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Anat",
+            coeffDict_,
+            200
+        )
+    ),
+    Ats_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Ats",
+            coeffDict_,
+            200
+        )
+    ),
+    CbpCrit_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CbpCrit",
+            coeffDict_,
+            1.2
+        )
+    ),
+    Cnc_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cnc",
+            coeffDict_,
+            0.1
+        )
+    ),
+    CnatCrit_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CnatCrit",
+            coeffDict_,
+            1250
+        )
+    ),
+    Cint_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cint",
+            coeffDict_,
+            0.75
+        )
+    ),
+    CtsCrit_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CtsCrit",
+            coeffDict_,
+            1000
+        )
+    ),
+    CrNat_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CrNat",
+            coeffDict_,
+            0.02
+        )
+    ),
+    C11_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "C11",
+            coeffDict_,
+            3.4e-6
+        )
+    ),
+    C12_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "C12",
+            coeffDict_,
+            1.0e-10
+        )
+    ),
+    CR_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CR",
+            coeffDict_,
+            0.12
+        )
+    ),
+    CalphaTheta_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CalphaTheta",
+            coeffDict_,
+            0.035
+        )
+    ),
+    Css_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Css",
+            coeffDict_,
+            1.5
+        )
+    ),
+    CtauL_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CtauL",
+            coeffDict_,
+            4360
+        )
+    ),
+    Cw1_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cw1",
+            coeffDict_,
+            0.44
+        )
+    ),
+    Cw2_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cw2",
+            coeffDict_,
+            0.92
+        )
+    ),
+    Cw3_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cw3",
+            coeffDict_,
+            0.3
+        )
+    ),
+    CwR_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CwR",
+            coeffDict_,
+            1.5
+        )
+    ),
+    Clambda_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Clambda",
+            coeffDict_,
+            2.495
+        )
+    ),
+    CmuStd_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CmuStd",
+            coeffDict_,
+            0.09
+        )
+    ),
+    Prtheta_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Prtheta",
+            coeffDict_,
+            0.85
+        )
+    ),
+    Sigmak_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Sigmak",
+            coeffDict_,
+            1
+        )
+    ),
+    Sigmaw_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Sigmaw",
+            coeffDict_,
+            1.17
+        )
+    ),
+    kt_
+    (
+        IOobject
+        (
+            "kt",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateK("kt", mesh_)
+    ),
+    omega_
+    (
+        IOobject
+        (
+            "omega",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateOmega("omega", mesh_)
+    ),
+    kl_
+    (
+        IOobject
+        (
+            "kl",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateK("kl", mesh_)
+    ),
+    nut_
+    (
+        IOobject
+        (
+            "nut",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateNut("nut", mesh_)
+    ),
+    y_(mesh_)
+{
+    bound(kt_, kMin_);
+    bound(kl_, kMin_);
+    bound(omega_, omegaMin_);
+
+    nut_ = kt_/(omega_ + omegaMin_);
+    nut_.correctBoundaryConditions();
+
+    printCoeffs();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+tmp<volSymmTensorField> kkLOmega::R() const
+{
+    return tmp<volSymmTensorField>
+    (
+        new volSymmTensorField
+        (
+            IOobject
+            (
+                "R",
+                runTime_.timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            ((2.0/3.0)*I)*(kt_) - nut_*twoSymm(fvc::grad(U_)),
+            kt_.boundaryField().types()
+        )
+    );
+}
+
+
+tmp<volSymmTensorField> kkLOmega::devReff() const
+{
+    return tmp<volSymmTensorField>
+    (
+        new volSymmTensorField
+        (
+            IOobject
+            (
+                "devRhoReff",
+                runTime_.timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+           -nuEff()*dev(twoSymm(fvc::grad(U_)))
+        )
+    );
+}
+
+
+tmp<fvVectorMatrix> kkLOmega::divDevReff(volVectorField& U) const
+{
+    return
+    (
+      - fvm::laplacian(nuEff(), U)
+      - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+    );
+}
+
+
+bool kkLOmega::read()
+{
+    if (RASModel::read())
+    {
+        A0_.readIfPresent(coeffDict());
+        As_.readIfPresent(coeffDict());
+        Av_.readIfPresent(coeffDict());
+        Abp_.readIfPresent(coeffDict());
+        Anat_.readIfPresent(coeffDict());
+        Abp_.readIfPresent(coeffDict());
+        Ats_.readIfPresent(coeffDict());
+        CbpCrit_.readIfPresent(coeffDict());
+        Cnc_.readIfPresent(coeffDict());
+        CnatCrit_.readIfPresent(coeffDict());
+        Cint_.readIfPresent(coeffDict());
+        CtsCrit_.readIfPresent(coeffDict());
+        CrNat_.readIfPresent(coeffDict());
+        C11_.readIfPresent(coeffDict());
+        C12_.readIfPresent(coeffDict());
+        CR_.readIfPresent(coeffDict());
+        CalphaTheta_.readIfPresent(coeffDict());
+        Css_.readIfPresent(coeffDict());
+        CtauL_.readIfPresent(coeffDict());
+        Cw1_.readIfPresent(coeffDict());
+        Cw2_.readIfPresent(coeffDict());
+        Cw3_.readIfPresent(coeffDict());
+        CwR_.readIfPresent(coeffDict());
+        Clambda_.readIfPresent(coeffDict());
+        CmuStd_.readIfPresent(coeffDict());
+        Prtheta_.readIfPresent(coeffDict());
+        Sigmak_.readIfPresent(coeffDict());
+        Sigmaw_.readIfPresent(coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+void kkLOmega::correct()
+{
+    RASModel::correct();
+
+    if (!turbulence_)
+    {
+        return;
+    }
+
+    if (mesh_.changing())
+    {
+        y_.correct();
+        y_.boundaryField() = max(y_.boundaryField(), VSMALL);
+    }
+
+
+    const volScalarField kT(kt_ + kl_);
+
+    const volScalarField lambdaT(sqrt(kT)/(omega_ + omegaMin_));
+
+    const volScalarField lambdaEff(min(Clambda_*y_, lambdaT));
+
+    const volScalarField fw
+    (
+        lambdaEff/(lambdaT + dimensionedScalar("SMALL", dimLength, ROOTVSMALL))
+    );
+
+    const volTensorField gradU(fvc::grad(U_));
+
+    const volScalarField omega(sqrt(2.0)*mag(skew(gradU)));
+
+    const volScalarField S2(2.0*magSqr(symm(gradU)));
+
+    const volScalarField ktS(fSS(omega)*fw*kt_);
+
+    const volScalarField nuts
+    (
+        fv(sqr(fw)*kt_/nu()/(omega_ + omegaMin_))
+       *fINT()
+       *Cmu(sqrt(S2))*sqrt(ktS)*lambdaEff
+    );
+    const volScalarField Pkt(nuts*S2);
+
+    const volScalarField ktL(kt_ - ktS);
+    const volScalarField ReOmega(sqr(y_)*omega/nu());
+    const volScalarField nutl
+    (
+        min
+        (
+            C11_*fTaul(lambdaEff, ktL)*omega*sqr(lambdaEff)
+          * sqrt(ktL)*lambdaEff/nu()
+          + C12_*BetaTS(ReOmega)*ReOmega*sqr(y_)*omega
+        ,
+            0.5*(kl_ + ktL)/sqrt(S2)
+        )
+    );
+
+    const volScalarField Pkl(nutl*S2);
+
+    const volScalarField alphaTEff
+    (
+        alphaT(lambdaEff, fv(sqr(fw)*kt_/nu()/(omega_ + omegaMin_)), ktS)
+    );
+
+    // By pass s0urce term divided by kl_
+
+    const dimensionedScalar fwMin("SMALL", dimless, ROOTVSMALL);
+
+    const volScalarField Rbp
+    (
+        CR_*(1.0 - exp(-gammaBP(omega)()/Abp_))*omega_
+      / (fw + fwMin)
+    );
+
+    const volScalarField fNatCrit(1.0 - exp(-Cnc_*sqrt(kl_)*y_/nu()));
+    // Natural source term divided by kl_
+    const volScalarField Rnat
+    (
+        CrNat_*(1.0 - exp(-gammaNAT(ReOmega, fNatCrit)/Anat_))*omega
+    );
+
+    const volScalarField Dt(nu()*magSqr(fvc::grad(sqrt(kt_))));
+
+    // Turbulent kinetic energy equation
+    tmp<fvScalarMatrix> ktEqn
+    (
+        fvm::ddt(kt_)
+      + fvm::div(phi_, kt_)
+      - fvm::Sp(fvc::div(phi_), kt_)
+      - fvm::laplacian(DkEff(alphaTEff), kt_, "laplacian(alphaTEff,kt)")
+     ==
+        Pkt
+      + (Rbp + Rnat)*kl_
+      - Dt
+      - fvm::Sp(omega_, kt_)
+    );
+
+    ktEqn().relax();
+    ktEqn().boundaryManipulate(kt_.boundaryField());
+
+    solve(ktEqn);
+    bound(kt_, kMin_);
+
+
+    const volScalarField Dl(nu()*magSqr(fvc::grad(sqrt(kl_))));
+
+    // Laminar kinetic energy equation
+    tmp<fvScalarMatrix> klEqn
+    (
+        fvm::ddt(kl_)
+      + fvm::div(phi_, kl_)
+      - fvm::Sp(fvc::div(phi_), kl_)
+      - fvm::laplacian(nu(), kl_, "laplacian(nu,kl)")
+     ==
+        Pkl
+      - fvm::Sp(Rbp, kl_)
+      - fvm::Sp(Rnat, kl_)
+      - Dl
+    );
+
+    klEqn().relax();
+    klEqn().boundaryManipulate(kl_.boundaryField());
+
+    solve(klEqn);
+    bound(kl_, kMin_);
+
+
+    omega_.boundaryField().updateCoeffs();
+    // Turbulence specific dissipation rate equation
+    tmp<fvScalarMatrix> omegaEqn
+    (
+        fvm::ddt(omega_)
+      + fvm::div(phi_, omega_)
+      - fvm::Sp(fvc::div(phi_), omega_)
+      - fvm::laplacian
+        (
+            DomegaEff(alphaTEff),
+            omega_,
+            "laplacian(alphaTEff,omega)"
+        )
+     ==
+        Cw1_*Pkt*omega_/(kt_ + kMin_)
+      + fvm::SuSp
+        (
+            (CwR_/(fw + fwMin) - 1.0)*kl_*(Rbp + Rnat)/(kt_ + kMin_)
+          , omega_
+        )
+      - fvm::Sp(Cw2_*omega_, omega_)
+      + Cw3_*fOmega(lambdaEff, lambdaT)*alphaTEff*sqr(fw)*sqrt(kt_)/pow3(y_)
+    );
+
+
+    omegaEqn().relax();
+    omegaEqn().boundaryManipulate(omega_.boundaryField());
+
+    solve(omegaEqn);
+    bound(omega_, omegaMin_);
+
+    // Re-calculate viscosity
+    nut_ = nuts + nutl;
+    nut_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.H b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.H
new file mode 100644
index 0000000000..b7d451b315
--- /dev/null
+++ b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.H
@@ -0,0 +1,296 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::incompressible::RASModels::kkLOmega
+
+Description
+    Low Reynolds-number k-kl-omega turbulence model for
+    incompressible flows.
+
+    Turbulence model described in:
+    \verbatim
+        D. Keith Walters, Davor Cokljat
+        "A Three-Equation Eddy-Viscosity Model for Reynold-Averaged
+        Navier-Stokes Simulations of Transitional Flow"
+    \endverbatim
+
+    The default model coefficients correspond to the following:
+    \verbatim
+        kkLOmegaCoeffs
+        {
+            A0             4.04
+            As             2.12
+            Av             6.75
+            Abp            0.6
+            Anat           200
+            Ats            200
+            CbpCrit        1.2
+            Cnc            0.1
+            CnatCrit       1250
+            Cint           0.75
+            CtsCrit        1000
+            CrNat          0.02
+            C11            3.4e-6
+            C12            1.0e-10
+            CR             0.12
+            CalphaTheta    0.035
+            Css            1.5
+            CtauL          4360
+            Cw1            0.44
+            Cw2            0.92
+            Cw3            0.3
+            CwR            1.5
+            Clambda        2.495
+            CmuStd         0.09
+            Prtheta        0.85
+            Sigmak         1
+            Sigmaw         1.17
+        }
+    \endverbatim
+
+SourceFiles
+    kkLOmega.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kkLOmega_H
+#define kkLOmega_H
+
+#include "RASModel.H"
+#include "wallDist.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class kkLOmega Declaration
+\*---------------------------------------------------------------------------*/
+
+class kkLOmega
+:
+    public RASModel
+{
+    // Private memmber functions
+
+        tmp<volScalarField> fv(const volScalarField& Ret) const;
+
+        tmp<volScalarField> fINT() const;
+
+        tmp<volScalarField> fSS(const volScalarField& omega) const;
+
+        tmp<volScalarField> Cmu(const volScalarField& S) const;
+
+        tmp<volScalarField> BetaTS(const volScalarField& Rew) const;
+
+        tmp<volScalarField> fTaul
+        (
+            const volScalarField& lambdaEff,
+            const volScalarField& ktL
+        ) const;
+
+        tmp<volScalarField> alphaT
+        (
+            const volScalarField& lambdaEff,
+            const volScalarField& fv,
+            const volScalarField& ktS
+        ) const;
+
+        tmp<volScalarField> fOmega
+        (
+            const volScalarField& lambdaEff,
+            const volScalarField& lambdaT
+        ) const;
+
+        tmp<volScalarField> gammaBP(const volScalarField& omega) const;
+
+        tmp<volScalarField> gammaNAT
+        (
+            const volScalarField& ReOmega,
+            const volScalarField& fNatCrit
+        ) const;
+
+protected:
+
+    // Protected data
+
+        // Model coefficients
+
+            dimensionedScalar A0_;
+            dimensionedScalar As_;
+            dimensionedScalar Av_;
+            dimensionedScalar Abp_;
+            dimensionedScalar Anat_;
+            dimensionedScalar Ats_;
+            dimensionedScalar CbpCrit_;
+            dimensionedScalar Cnc_;
+            dimensionedScalar CnatCrit_;
+            dimensionedScalar Cint_;
+            dimensionedScalar CtsCrit_;
+            dimensionedScalar CrNat_;
+            dimensionedScalar C11_;
+            dimensionedScalar C12_;
+            dimensionedScalar CR_;
+            dimensionedScalar CalphaTheta_;
+            dimensionedScalar Css_;
+            dimensionedScalar CtauL_;
+            dimensionedScalar Cw1_;
+            dimensionedScalar Cw2_;
+            dimensionedScalar Cw3_;
+            dimensionedScalar CwR_;
+            dimensionedScalar Clambda_;
+            dimensionedScalar CmuStd_;
+            dimensionedScalar Prtheta_;
+            dimensionedScalar Sigmak_;
+            dimensionedScalar Sigmaw_;
+
+
+        // Fields
+
+            volScalarField kt_;
+            volScalarField omega_;
+            volScalarField kl_;
+            volScalarField nut_;
+            wallDist y_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("kkLOmega");
+
+    // Constructors
+
+        //- Construct from components
+        kkLOmega
+        (
+            const volVectorField& U,
+            const surfaceScalarField& phi,
+            transportModel& transport,
+            const word& turbulenceModelName = turbulenceModel::typeName,
+            const word& modelName = typeName
+        );
+
+
+    //- Destructor
+    virtual ~kkLOmega()
+    {}
+
+
+    // Member Functions
+
+        //- Return the turbulence viscosity
+        virtual tmp<volScalarField> nut() const
+        {
+            return nut_;
+        }
+
+        //- Return the effective diffusivity for k
+        tmp<volScalarField> DkEff(const volScalarField& alphaT) const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField("DkEff", alphaT/Sigmak_ + nu())
+            );
+        }
+
+        //- Return the effective diffusivity for omega
+        tmp<volScalarField> DomegaEff(const volScalarField& alphaT) const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField("DomegaEff", alphaT/Sigmaw_ + nu())
+            );
+        }
+
+        //- Return the laminar kinetic energy
+        virtual tmp<volScalarField> kl() const
+        {
+            return kl_;
+        }
+
+        //- Return the turbulence kinetic energy
+        virtual tmp<volScalarField> k() const
+        {
+            return kt_;
+        }
+
+        //- Return the turbulence specific dissipation rate
+        virtual tmp<volScalarField> omega() const
+        {
+            return omega_;
+        }
+
+        //- Return the turbulence kinetic energy dissipation rate
+        virtual tmp<volScalarField> epsilon() const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        "epsilon",
+                        mesh_.time().timeName(),
+                        mesh_
+                    ),
+                    kt_*omega_,
+                    omega_.boundaryField().types()
+                )
+            );
+        }
+
+        //- Return the Reynolds stress tensor
+        virtual tmp<volSymmTensorField> R() const;
+
+        //- Return the effective stress tensor including the laminar stress
+        virtual tmp<volSymmTensorField> devReff() const;
+
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Solve the turbulence equations and correct the turbulence viscosity
+        virtual void correct();
+
+        //- Read RASProperties dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/Allrun b/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/Allrun
index 9210b03c39..caad132d7b 100755
--- a/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/Allrun
+++ b/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/Allrun
@@ -1,7 +1,6 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # run from this directory
 
-
 # Source tutorial run functions
 . $WM_PROJECT_DIR/bin/tools/RunFunctions
 
diff --git a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/changeDictionaryDict b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/changeDictionaryDict
index 971bafaa29..4ffe687151 100644
--- a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/changeDictionaryDict
+++ b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/changeDictionaryDict
@@ -27,7 +27,7 @@ dictionaryReplacement
                 D               700;
                 I               500;
                 length          1.05;
-                jump            uniform 0
+                jump            uniform 0;
                 value           uniform 0;
             }
             porous_half1
diff --git a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/mechanicalProperties b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/mechanicalProperties
index c92d02a0f7..377c6cee24 100644
--- a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/mechanicalProperties
+++ b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/mechanicalProperties
@@ -17,20 +17,20 @@ FoamFile
 
 rho
 {
-    rho         rhoInf;
-    rhoInf      rhoInf [ 1 -3 0 0 0 0 0 ] 7854;
+    type        uniform;
+    value       7854;
 }
 
 nu
 {
-    nu          nuInf;
-    nuInf       nuInf [ 0 0 0 0 0 0 0 ] 0.3;
+    type        uniform;
+    value       0.3;
 }
 
 E
 {
-    E           EInf;
-    EInf        EInf [ 1 -1 -2 0 0 0 0 ] 2e+11;
+    type        uniform;
+    value       2e+11;
 }
 
 planeStress     yes;
diff --git a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/thermalProperties b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/thermalProperties
index f68549dbba..77f422aca4 100644
--- a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/thermalProperties
+++ b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/thermalProperties
@@ -17,20 +17,20 @@ FoamFile
 
 C
 {
-    C           CInf;
-    CInf        CInf [ 0 2 -2 -1 0 0 0 ] 434;
+    type        uniform;
+    value       434;
 }
 
 k
 {
-    k           kInf;
-    kInf        kInf [ 1 1 -3 -1 0 0 0 ] 60.5;
+    type        uniform;
+    value       60.5;
 }
 
 alpha
 {
-    alpha       alphaInf;
-    alphaInf    alphaInf [ 0 0 0 -1 0 0 0 ] 1.1e-05;
+    type        uniform;
+    value       1.1e-05;
 }
 
 thermalStress   no;
diff --git a/tutorials/stressAnalysis/solidEquilibriumDisplacementFoam/beamEndLoad/constant/mechanicalProperties b/tutorials/stressAnalysis/solidEquilibriumDisplacementFoam/beamEndLoad/constant/mechanicalProperties
index f14ecc2556..8e28bddeb6 100644
--- a/tutorials/stressAnalysis/solidEquilibriumDisplacementFoam/beamEndLoad/constant/mechanicalProperties
+++ b/tutorials/stressAnalysis/solidEquilibriumDisplacementFoam/beamEndLoad/constant/mechanicalProperties
@@ -15,13 +15,24 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-rho             rho [ 1 -3 0 0 0 0 0 ] 7854;
+rho
+{
+    type        uniform;
+    value       7854;
+}
 
-E               E [ 1 -1 -2 0 0 0 0 ] 2e+11;
+nu
+{
+    type        uniform;
+    value       0.0;
+}
 
-nu              nu [ 0 0 0 0 0 0 0 ] 0;
+E
+{
+    type        uniform;
+    value       2e+11;
+}
 
 planeStress     yes;
 
-
 // ************************************************************************* //
diff --git a/tutorials/stressAnalysis/solidEquilibriumDisplacementFoam/beamEndLoad/constant/thermalProperties b/tutorials/stressAnalysis/solidEquilibriumDisplacementFoam/beamEndLoad/constant/thermalProperties
index 9d09fa7901..77f422aca4 100644
--- a/tutorials/stressAnalysis/solidEquilibriumDisplacementFoam/beamEndLoad/constant/thermalProperties
+++ b/tutorials/stressAnalysis/solidEquilibriumDisplacementFoam/beamEndLoad/constant/thermalProperties
@@ -15,11 +15,23 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-C               C [ 0 2 -2 -1 0 0 0 ] 434;
+C
+{
+    type        uniform;
+    value       434;
+}
 
-k               k [ 1 1 -3 -1 0 0 0 ] 60.5;
+k
+{
+    type        uniform;
+    value       60.5;
+}
 
-alpha           alpha [ 0 0 0 -1 0 0 0 ] 1.1e-05;
+alpha
+{
+    type        uniform;
+    value       1.1e-05;
+}
 
 thermalStress   no;
 
diff --git a/wmake/rules/linux64Gcc/c++ b/wmake/rules/linux64Gcc/c++
index 520e39c3b1..3dfb033360 100644
--- a/wmake/rules/linux64Gcc/c++
+++ b/wmake/rules/linux64Gcc/c++
@@ -17,5 +17,5 @@ cpptoo      = $(Ctoo)
 
 LINK_LIBS   = $(c++DBUG)
 
-LINKLIBSO   = $(CC) $(c++FLAGS) -shared
-LINKEXE     = $(CC) $(c++FLAGS) -Xlinker --add-needed
+LINKLIBSO   = $(CC) $(c++FLAGS) -shared -Xlinker --add-needed -Xlinker --no-as-needed
+LINKEXE     = $(CC) $(c++FLAGS) -Xlinker --add-needed -Xlinker --no-as-needed
diff --git a/wmake/rules/linux64Gcc45/c++ b/wmake/rules/linux64Gcc45/c++
index 3ca193e296..98b25ed1fe 100644
--- a/wmake/rules/linux64Gcc45/c++
+++ b/wmake/rules/linux64Gcc45/c++
@@ -17,5 +17,5 @@ cpptoo      = $(Ctoo)
 
 LINK_LIBS   = $(c++DBUG)
 
-LINKLIBSO   = $(CC) $(c++FLAGS) -shared
-LINKEXE     = $(CC) $(c++FLAGS) -Xlinker --add-needed
+LINKLIBSO   = $(CC) $(c++FLAGS) -shared -Xlinker --add-needed -Xlinker --no-as-needed
+LINKEXE     = $(CC) $(c++FLAGS) -Xlinker --add-needed -Xlinker --no-as-needed
diff --git a/wmake/rules/linux64Gcc46/c++ b/wmake/rules/linux64Gcc46/c++
index 3ca193e296..98b25ed1fe 100644
--- a/wmake/rules/linux64Gcc46/c++
+++ b/wmake/rules/linux64Gcc46/c++
@@ -17,5 +17,5 @@ cpptoo      = $(Ctoo)
 
 LINK_LIBS   = $(c++DBUG)
 
-LINKLIBSO   = $(CC) $(c++FLAGS) -shared
-LINKEXE     = $(CC) $(c++FLAGS) -Xlinker --add-needed
+LINKLIBSO   = $(CC) $(c++FLAGS) -shared -Xlinker --add-needed -Xlinker --no-as-needed
+LINKEXE     = $(CC) $(c++FLAGS) -Xlinker --add-needed -Xlinker --no-as-needed
diff --git a/wmake/rules/linuxGcc/c++ b/wmake/rules/linuxGcc/c++
index e862181fc5..357f4106e1 100644
--- a/wmake/rules/linuxGcc/c++
+++ b/wmake/rules/linuxGcc/c++
@@ -17,5 +17,5 @@ cpptoo      = $(Ctoo)
 
 LINK_LIBS   = $(c++DBUG)
 
-LINKLIBSO   = $(CC) $(c++FLAGS) -shared
-LINKEXE     = $(CC) $(c++FLAGS) -Xlinker --add-needed
+LINKLIBSO   = $(CC) $(c++FLAGS) -shared -Xlinker --add-needed -Xlinker --no-as-needed
+LINKEXE     = $(CC) $(c++FLAGS) -Xlinker --add-needed -Xlinker --no-as-needed
diff --git a/wmake/rules/linuxGcc45/c++ b/wmake/rules/linuxGcc45/c++
index e862181fc5..357f4106e1 100644
--- a/wmake/rules/linuxGcc45/c++
+++ b/wmake/rules/linuxGcc45/c++
@@ -17,5 +17,5 @@ cpptoo      = $(Ctoo)
 
 LINK_LIBS   = $(c++DBUG)
 
-LINKLIBSO   = $(CC) $(c++FLAGS) -shared
-LINKEXE     = $(CC) $(c++FLAGS) -Xlinker --add-needed
+LINKLIBSO   = $(CC) $(c++FLAGS) -shared -Xlinker --add-needed -Xlinker --no-as-needed
+LINKEXE     = $(CC) $(c++FLAGS) -Xlinker --add-needed -Xlinker --no-as-needed
diff --git a/wmake/rules/linuxGcc46/c++ b/wmake/rules/linuxGcc46/c++
index e862181fc5..357f4106e1 100644
--- a/wmake/rules/linuxGcc46/c++
+++ b/wmake/rules/linuxGcc46/c++
@@ -17,5 +17,5 @@ cpptoo      = $(Ctoo)
 
 LINK_LIBS   = $(c++DBUG)
 
-LINKLIBSO   = $(CC) $(c++FLAGS) -shared
-LINKEXE     = $(CC) $(c++FLAGS) -Xlinker --add-needed
+LINKLIBSO   = $(CC) $(c++FLAGS) -shared -Xlinker --add-needed -Xlinker --no-as-needed
+LINKEXE     = $(CC) $(c++FLAGS) -Xlinker --add-needed -Xlinker --no-as-needed