diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
index 52f23658a8..0adcd30927 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
@@ -23,16 +23,13 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
     scalar minDistSqr = magSqr(1e-6 * globalBb.span());
 
     // Non-empty directions
-    const Vector<label> validDirs = (mesh.directions() + Vector<label>::one)/2;
+    const Vector<label> validDirs = (mesh.geometricD() + Vector<label>::one)/2;
+    Info<< "    Mesh (non-empty, non-wedge) directions " << validDirs << endl;
 
-    Info<< "    Mesh (non-empty) directions " << validDirs << endl;
+    const Vector<label> solDirs = (mesh.solutionD() + Vector<label>::one)/2;
+    Info<< "    Mesh (non-empty) directions " << solDirs << endl;
 
-    scalar nGeomDims = mesh.nGeometricD();
-
-    Info<< "    Mesh (non-empty, non-wedge) dimensions "
-        << nGeomDims << endl;
-
-    if (nGeomDims < 3)
+    if (mesh.nGeometricD() < 3)
     {
         pointSet nonAlignedPoints(mesh, "nonAlignedEdges", mesh.nPoints()/100);
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
index 0428b252d7..e0b4c0536c 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
@@ -54,40 +54,79 @@ void Foam::polyMesh::calcDirections() const
 {
     for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
     {
-        directions_[cmpt] = 1;
+        solutionD_[cmpt] = 1;
     }
 
-    label nEmptyPatches = 0;
+    // Knock out empty and wedge directions. Note:they will be present on all
+    // domains.
 
-    vector dirVec = vector::zero;
+    label nEmptyPatches = 0;
+    label nWedgePatches = 0;
+
+    vector emptyDirVec = vector::zero;
+    vector wedgeDirVec = vector::zero;
 
     forAll(boundaryMesh(), patchi)
     {
-        if (isA<emptyPolyPatch>(boundaryMesh()[patchi]))
+        if (boundaryMesh()[patchi].size())
         {
-            if (boundaryMesh()[patchi].size())
+            if (isA<emptyPolyPatch>(boundaryMesh()[patchi]))
             {
                 nEmptyPatches++;
-                dirVec += sum(cmptMag(boundaryMesh()[patchi].faceAreas()));
+                emptyDirVec += sum(cmptMag(boundaryMesh()[patchi].faceAreas()));
+            }
+            else if (isA<wedgePolyPatch>(boundaryMesh()[patchi]))
+            {
+                const wedgePolyPatch& wpp = refCast<const wedgePolyPatch>
+                (
+                    boundaryMesh()[patchi]
+                );
+
+                nWedgePatches++;
+                wedgeDirVec += cmptMag(wpp.centreNormal());
             }
         }
     }
 
     if (nEmptyPatches)
     {
-        reduce(dirVec, sumOp<vector>());
+        reduce(emptyDirVec, sumOp<vector>());
 
-        dirVec /= mag(dirVec);
+        emptyDirVec /= mag(emptyDirVec);
 
         for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
         {
-            if (dirVec[cmpt] > 1e-6)
+            if (emptyDirVec[cmpt] > 1e-6)
             {
-                directions_[cmpt] = -1;
+                solutionD_[cmpt] = -1;
             }
             else
             {
-                directions_[cmpt] = 1;
+                solutionD_[cmpt] = 1;
+            }
+        }
+    }
+
+
+    // Knock out wedge directions
+
+    geometricD_ = solutionD_;
+
+    if (nWedgePatches)
+    {
+        reduce(wedgeDirVec, sumOp<vector>());
+
+        wedgeDirVec /= mag(wedgeDirVec);
+
+        for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+        {
+            if (wedgeDirVec[cmpt] > 1e-6)
+            {
+                geometricD_[cmpt] = -1;
+            }
+            else
+            {
+                geometricD_[cmpt] = 1;
             }
         }
     }
@@ -163,7 +202,8 @@ Foam::polyMesh::polyMesh(const IOobject& io)
         *this
     ),
     bounds_(points_),
-    directions_(Vector<label>::zero),
+    geometricD_(Vector<label>::zero),
+    solutionD_(Vector<label>::zero),
     pointZones_
     (
         IOobject
@@ -350,7 +390,8 @@ Foam::polyMesh::polyMesh
         0
     ),
     bounds_(points_, syncPar),
-    directions_(Vector<label>::zero),
+    geometricD_(Vector<label>::zero),
+    solutionD_(Vector<label>::zero),
     pointZones_
     (
         IOobject
@@ -505,7 +546,8 @@ Foam::polyMesh::polyMesh
         0
     ),
     bounds_(points_, syncPar),
-    directions_(Vector<label>::zero),
+    geometricD_(Vector<label>::zero),
+    solutionD_(Vector<label>::zero),
     pointZones_
     (
         IOobject
@@ -766,44 +808,37 @@ const Foam::fileName& Foam::polyMesh::facesInstance() const
 }
 
 
-const Foam::Vector<Foam::label>& Foam::polyMesh::directions() const
+const Foam::Vector<Foam::label>& Foam::polyMesh::geometricD() const
 {
-    if (directions_.x() == 0)
+    if (geometricD_.x() == 0)
     {
         calcDirections();
     }
 
-    return directions_;
+    return geometricD_;
 }
 
 
 Foam::label Foam::polyMesh::nGeometricD() const
 {
-    label nWedges = 0;
+    return cmptSum(geometricD() + Vector<label>::one)/2;
+}
 
-    forAll(boundary_, patchi)
+
+const Foam::Vector<Foam::label>& Foam::polyMesh::solutionD() const
+{
+    if (solutionD_.x() == 0)
     {
-        if (isA<wedgePolyPatch>(boundary_[patchi]))
-        {
-            nWedges++;
-        }
+        calcDirections();
     }
 
-    if (nWedges != 0 && nWedges != 2 && nWedges != 4)
-    {
-        FatalErrorIn("label polyMesh::nGeometricD() const")
-            << "Number of wedge patches " << nWedges << " is incorrect, "
-               "should be 0, 2 or 4"
-            << exit(FatalError);
-    }
-
-    return nSolutionD() - nWedges/2;
+    return solutionD_;
 }
 
 
 Foam::label Foam::polyMesh::nSolutionD() const
 {
-    return cmptSum(directions() + Vector<label>::one)/2;
+    return cmptSum(solutionD() + Vector<label>::one)/2;
 }
 
 
@@ -823,6 +858,10 @@ void Foam::polyMesh::addPatches
             << abort(FatalError);
     }
 
+    // Reset valid directions
+    geometricD_ = Vector<label>::zero;
+    solutionD_ = Vector<label>::zero;
+
     boundary_.setSize(p.size());
 
     // Copy the patch pointers
@@ -1037,6 +1076,10 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
     faceZones_.movePoints(points_);
     cellZones_.movePoints(points_);
 
+    // Reset valid directions (could change with rotation)
+    geometricD_ = Vector<label>::zero;
+    solutionD_ = Vector<label>::zero;
+
 
     // Hack until proper callbacks. Below are all the polyMeh MeshObjects with a
     // movePoints function.
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.H b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
index 588394d1c2..ecb805ad7f 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
@@ -120,9 +120,13 @@ private:
             //  Created from points on construction, updated when the mesh moves
             boundBox bounds_;
 
+            //- vector of non-constrained directions in mesh
+            //  defined according to the presence of empty and wedge patches
+            mutable Vector<label> geometricD_;
+
             //- vector of valid directions in mesh
             //  defined according to the presence of empty patches
-            mutable Vector<label> directions_;
+            mutable Vector<label> solutionD_;
 
 
         // Zoning information
@@ -309,17 +313,22 @@ public:
                 return bounds_;
             }
 
-            //- Return the vector of valid directions in mesh.
-            //  Defined according to the presence of empty patches.
-            //  1 indicates valid direction and -1 an invalid direction.
-            const Vector<label>& directions() const;
+            //- Return the vector of geometric directions in mesh.
+            //  Defined according to the presence of empty and wedge patches.
+            //  1 indicates unconstrained direction and -1 a constrained
+            //  direction.
+            const Vector<label>& geometricD() const;
 
             //- Return the number of valid geometric dimensions in the mesh
             label nGeometricD() const;
 
-            //- Return the number of valid solution dimensions in the mesh.
-            //  For wedge cases this includes the circumferential direction
-            //  in case of swirl.
+            //- Return the vector of solved-for directions in mesh.
+            //  Differs from geometricD in that it includes for wedge cases
+            //  the circumferential direction in case of swirl.
+            //  1 indicates valid direction and -1 an invalid direction.
+            const Vector<label>& solutionD() const;
+
+            //- Return the number of valid solved-for dimensions in the mesh
             label nSolutionD() const;
 
             //- Return point zone mesh
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C b/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
index 52f4446444..7183120984 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
@@ -68,6 +68,10 @@ void Foam::polyMesh::clearGeom()
         boundary_[patchI].clearGeom();
     }
 
+    // Reset valid directions (could change with rotation)
+    geometricD_ = Vector<label>::zero;
+    solutionD_ = Vector<label>::zero;
+
     pointMesh::Delete(*this);
 }
 
@@ -87,6 +91,10 @@ void Foam::polyMesh::clearAddressing()
     // recalculation
     deleteDemandDrivenData(globalMeshDataPtr_);
 
+    // Reset valid directions
+    geometricD_ = Vector<label>::zero;
+    solutionD_ = Vector<label>::zero;
+
     pointMesh::Delete(*this);
 }
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C
index d12dffff62..f1cc51cb9e 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C
@@ -217,7 +217,8 @@ Foam::polyMesh::polyMesh
         boundaryFaces.size() + 1    // add room for a default patch
     ),
     bounds_(points_, syncPar),
-    directions_(Vector<label>::zero),
+    geometricD_(Vector<label>::zero),
+    solutionD_(Vector<label>::zero),
     pointZones_
     (
         IOobject
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C b/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C
index 24d9ebbb13..7dcd3ba6b4 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C
@@ -68,6 +68,11 @@ void Foam::polyMesh::updateMesh(const mapPolyMesh& mpm)
         newMotionPoints.map(oldMotionPoints, mpm.pointMap());
     }
 
+    // Reset valid directions (could change by faces put into empty patches)
+    geometricD_ = Vector<label>::zero;
+    solutionD_ = Vector<label>::zero;
+
+
     // Hack until proper callbacks. Below are all the polyMesh-MeshObjects.
 
     // pointMesh
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C
index 133f1e7933..6bef8de896 100644
--- a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C
@@ -140,11 +140,11 @@ void Foam::quadraticFitSnGradData::findFaceDirs
     #ifndef SPHERICAL_GEOMETRY
         if (mesh.nGeometricD() <= 2) // find the normal direcion
         {
-            if (mesh.directions()[0] == -1)
+            if (mesh.geometricD()[0] == -1)
             {
                 kdir = vector(1, 0, 0);
             }
-            else if (mesh.directions()[1] == -1)
+            else if (mesh.geometricD()[1] == -1)
             {
                 kdir = vector(0, 1, 0);
             }
@@ -153,7 +153,7 @@ void Foam::quadraticFitSnGradData::findFaceDirs
                 kdir = vector(0, 0, 1);
             }
         }
-        else // 3D so find a direction in the place of the face
+        else // 3D so find a direction in the plane of the face
         {
             const face& f = mesh.faces()[faci];
             kdir = mesh.points()[f[0]] - mesh.points()[f[1]];
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index e1ef6135a6..2134420ead 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -713,7 +713,7 @@ Foam::fvMatrix<Type>::H() const
     (
         pow
         (
-            psi_.mesh().directions(),
+            psi_.mesh().solutionD(),
             pTraits<typename powProduct<Vector<label>, Type::rank>::type>::zero
         )
     );
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
index e8f4e6de49..21f3738a02 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
@@ -82,7 +82,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
     (
         pow
         (
-            psi_.mesh().directions(),
+            psi_.mesh().solutionD(),
             pTraits<typename powProduct<Vector<label>, Type::rank>::type>::zero
         )
     );
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C
index 0872cb81a9..28a3c41e64 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C
@@ -83,11 +83,11 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::findFaceDirs
 #   ifndef SPHERICAL_GEOMETRY
     if (mesh.nGeometricD() <= 2) // find the normal direction
     {
-        if (mesh.directions()[0] == -1)
+        if (mesh.geometricD()[0] == -1)
         {
             kdir = vector(1, 0, 0);
         }
-        else if (mesh.directions()[1] == -1)
+        else if (mesh.geometricD()[1] == -1)
         {
             kdir = vector(0, 1, 0);
         }
@@ -115,7 +115,7 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::findFaceDirs
 
         if (magk < SMALL)
         {
-            FatalErrorIn("findFaceDirs") << " calculated kdir = zero"
+            FatalErrorIn("findFaceDirs(..)") << " calculated kdir = zero"
                 << exit(FatalError);
         }
         else
diff --git a/src/meshTools/meshTools/meshTools.C b/src/meshTools/meshTools/meshTools.C
index a4f57a8de8..f826675b30 100644
--- a/src/meshTools/meshTools/meshTools.C
+++ b/src/meshTools/meshTools/meshTools.C
@@ -25,7 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "meshTools.H"
-#include "primitiveMesh.H"
+#include "polyMesh.H"
 #include "hexMatcher.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -635,6 +635,103 @@ Foam::label Foam::meshTools::walkFace
 }
 
 
+void Foam::meshTools::constrainToMeshCentre(const polyMesh& mesh, point& pt)
+{
+    const Vector<label>& dirs = mesh.geometricD();
+    const point& min = mesh.bounds().min();
+    const point& max = mesh.bounds().max();
+
+    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+    {
+        if (dirs[cmpt] == -1)
+        {
+            pt[cmpt] = 0.5*(min[cmpt]+max[cmpt]);
+        }
+    }
+}
+
+
+void Foam::meshTools::constrainToMeshCentre
+(
+    const polyMesh& mesh,
+    pointField& pts
+)
+{
+    const Vector<label>& dirs = mesh.geometricD();
+    const point& min = mesh.bounds().min();
+    const point& max = mesh.bounds().max();
+
+    bool isConstrained = false;
+    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+    {
+        if (dirs[cmpt] == -1)
+        {
+            isConstrained = true;
+            break;
+        }
+    }
+
+    if (isConstrained)
+    {
+        forAll(pts, i)
+        {
+            for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+            {
+                if (dirs[cmpt] == -1)
+                {
+                    pts[i][cmpt] = 0.5*(min[cmpt]+max[cmpt]);
+                }
+            }
+        }
+    }
+}
+
+
+//- Set the constrained components of directions/velocity to zero
+void Foam::meshTools::constrainDirection(const polyMesh& mesh, vector& d)
+{
+    const Vector<label>& dirs = mesh.geometricD();
+
+    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+    {
+        if (dirs[cmpt] == -1)
+        {
+            d[cmpt] = 0.0;
+        }
+    }
+}
+
+
+void Foam::meshTools::constrainDirection(const polyMesh& mesh, vectorField& d)
+{
+    const Vector<label>& dirs = mesh.geometricD();
+
+    bool isConstrained = false;
+    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+    {
+        if (dirs[cmpt] == -1)
+        {
+            isConstrained = true;
+            break;
+        }
+    }
+
+    if (isConstrained)
+    {
+        forAll(d, i)
+        {
+            for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+            {
+                if (dirs[cmpt] == -1)
+                {
+                    d[i][cmpt] = 0.0;
+                }
+            }
+        }
+    }
+}
+
+
 void Foam::meshTools::getParallelEdges
 (
     const primitiveMesh& mesh,
diff --git a/src/meshTools/meshTools/meshTools.H b/src/meshTools/meshTools/meshTools.H
index 304159e045..3710d5af06 100644
--- a/src/meshTools/meshTools/meshTools.H
+++ b/src/meshTools/meshTools/meshTools.H
@@ -50,7 +50,7 @@ namespace Foam
 {
 
 class primitiveMesh;
-
+class polyMesh;
 
 /*---------------------------------------------------------------------------*\
                         Namespace meshTools Declaration
@@ -81,196 +81,212 @@ namespace meshTools
     static const label pXpYpZMask = 1 << pXpYpZ;
 
 
-    //- Check if n is in same direction as normals of all faceLabels
-    bool visNormal
-    (
-        const vector& n,
-        const vectorField& faceNormals,
-        const labelList& faceLabels
-    );
+    // Normal handling
 
-    //- Calculate point normals on a 'box' mesh (all edges aligned with
-    //  coordinate axes)
-    vectorField calcBoxPointNormals(const primitivePatch& pp);
+        //- Check if n is in same direction as normals of all faceLabels
+        bool visNormal
+        (
+            const vector& n,
+            const vectorField& faceNormals,
+            const labelList& faceLabels
+        );
 
-    //- Normalized edge vector
-    vector normEdgeVec(const primitiveMesh&, const label edgeI);
+        //- Calculate point normals on a 'box' mesh (all edges aligned with
+        //  coordinate axes)
+        vectorField calcBoxPointNormals(const primitivePatch& pp);
 
-    //- Write obj representation of point
-    void writeOBJ
-    (
-        Ostream& os,
-        const point& pt
-    );
+        //- Normalized edge vector
+        vector normEdgeVec(const primitiveMesh&, const label edgeI);
 
 
-    //- Write obj representation of faces subset
-    void writeOBJ
-    (
-        Ostream& os,
-        const faceList&,
-        const pointField&,
-        const labelList& faceLabels
-    );
+    // OBJ writing
 
-    //- Write obj representation of faces
-    void writeOBJ
-    (
-        Ostream& os,
-        const faceList&,
-        const pointField&
-    );
+        //- Write obj representation of point
+        void writeOBJ
+        (
+            Ostream& os,
+            const point& pt
+        );
 
-    //- Write obj representation of cell subset
-    void writeOBJ
-    (
-        Ostream& os,
-        const cellList&,
-        const faceList&,
-        const pointField&,
-        const labelList& cellLabels
-    );
+        //- Write obj representation of faces subset
+        void writeOBJ
+        (
+            Ostream& os,
+            const faceList&,
+            const pointField&,
+            const labelList& faceLabels
+        );
 
-    //- Is edge used by cell
-    bool edgeOnCell
-    (
-        const primitiveMesh&,
-        const label cellI,
-        const label edgeI
-    );
+        //- Write obj representation of faces
+        void writeOBJ
+        (
+            Ostream& os,
+            const faceList&,
+            const pointField&
+        );
 
-    //- Is edge used by face
-    bool edgeOnFace
-    (
-        const primitiveMesh&,
-        const label faceI,
-        const label edgeI
-    );
-    
-    //- Is face used by cell
-    bool faceOnCell
-    (
-        const primitiveMesh&,
-        const label cellI,
-        const label faceI
-    );
+        //- Write obj representation of cell subset
+        void writeOBJ
+        (
+            Ostream& os,
+            const cellList&,
+            const faceList&,
+            const pointField&,
+            const labelList& cellLabels
+        );
 
-    //- Return edge among candidates that uses the two vertices.
-    label findEdge
-    (
-        const edgeList& edges,
-        const labelList& candidates,
-        const label v0,
-        const label v1
-    );
 
-    //- Return edge between two vertices. Returns -1 if no edge.
-    label findEdge
-    (
-        const primitiveMesh&,
-        const label v0,
-        const label v1
-    );
+    // Cell/face/edge walking
 
-    //- Return edge shared by two faces. Throws error if no edge found.
-    label getSharedEdge
-    (
-        const primitiveMesh&,
-        const label f0,
-        const label f1
-    );
-    
-    //- Return face shared by two cells. Throws error if none found.
-    label getSharedFace
-    (
-        const primitiveMesh&,
-        const label cell0,
-        const label cell1
-    );
+        //- Is edge used by cell
+        bool edgeOnCell
+        (
+            const primitiveMesh&,
+            const label cellI,
+            const label edgeI
+        );
 
-    //- Get faces on cell using edgeI. Throws error if no two found.
-    void getEdgeFaces
-    (
-        const primitiveMesh&,
-        const label cellI,
-        const label edgeI,
-        label& face0,
-        label& face1
-    );
+        //- Is edge used by face
+        bool edgeOnFace
+        (
+            const primitiveMesh&,
+            const label faceI,
+            const label edgeI
+        );
 
-    //- Return label of other edge (out of candidates edgeLabels)
-    //  connected to vertex but not edgeI. Throws error if none found.
-    label otherEdge
-    (
-        const primitiveMesh&,
-        const labelList& edgeLabels,
-        const label edgeI,
-        const label vertI
-    );
-    
-    //- Return face on cell using edgeI but not faceI. Throws error 
-    //  if none found.
-    label otherFace
-    (
-        const primitiveMesh&,
-        const label cellI,
-        const label faceI,
-        const label edgeI
-    );
+        //- Is face used by cell
+        bool faceOnCell
+        (
+            const primitiveMesh&,
+            const label cellI,
+            const label faceI
+        );
 
-    //- Return cell on other side of face. Throws error 
-    //  if face not internal.
-    label otherCell
-    (
-        const primitiveMesh&,
-        const label cellI,
-        const label faceI
-    );
-    
-    //- Returns label of edge nEdges away from startEdge (in the direction
-    // of startVertI)
-    label walkFace
-    (
-        const primitiveMesh&,
-        const label faceI,
-        const label startEdgeI,
-        const label startVertI,
-        const label nEdges
-    );
+        //- Return edge among candidates that uses the two vertices.
+        label findEdge
+        (
+            const edgeList& edges,
+            const labelList& candidates,
+            const label v0,
+            const label v1
+        );
+
+        //- Return edge between two vertices. Returns -1 if no edge.
+        label findEdge
+        (
+            const primitiveMesh&,
+            const label v0,
+            const label v1
+        );
+
+        //- Return edge shared by two faces. Throws error if no edge found.
+        label getSharedEdge
+        (
+            const primitiveMesh&,
+            const label f0,
+            const label f1
+        );
+
+        //- Return face shared by two cells. Throws error if none found.
+        label getSharedFace
+        (
+            const primitiveMesh&,
+            const label cell0,
+            const label cell1
+        );
+
+        //- Get faces on cell using edgeI. Throws error if no two found.
+        void getEdgeFaces
+        (
+            const primitiveMesh&,
+            const label cellI,
+            const label edgeI,
+            label& face0,
+            label& face1
+        );
+
+        //- Return label of other edge (out of candidates edgeLabels)
+        //  connected to vertex but not edgeI. Throws error if none found.
+        label otherEdge
+        (
+            const primitiveMesh&,
+            const labelList& edgeLabels,
+            const label edgeI,
+            const label vertI
+        );
+
+        //- Return face on cell using edgeI but not faceI. Throws error 
+        //  if none found.
+        label otherFace
+        (
+            const primitiveMesh&,
+            const label cellI,
+            const label faceI,
+            const label edgeI
+        );
+
+        //- Return cell on other side of face. Throws error 
+        //  if face not internal.
+        label otherCell
+        (
+            const primitiveMesh&,
+            const label cellI,
+            const label faceI
+        );
+
+        //- Returns label of edge nEdges away from startEdge (in the direction
+        // of startVertI)
+        label walkFace
+        (
+            const primitiveMesh&,
+            const label faceI,
+            const label startEdgeI,
+            const label startVertI,
+            const label nEdges
+        );
+
+
+    // Constraints on position
+
+        //- Set the constrained components of position to mesh centre
+        void constrainToMeshCentre(const polyMesh&, point&);
+        void constrainToMeshCentre(const polyMesh&, pointField&);
+
+        //- Set the constrained components of directions/velocity to zero
+        void constrainDirection(const polyMesh&, vector&);
+        void constrainDirection(const polyMesh&, vectorField&);
 
     
-    //
     // Hex only functionality.
-    //
 
-    //- Given edge on hex find other 'parallel', non-connected edges.
-    void getParallelEdges
-    (
-        const primitiveMesh&,
-        const label cellI,
-        const label e0,
-        label&,
-        label&,
-        label&
-    );
+        //- Given edge on hex find other 'parallel', non-connected edges.
+        void getParallelEdges
+        (
+            const primitiveMesh&,
+            const label cellI,
+            const label e0,
+            label&,
+            label&,
+            label&
+        );
 
-    //- Given edge on hex find all 'parallel' (i.e. non-connected)
-    //  edges and average direction of them
-    vector edgeToCutDir
-    (
-        const primitiveMesh&,
-        const label cellI,
-        const label edgeI
-    );
+        //- Given edge on hex find all 'parallel' (i.e. non-connected)
+        //  edges and average direction of them
+        vector edgeToCutDir
+        (
+            const primitiveMesh&,
+            const label cellI,
+            const label edgeI
+        );
 
-    //- Reverse of edgeToCutDir: given direction find edge bundle and
-    //  return one of them.
-    label cutDirToEdge
-    (
-        const primitiveMesh&,
-        const label cellI,
-        const vector& cutDir
-    );
+        //- Reverse of edgeToCutDir: given direction find edge bundle and
+        //  return one of them.
+        label cutDirToEdge
+        (
+            const primitiveMesh&,
+            const label cellI,
+            const vector& cutDir
+        );
     
 } // End namespace meshTools
 
diff --git a/src/turbulenceModels/LES/LESdeltas/cubeRootVolDelta/cubeRootVolDelta.C b/src/turbulenceModels/LES/LESdeltas/cubeRootVolDelta/cubeRootVolDelta.C
index a3ce3ab3a9..2002deaa21 100644
--- a/src/turbulenceModels/LES/LESdeltas/cubeRootVolDelta/cubeRootVolDelta.C
+++ b/src/turbulenceModels/LES/LESdeltas/cubeRootVolDelta/cubeRootVolDelta.C
@@ -42,8 +42,7 @@ addToRunTimeSelectionTable(LESdelta, cubeRootVolDelta, dictionary);
 
 void cubeRootVolDelta::calcDelta()
 {
-    const Vector<label>& directions = mesh().directions();
-    label nD = (directions.nComponents + cmptSum(directions))/2;
+    label nD = mesh().nGeometricD();
 
     if (nD == 3)
     {
@@ -55,14 +54,15 @@ void cubeRootVolDelta::calcDelta()
             << "Case is 2D, LES is not strictly applicable\n"
             << endl;
 
+        const Vector<label>& directions = mesh().geometricD();
+
         scalar thickness = 0.0;
         for (direction dir=0; dir<directions.nComponents; dir++)
         {
             if (directions[dir] == -1)
             {
-                boundBox bb(mesh().points(), false);
-
-                thickness = bb.span()[dir];
+                thickness = mesh().bounds().span()[dir];
+                break;
             }
         }
 
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.C b/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.C
index 3f3aa7699b..005ee2eb16 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.C
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.C
@@ -41,8 +41,7 @@ namespace Foam
 
 void Foam::IDDESDelta::calcDelta()
 {
-    const Vector<label>& directions = mesh().directions();
-    label nD = (directions.nComponents + cmptSum(directions))/2;
+    label nD = mesh().nGeometricD();
 
     // initialise hwn as wall distance
     volScalarField hwn = wallDist(mesh()).y();