From ceecb84573fd4c53f1b48cd9b6600a24089e9a35 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 17 Dec 2012 13:45:19 +0000
Subject: [PATCH] ENH: checkMesh: moved some checking functionality to polyMesh

---
 .../manipulation/checkMesh/checkGeometry.C    |    2 +-
 src/OpenFOAM/Make/files                       |    6 +-
 src/OpenFOAM/meshes/polyMesh/polyMesh.C       |    2 +-
 src/OpenFOAM/meshes/polyMesh/polyMesh.H       |   81 ++
 .../polyMesh/polyMeshCheck/polyMeshCheck.C    |  679 ++++++++++
 .../polyMesh/polyMeshCheck/polyMeshTools.C    |  190 +++
 .../polyMesh/polyMeshCheck/polyMeshTools.H    |  109 ++
 .../meshes/primitiveMesh/primitiveMesh.H      |  223 ++--
 .../primitiveMeshCheck/primitiveMeshCheck.C   | 1178 +++++++----------
 .../primitiveMeshCheckMotion.C                |  266 ----
 .../primitiveMeshCheck/primitiveMeshTools.C   |  638 +++++++++
 .../primitiveMeshCheck/primitiveMeshTools.H   |  168 +++
 12 files changed, 2467 insertions(+), 1075 deletions(-)
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H
 delete mode 100644 src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheckMotion.C
 create mode 100644 src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C
 create mode 100644 src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H

diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
index d964ac12cc..ef521a257d 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
@@ -803,7 +803,7 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
     if (allGeometry)
     {
         cellSet cells(mesh, "underdeterminedCells", mesh.nCells()/100);
-        if (mesh.checkCellDeterminant(true, &cells, mesh.geometricD()))
+        if (mesh.checkCellDeterminant(true, &cells))
         {
             noFailedChecks++;
 
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index 1cd7f37fe8..55208034a1 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -425,6 +425,10 @@ $(polyMesh)/polyMeshInitMesh.C
 $(polyMesh)/polyMeshClear.C
 $(polyMesh)/polyMeshUpdate.C
 
+polyMeshCheck = $(polyMesh)/polyMeshCheck
+$(polyMeshCheck)/polyMeshCheck.C
+$(polyMeshCheck)/polyMeshTools.C
+
 primitiveMesh = meshes/primitiveMesh
 $(primitiveMesh)/primitiveMesh.C
 $(primitiveMesh)/primitiveMeshCellCells.C
@@ -445,9 +449,9 @@ $(primitiveMesh)/primitiveMeshCalcCellShapes.C
 
 primitiveMeshCheck = $(primitiveMesh)/primitiveMeshCheck
 $(primitiveMeshCheck)/primitiveMeshCheck.C
-$(primitiveMeshCheck)/primitiveMeshCheckMotion.C
 $(primitiveMeshCheck)/primitiveMeshCheckPointNearness.C
 $(primitiveMeshCheck)/primitiveMeshCheckEdgeLength.C
+$(primitiveMeshCheck)/primitiveMeshTools.C
 
 primitivePatch = $(primitiveMesh)/primitivePatch
 $(primitivePatch)/patchZones.C
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
index 6388190bf6..abe537320f 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
@@ -1129,7 +1129,7 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
     if (debug)
     {
         // Check mesh motion
-        if (primitiveMesh::checkMeshMotion(points_, true))
+        if (checkMeshMotion(points_, true))
         {
             Info<< "tmp<scalarField> polyMesh::movePoints"
                 << "(const pointField&) : "
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.H b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
index 09c046feb8..e4c1897671 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
@@ -34,6 +34,7 @@ SourceFiles
     polyMeshFromShapeMesh.C
     polyMeshIO.C
     polyMeshUpdate.C
+    polyMeshCheck.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -224,7 +225,45 @@ private:
                 cellList& cells
             );
 
+        // Geometry checks
 
+            //- Check non-orthogonality
+            bool checkFaceOrthogonality
+            (
+                const vectorField& fAreas,
+                const vectorField& cellCtrs,
+                const bool report,
+                const bool detailedReport,
+                labelHashSet* setPtr
+            ) const;
+
+            //- Check face skewness
+            bool checkFaceSkewness
+            (
+                const pointField& points,
+                const vectorField& fCtrs,
+                const vectorField& fAreas,
+                const vectorField& cellCtrs,
+                const bool report,
+                const bool detailedReport,
+                labelHashSet* setPtr
+            ) const;
+
+            bool checkEdgeAlignment
+            (
+                const pointField& p,
+                const bool report,
+                const Vector<label>& directions,
+                labelHashSet* setPtr
+            ) const;
+
+            bool checkCellDeterminant
+            (
+                const vectorField& faceAreas,
+                const bool report,
+                labelHashSet* setPtr,
+                const Vector<label>& meshD
+            ) const;
 
 public:
 
@@ -535,6 +574,48 @@ public:
             void removeFiles() const;
 
 
+        // Geometric checks. Selectively override primitiveMesh functionality.
+
+            //- Check boundary for closedness
+            virtual bool checkClosedBoundary(const bool report = false) const;
+
+            //- Check non-orthogonality
+            virtual bool checkFaceOrthogonality
+            (
+                const bool report = false,
+                labelHashSet* setPtr = NULL
+            ) const;
+
+            //- Check face skewness
+            virtual bool checkFaceSkewness
+            (
+                const bool report = false,
+                labelHashSet* setPtr = NULL
+            ) const;
+
+            //- Check edge alignment for 1D/2D cases
+            virtual bool checkEdgeAlignment
+            (
+                const bool report,
+                const Vector<label>& directions,
+                labelHashSet* setPtr
+            ) const;
+
+            virtual bool checkCellDeterminant
+            (
+                const bool report,
+                labelHashSet* setPtr
+            ) const;
+
+            //- Check mesh motion for correctness given motion points
+            virtual bool checkMeshMotion
+            (
+                const pointField& newPoints,
+                const bool report = false,
+                const bool detailedReport = false
+            ) const;
+
+
         // Helper functions
 
             //- Find the cell, tetFaceI and tetPtI for the given position
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C
new file mode 100644
index 0000000000..748c7da21a
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C
@@ -0,0 +1,679 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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 "polyMesh.H"
+#include "polyMeshTools.H"
+#include "unitConversion.H"
+#include "syncTools.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::polyMesh::checkFaceOrthogonality
+(
+    const vectorField& fAreas,
+    const vectorField& cellCtrs,
+    const bool report,
+    const bool detailedReport,
+    labelHashSet* setPtr
+) const
+{
+    if (debug)
+    {
+        Info<< "bool polyMesh::checkFaceOrthogonality("
+            << "const bool, labelHashSet*) const: "
+            << "checking mesh non-orthogonality" << endl;
+    }
+
+    const labelList& own = faceOwner();
+    const labelList& nei = faceNeighbour();
+
+    // Calculate orthogonality for all internal and coupled boundary faces
+    // (1 for uncoupled boundary faces)
+    tmp<scalarField> tortho = polyMeshTools::faceOrthogonality
+    (
+        *this,
+        fAreas,
+        cellCtrs
+    );
+    const scalarField& ortho = tortho();
+
+    // Severe nonorthogonality threshold
+    const scalar severeNonorthogonalityThreshold =
+        ::cos(degToRad(primitiveMesh::nonOrthThreshold_));
+
+
+    scalar minDDotS = GREAT;
+    scalar sumDDotS = 0.0;
+    label nSummed = 0;
+    label severeNonOrth = 0;
+    label errorNonOrth = 0;
+
+
+    // Statistics only for internal and masters of coupled faces
+    PackedBoolList isMasterFace(syncTools::getInternalOrMasterFaces(*this));
+
+    forAll(ortho, faceI)
+    {
+        if (ortho[faceI] < severeNonorthogonalityThreshold)
+        {
+            if (ortho[faceI] > SMALL)
+            {
+                if (setPtr)
+                {
+                    setPtr->insert(faceI);
+                }
+
+                severeNonOrth++;
+            }
+            else
+            {
+                // Error : non-ortho too large
+                if (setPtr)
+                {
+                    setPtr->insert(faceI);
+                }
+                if (detailedReport && errorNonOrth == 0)
+                {
+                    // Non-orthogonality greater than 90 deg
+                    WarningIn
+                    (
+                        "polyMesh::checkFaceOrthogonality"
+                        "(const pointField&, const bool) const"
+                    )   << "Severe non-orthogonality for face "
+                        << faceI
+                        << " between cells " << own[faceI]
+                        << " and " << nei[faceI]
+                        << ": Angle = " << radToDeg(::acos(ortho[faceI]))
+                        << " deg." << endl;
+                }
+
+                errorNonOrth++;
+            }
+        }
+
+        if (isMasterFace[faceI])
+        {
+            minDDotS = min(minDDotS, ortho[faceI]);
+            sumDDotS += ortho[faceI];
+            nSummed++;
+        }
+    }
+
+    reduce(minDDotS, minOp<scalar>());
+    reduce(sumDDotS, sumOp<scalar>());
+    reduce(nSummed, sumOp<label>());
+    reduce(severeNonOrth, sumOp<label>());
+    reduce(errorNonOrth, sumOp<label>());
+
+    if (debug || report)
+    {
+        if (nSummed > 0)
+        {
+            if (debug || report)
+            {
+                Info<< "    Mesh non-orthogonality Max: "
+                    << radToDeg(::acos(minDDotS))
+                    << " average: " << radToDeg(::acos(sumDDotS/nSummed))
+                    << endl;
+            }
+        }
+
+        if (severeNonOrth > 0)
+        {
+            Info<< "   *Number of severely non-orthogonal faces: "
+                << severeNonOrth << "." << endl;
+        }
+    }
+
+    if (errorNonOrth > 0)
+    {
+        if (debug || report)
+        {
+            Info<< " ***Number of non-orthogonality errors: "
+                << errorNonOrth << "." << endl;
+        }
+
+        return true;
+    }
+    else
+    {
+        if (debug || report)
+        {
+            Info<< "    Non-orthogonality check OK." << endl;
+        }
+
+        return false;
+    }
+}
+
+
+bool Foam::polyMesh::checkFaceSkewness
+(
+    const pointField& points,
+    const vectorField& fCtrs,
+    const vectorField& fAreas,
+    const vectorField& cellCtrs,
+    const bool report,
+    const bool detailedReport,
+    labelHashSet* setPtr
+) const
+{
+    if (debug)
+    {
+        Info<< "bool polyMesh::checkFaceSkewnesss("
+            << "const bool, labelHashSet*) const: "
+            << "checking face skewness" << endl;
+    }
+
+    const labelList& own = faceOwner();
+    const labelList& nei = faceNeighbour();
+
+    // Warn if the skew correction vector is more than skewWarning times
+    // larger than the face area vector
+
+    tmp<scalarField> tskew = polyMeshTools::faceSkewness
+    (
+        *this,
+        points,
+        fCtrs,
+        fAreas,
+        cellCtrs
+    );
+    const scalarField& skew = tskew();
+
+    scalar maxSkew = max(skew);
+    label nWarnSkew = 0;
+
+    // Statistics only for all faces except slave coupled faces
+    PackedBoolList isMasterFace(syncTools::getMasterFaces(*this));
+
+    forAll(skew, faceI)
+    {
+        // Check if the skewness vector is greater than the PN vector.
+        // This does not cause trouble but is a good indication of a poor mesh.
+        if (skew[faceI] > skewThreshold_)
+        {
+            if (setPtr)
+            {
+                setPtr->insert(faceI);
+            }
+            if (detailedReport && nWarnSkew == 0)
+            {
+                // Non-orthogonality greater than 90 deg
+                if (isInternalFace(faceI))
+                {
+                    WarningIn
+                    (
+                        "polyMesh::checkFaceSkewnesss"
+                        "(const pointField&, const bool) const"
+                    )   << "Severe skewness " << skew[faceI]
+                        << " for face " << faceI
+                        << " between cells " << own[faceI]
+                        << " and " << nei[faceI];
+                }
+                else
+                {
+                    WarningIn
+                    (
+                        "polyMesh::checkFaceSkewnesss"
+                        "(const pointField&, const bool) const"
+                    )   << "Severe skewness " << skew[faceI]
+                        << " for boundary face " << faceI
+                        << " on cell " << own[faceI];
+                }
+            }
+
+            if (isMasterFace[faceI])
+            {
+                nWarnSkew++;
+            }
+        }
+    }
+
+    reduce(maxSkew, maxOp<scalar>());
+    reduce(nWarnSkew, sumOp<label>());
+
+    if (nWarnSkew > 0)
+    {
+        if (debug || report)
+        {
+            Info<< " ***Max skewness = " << maxSkew
+                << ", " << nWarnSkew << " highly skew faces detected"
+                   " which may impair the quality of the results"
+                << endl;
+        }
+
+        return true;
+    }
+    else
+    {
+        if (debug || report)
+        {
+            Info<< "    Max skewness = " << maxSkew << " OK." << endl;
+        }
+
+        return false;
+    }
+}
+
+
+// Check 1D/2Dness of edges. Gets passed the non-empty directions and
+// checks all edges in the mesh whether they:
+// - have no component in a non-empty direction or
+// - are only in a singe non-empty direction.
+// Empty direction info is passed in as a vector of labels (synchronised)
+// which are 1 if the direction is non-empty, 0 if it is.
+bool Foam::polyMesh::checkEdgeAlignment
+(
+    const pointField& p,
+    const bool report,
+    const Vector<label>& directions,
+    labelHashSet* setPtr
+) const
+{
+    if (debug)
+    {
+        Info<< "bool polyMesh::checkEdgeAlignment("
+            << "const bool, const Vector<label>&, labelHashSet*) const: "
+            << "checking edge alignment" << endl;
+    }
+
+    label nDirs = 0;
+    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+    {
+        if (directions[cmpt] == 1)
+        {
+            nDirs++;
+        }
+        else if (directions[cmpt] != 0)
+        {
+            FatalErrorIn
+            (
+                "polyMesh::checkEdgeAlignment"
+                "(const bool, const Vector<label>&, labelHashSet*)"
+            )   << "directions should contain 0 or 1 but is now " << directions
+                << exit(FatalError);
+        }
+    }
+
+    if (nDirs == vector::nComponents)
+    {
+        return false;
+    }
+
+
+    const faceList& fcs = faces();
+
+    EdgeMap<label> edgesInError;
+
+    forAll(fcs, faceI)
+    {
+        const face& f = fcs[faceI];
+
+        forAll(f, fp)
+        {
+            label p0 = f[fp];
+            label p1 = f.nextLabel(fp);
+            if (p0 < p1)
+            {
+                vector d(p[p1]-p[p0]);
+                scalar magD = mag(d);
+
+                if (magD > ROOTVSMALL)
+                {
+                    d /= magD;
+
+                    // Check how many empty directions are used by the edge.
+                    label nEmptyDirs = 0;
+                    label nNonEmptyDirs = 0;
+                    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+                    {
+                        if (mag(d[cmpt]) > 1e-6)
+                        {
+                            if (directions[cmpt] == 0)
+                            {
+                                nEmptyDirs++;
+                            }
+                            else
+                            {
+                                nNonEmptyDirs++;
+                            }
+                        }
+                    }
+
+                    if (nEmptyDirs == 0)
+                    {
+                        // Purely in ok directions.
+                    }
+                    else if (nEmptyDirs == 1)
+                    {
+                        // Ok if purely in empty directions.
+                        if (nNonEmptyDirs > 0)
+                        {
+                            edgesInError.insert(edge(p0, p1), faceI);
+                        }
+                    }
+                    else if (nEmptyDirs > 1)
+                    {
+                        // Always an error
+                        edgesInError.insert(edge(p0, p1), faceI);
+                    }
+                }
+            }
+        }
+    }
+
+    label nErrorEdges = returnReduce(edgesInError.size(), sumOp<label>());
+
+    if (nErrorEdges > 0)
+    {
+        if (debug || report)
+        {
+            Info<< " ***Number of edges not aligned with or perpendicular to "
+                << "non-empty directions: " << nErrorEdges << endl;
+        }
+
+        if (setPtr)
+        {
+            setPtr->resize(2*edgesInError.size());
+            forAllConstIter(EdgeMap<label>, edgesInError, iter)
+            {
+                setPtr->insert(iter.key()[0]);
+                setPtr->insert(iter.key()[1]);
+            }
+        }
+
+        return true;
+    }
+    else
+    {
+        if (debug || report)
+        {
+            Info<< "    All edges aligned with or perpendicular to "
+                << "non-empty directions." << endl;
+        }
+        return false;
+    }
+}
+
+
+bool Foam::polyMesh::checkCellDeterminant
+(
+    const vectorField& faceAreas,
+    const bool report,
+    labelHashSet* setPtr,
+    const Vector<label>& meshD
+) const
+{
+    if (debug)
+    {
+        Info<< "bool polyMesh::checkCellDeterminant(const bool"
+            << ", labelHashSet*) const: "
+            << "checking for under-determined cells" << endl;
+    }
+
+    // Determine number of dimensions and (for 2D) missing dimension
+    label nDims = 0;
+    label twoD = -1;
+    for (direction dir = 0; dir < vector::nComponents; dir++)
+    {
+        if (meshD[dir] == 1)
+        {
+            nDims++;
+        }
+        else
+        {
+            twoD = dir;
+        }
+    }
+
+
+    tmp<scalarField> tcellDeterminant = primitiveMeshTools::cellDeterminant
+    (
+        *this,
+        meshD,
+        faceAreas,
+        syncTools::getInternalOrCoupledFaces(*this)
+    );
+    scalarField& cellDeterminant = tcellDeterminant();
+
+
+    label nErrorCells = 0;
+    scalar minDet = min(cellDeterminant);
+    scalar sumDet = sum(cellDeterminant);
+
+    forAll (cellDeterminant, cellI)
+    {
+        if (cellDeterminant[cellI] < 1e-3)
+        {
+            if (setPtr)
+            {
+                setPtr->insert(cellI);
+            }
+
+            nErrorCells++;
+        }
+    }
+
+    reduce(nErrorCells, sumOp<label>());
+    reduce(minDet, minOp<scalar>());
+    reduce(sumDet, sumOp<scalar>());
+    label nSummed = returnReduce(cellDeterminant.size(), sumOp<label>());
+
+    if (debug || report)
+    {
+        if (nSummed > 0)
+        {
+            Info<< "    Cell determinant (wellposedness) : minimum: " << minDet
+                << " average: " << sumDet/nSummed
+                << endl;
+        }
+    }
+
+    if (nErrorCells > 0)
+    {
+        if (debug || report)
+        {
+            Info<< " ***Cells with small determinant found, number of cells: "
+                << nErrorCells << endl;
+        }
+
+        return true;
+    }
+    else
+    {
+        if (debug || report)
+        {
+            Info<< "    Cell determinant check OK." << endl;
+        }
+
+        return false;
+    }
+
+    return false;
+}
+
+
+bool Foam::polyMesh::checkClosedBoundary(const bool report) const
+{
+    return primitiveMesh::checkClosedBoundary
+    (
+        faceAreas(),
+        report,
+        syncTools::getInternalOrCoupledFaces(*this)
+    );
+}
+
+
+bool Foam::polyMesh::checkFaceOrthogonality
+(
+    const bool report,
+    labelHashSet* setPtr
+) const
+{
+    return checkFaceOrthogonality
+    (
+        faceAreas(),
+        cellCentres(),
+        report,
+        false,  // detailedReport
+        setPtr
+    );
+}
+
+
+bool Foam::polyMesh::checkFaceSkewness
+(
+    const bool report,
+    labelHashSet* setPtr
+) const
+{
+    return checkFaceSkewness
+    (
+        points(),
+        faceCentres(),
+        faceAreas(),
+        cellCentres(),
+        report,
+        false,  // detailedReport
+        setPtr
+    );
+}
+
+
+bool Foam::polyMesh::checkEdgeAlignment
+(
+    const bool report,
+    const Vector<label>& directions,
+    labelHashSet* setPtr
+) const
+{
+    return checkEdgeAlignment
+    (
+        points(),
+        report,
+        directions,
+        setPtr
+    );
+}
+
+
+bool Foam::polyMesh::checkCellDeterminant
+(
+    const bool report,
+    labelHashSet* setPtr
+) const
+{
+    return checkCellDeterminant
+    (
+        faceAreas(),
+        report,
+        setPtr,
+        geometricD()
+    );
+}
+
+
+bool Foam::polyMesh::checkMeshMotion
+(
+    const pointField& newPoints,
+    const bool report,
+    const bool detailedReport
+) const
+{
+    if (debug || report)
+    {
+        Pout<< "bool polyMesh::checkMeshMotion("
+            << "const pointField&, const bool, const bool) const: "
+            << "checking mesh motion" << endl;
+    }
+
+    vectorField fCtrs(nFaces());
+    vectorField fAreas(nFaces());
+
+    makeFaceCentresAndAreas(newPoints, fCtrs, fAreas);
+
+    // Check cell volumes and calculate new cell centres
+    vectorField cellCtrs(nCells());
+    scalarField cellVols(nCells());
+
+    makeCellCentresAndVols(fCtrs, fAreas, cellCtrs, cellVols);
+
+    // Check cell volumes
+    bool error = checkCellVolumes
+    (
+        cellVols,       // vols
+        report,         // report
+        detailedReport, // detailedReport
+        NULL            // setPtr
+    );
+
+
+    // Check face areas
+    bool areaError = checkFaceAreas
+    (
+        faceAreas(),
+        report,         // report
+        detailedReport, // detailedReport,
+        NULL            // setPtr
+    );
+    error = error || areaError;
+
+
+    // Check pyramid volumes
+    bool pyrVolError = checkFacePyramids
+    (
+        newPoints,
+        cellCtrs,
+        report,         // report,
+        detailedReport, // detailedReport,
+        -SMALL,         // minPyrVol
+        NULL            // setPtr
+    );
+    error = error || pyrVolError;
+
+
+    // Check face non-orthogonality
+    bool nonOrthoError = checkFaceOrthogonality
+    (
+        fAreas,
+        cellCtrs,
+        report,         // report
+        detailedReport, // detailedReport
+        NULL            // setPtr
+    );
+    error = error || nonOrthoError;
+
+
+    if (!error && (debug || report))
+    {
+        Pout<< "Mesh motion check OK." << endl;
+    }
+
+    return error;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C
new file mode 100644
index 0000000000..fc0739e79f
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C
@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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 "polyMeshTools.H"
+#include "syncTools.H"
+#include "pyramidPointFaceRef.H"
+#include "primitiveMeshTools.H"
+#include "polyMeshTools.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceOrthogonality
+(
+    const polyMesh& mesh,
+    const vectorField& areas,
+    const vectorField& cc
+)
+{
+    const labelList& own = mesh.faceOwner();
+    const labelList& nei = mesh.faceNeighbour();
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+    tmp<scalarField> tortho(new scalarField(mesh.nFaces(), 1.0));
+    scalarField& ortho = tortho();
+
+    // Internal faces
+    forAll(nei, faceI)
+    {
+        ortho[faceI] = primitiveMeshTools::faceOrthogonality
+        (
+            cc[own[faceI]],
+            cc[nei[faceI]],
+            areas[faceI]
+        );
+    }
+
+
+    // Coupled faces
+
+    pointField neighbourCc;
+    syncTools::swapBoundaryCellList(mesh, cc, neighbourCc);
+
+    forAll(pbm, patchI)
+    {
+        const polyPatch& pp = pbm[patchI];
+        if (pp.coupled())
+        {
+            forAll(pp, i)
+            {
+                label faceI = pp.start() + i;
+                label bFaceI = faceI - mesh.nInternalFaces();
+
+                ortho[faceI] = primitiveMeshTools::faceOrthogonality
+                (
+                    cc[own[faceI]],
+                    neighbourCc[bFaceI],
+                    areas[faceI]
+                );
+            }
+        }
+    }
+
+    return tortho;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceSkewness
+(
+    const polyMesh& mesh,
+    const pointField& p,
+    const vectorField& fCtrs,
+    const vectorField& fAreas,
+    const vectorField& cellCtrs
+)
+{
+    const labelList& own = mesh.faceOwner();
+    const labelList& nei = mesh.faceNeighbour();
+    const faceList& fcs = mesh.faces();
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+    tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
+    scalarField& skew = tskew();
+
+    forAll(nei, faceI)
+    {
+        skew[faceI] = primitiveMeshTools::faceSkewness
+        (
+            mesh,
+            p,
+            fCtrs,
+            fAreas,
+
+            faceI,
+            cellCtrs[own[faceI]],
+            cellCtrs[nei[faceI]]
+        );
+    }
+
+
+    // Boundary faces: consider them to have only skewness error.
+    // (i.e. treat as if mirror cell on other side)
+
+    pointField neighbourCc;
+    syncTools::swapBoundaryCellList(mesh, cellCtrs, neighbourCc);
+
+    forAll(pbm, patchI)
+    {
+        const polyPatch& pp = pbm[patchI];
+        if (pp.coupled())
+        {
+            forAll(pp, i)
+            {
+                label faceI = pp.start() + i;
+                label bFaceI = faceI - mesh.nInternalFaces();
+
+                skew[faceI] = primitiveMeshTools::faceSkewness
+                (
+                    mesh,
+                    p,
+                    fCtrs,
+                    fAreas,
+
+                    faceI,
+                    cellCtrs[own[faceI]],
+                    neighbourCc[bFaceI]
+                );
+            }
+        }
+        else
+        {
+            forAll(pp, i)
+            {
+                label faceI = pp.start() + i;
+
+                vector Cpf = fCtrs[faceI] - cellCtrs[own[faceI]];
+
+                vector normal = fAreas[faceI];
+                normal /= mag(normal) + VSMALL;
+                vector d = normal*(normal & Cpf);
+
+
+                // Skewness vector
+                vector sv =
+                    Cpf
+                  - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d;
+                vector svHat = sv/(mag(sv) + VSMALL);
+
+                // Normalisation distance calculated as the approximate distance
+                // from the face centre to the edge of the face in the direction
+                // of the skewness
+                scalar fd = 0.4*mag(d) + VSMALL;
+                const face& f = fcs[faceI];
+                forAll(f, pi)
+                {
+                    fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI])));
+                }
+
+                // Normalised skewness
+                 skew[faceI] = mag(sv)/fd;
+            }
+        }
+    }
+
+    return tskew;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H
new file mode 100644
index 0000000000..53c758c3be
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H
@@ -0,0 +1,109 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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/>.
+
+Namespace
+    Foam::polyMeshTools
+
+Description
+    Collection of static functions operating on polyMesh (mainly checks) so
+    that need access to patch information.
+
+SourceFiles
+    polyMeshTools.C
+
+\*---------------------------------------------------------------------------*/
+#ifndef polyMeshTools_H
+#define polyMeshTools_H
+
+#include "polyMesh.H"
+#include "primitiveMeshTools.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Namespace polyMeshTools Declaration
+\*---------------------------------------------------------------------------*/
+
+class polyMeshTools
+:
+    public primitiveMeshTools
+{
+
+public:
+
+    //- Generate orthogonality field. (1 for fully orthogonal, < 1 for
+    //  non-orthogonal)
+    static tmp<scalarField> faceOrthogonality
+    (
+        const polyMesh& mesh,
+        const vectorField& fAreas,
+        const vectorField& cellCtrs
+    );
+
+//    static tmp<scalarField> faceOrthogonality(const polyMesh& mesh)
+//    {
+//        return faceOrthogonality
+//        (
+//            mesh,
+//            mesh.faceAreas(),
+//            mesh.cellCentres()
+//        );
+//    }
+
+    //- Generate skewness field
+    static tmp<scalarField> faceSkewness
+    (
+        const polyMesh& mesh,
+        const pointField& points,
+        const vectorField& fCtrs,
+        const vectorField& fAreas,
+        const vectorField& cellCtrs
+    );
+
+//    static tmp<scalarField> faceSkewness(const polyMesh& mesh)
+//    {
+//        return faceSkewness
+//        (
+//            mesh,
+//            mesh.points(),
+//            mesh.faceCentres(),
+//            mesh.faceAreas(),
+//            mesh.cellCentres()
+//        );
+//    }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
index 8c647d8454..c4ab139ca3 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,9 +43,6 @@ SourceFiles
     primitiveMeshEdges.C
     primitiveMeshCellCentresAndVols.C
     primitiveMeshFaceCentresAndAreas.C
-    primitiveMeshEdgeVectors.C
-    primitiveMeshCheck.C
-    primitiveMeshCheckMotion.C
     primitiveMeshFindCell.C
 
 \*---------------------------------------------------------------------------*/
@@ -69,6 +66,8 @@ SourceFiles
 namespace Foam
 {
 
+class PackedBoolList;
+
 /*---------------------------------------------------------------------------*\
                       Class primitiveMesh Declaration
 \*---------------------------------------------------------------------------*/
@@ -225,6 +224,28 @@ class primitiveMesh
                 const labelList&
             );
 
+protected:
+
+    // Static data members
+
+        //- Static data to control mesh checking
+
+            //- Cell closedness warning threshold
+            //  set as the fraction of un-closed area to closed area
+            static scalar closedThreshold_;
+
+            //- Aspect ratio warning threshold
+            static scalar aspectThreshold_;
+
+            //- Non-orthogonality warning threshold in deg
+            static scalar nonOrthThreshold_;
+
+            //- Skewness warning threshold
+            static scalar skewThreshold_;
+
+            //- Threshold where faces are considered coplanar
+            static scalar planarCosAngle_;
+
 
         // Geometrical calculations
 
@@ -251,7 +272,7 @@ class primitiveMesh
             void calcEdgeVectors() const;
 
 
-        // Helper functions for mesh checking
+        // Mesh checking
 
             //- Check if all points on face are shared with another face.
             bool checkDuplicateFaces
@@ -270,29 +291,104 @@ class primitiveMesh
                 labelHashSet*
             ) const;
 
+            //- Check boundary for closedness
+            bool checkClosedBoundary
+            (
+                const vectorField&,
+                const bool,
+                const PackedBoolList&
+            ) const;
 
-    // Static data members
+            //- Check cells for closedness
+            bool checkClosedCells
+            (
+                const vectorField& faceAreas,
+                const scalarField& cellVolumes,
+                const bool report,
+                labelHashSet* setPtr,
+                labelHashSet* aspectSetPtr,
+                const Vector<label>& meshD
+            ) const;
 
-        //- Static data to control mesh checking
+            //- Check for negative face areas
+            bool checkFaceAreas
+            (
+                const vectorField& faceAreas,
+                const bool report,
+                const bool detailedReport,
+                labelHashSet* setPtr
+            ) const;
 
-            //- Cell closedness warning threshold
-            //  set as the fraction of un-closed area to closed area
-            static scalar closedThreshold_;
+            //- Check for negative cell volumes
+            bool checkCellVolumes
+            (
+                const scalarField& vols,
+                const bool report,
+                const bool detailedReport,
+                labelHashSet* setPtr
+            ) const;
 
-            //- Aspect ratio warning threshold
-            static scalar aspectThreshold_;
+            //- Check for non-orthogonality
+            bool checkFaceOrthogonality
+            (
+                const vectorField& fAreas,
+                const vectorField& cellCtrs,
+                const bool report,
+                labelHashSet* setPtr
+            ) const;
 
-            //- Non-orthogonality warning threshold in deg
-            static scalar nonOrthThreshold_;
+            //- Check face pyramid volume
+            bool checkFacePyramids
+            (
+                const pointField& points,
+                const vectorField& ctrs,
+                const bool report,
+                const bool detailedReport,
+                const scalar minPyrVol,
+                labelHashSet* setPtr
+            ) const;
 
-            //- Skewness warning threshold
-            static scalar skewThreshold_;
+            //- Check face skewness
+            bool checkFaceSkewness
+            (
+                const pointField& points,
+                const vectorField& fCtrs,
+                const vectorField& fAreas,
+                const vectorField& cellCtrs,
+                const bool report,
+                labelHashSet* setPtr
+            ) const;
 
-            //- Threshold where faces are considered coplanar
-            static scalar planarCosAngle_;
+            //- Check face angles
+            bool checkFaceAngles
+            (
+                const pointField& points,
+                const vectorField& faceAreas,
+                const bool report,
+                const scalar maxDeg,
+                labelHashSet* setPtr
+            ) const;
 
+            //- Check face warpage
+            bool checkFaceFlatness
+            (
+                const pointField& points,
+                const vectorField& faceCentres,
+                const vectorField& faceAreas,
+                const bool report,
+                const scalar warnFlatness,
+                labelHashSet* setPtr
+            ) const;
+
+            //- Check for concave cells by the planes of faces
+            bool checkConcaveCells
+            (
+                const vectorField& fAreas,
+                const pointField& fCentres,
+                const bool report,
+                labelHashSet* setPtr
+            ) const;
 
-protected:
 
         //- Construct null
         primitiveMesh();
@@ -502,29 +598,36 @@ public:
 
             // Topological checks
 
+                //- Check face ordering
+                virtual bool checkUpperTriangular
+                (
+                    const bool report = false,
+                    labelHashSet* setPtr = NULL
+                ) const;
+
                 //- Check cell zip-up
-                bool checkCellsZipUp
+                virtual bool checkCellsZipUp
                 (
                     const bool report = false,
                     labelHashSet* setPtr = NULL
                 ) const;
 
                 //- Check uniqueness of face vertices
-                bool checkFaceVertices
+                virtual bool checkFaceVertices
+                (
+                    const bool report = false,
+                    labelHashSet* setPtr = NULL
+                ) const;
+
+                //- Check for unused points
+                virtual bool checkPoints
                 (
                     const bool report = false,
                     labelHashSet* setPtr = NULL
                 ) const;
 
                 //- Check face-face connectivity
-                bool checkFaceFaces
-                (
-                    const bool report = false,
-                    labelHashSet* setPtr = NULL
-                ) const;
-
-                //- Check face ordering
-                bool checkUpperTriangular
+                virtual bool checkFaceFaces
                 (
                     const bool report = false,
                     labelHashSet* setPtr = NULL
@@ -534,10 +637,11 @@ public:
             // Geometric checks
 
                 //- Check boundary for closedness
-                bool checkClosedBoundary(const bool report = false) const;
+                virtual bool checkClosedBoundary(const bool report = false)
+                const;
 
                 //- Check cells for closedness
-                bool checkClosedCells
+                virtual bool checkClosedCells
                 (
                     const bool report = false,
                     labelHashSet* setPtr = NULL,
@@ -546,28 +650,28 @@ public:
                 ) const;
 
                 //- Check for negative face areas
-                bool checkFaceAreas
+                virtual bool checkFaceAreas
                 (
                     const bool report = false,
                     labelHashSet* setPtr = NULL
                 ) const;
 
                 //- Check for negative cell volumes
-                bool checkCellVolumes
+                virtual bool checkCellVolumes
                 (
                     const bool report = false,
                     labelHashSet* setPtr = NULL
                 ) const;
 
                 //- Check for non-orthogonality
-                bool checkFaceOrthogonality
+                virtual bool checkFaceOrthogonality
                 (
                     const bool report = false,
                     labelHashSet* setPtr = NULL
                 ) const;
 
                 //- Check face pyramid volume
-                bool checkFacePyramids
+                virtual bool checkFacePyramids
                 (
                     const bool report = false,
                     const scalar minPyrVol = -SMALL,
@@ -575,14 +679,14 @@ public:
                 ) const;
 
                 //- Check face skewness
-                bool checkFaceSkewness
+                virtual bool checkFaceSkewness
                 (
                     const bool report = false,
                     labelHashSet* setPtr = NULL
                 ) const;
 
                 //- Check face angles
-                bool checkFaceAngles
+                virtual bool checkFaceAngles
                 (
                     const bool report = false,
                     const scalar maxSin = 10,    // In degrees
@@ -592,31 +696,16 @@ public:
                 //- Check face warpage: decompose face and check ratio between
                 //  magnitude of sum of triangle areas and sum of magnitude of
                 //  triangle areas.
-                bool checkFaceFlatness
+                virtual bool checkFaceFlatness
                 (
                     const bool report,
                     const scalar warnFlatness,  // When to include in set.
                     labelHashSet* setPtr
                 ) const;
 
-                //- Check edge alignment for 1D/2D cases
-                bool checkEdgeAlignment
-                (
-                    const bool report,
-                    const Vector<label>& directions,
-                    labelHashSet* setPtr = NULL
-                ) const;
-
-                //- Check for unused points
-                bool checkPoints
-                (
-                    const bool report = false,
-                    labelHashSet* setPtr = NULL
-                ) const;
-
                 //- Check for point-point-nearness,
                 //  e.g. colocated points which may be part of baffles.
-                bool checkPointNearness
+                virtual bool checkPointNearness
                 (
                     const bool report,
                     const scalar reportDistSqr,
@@ -624,23 +713,15 @@ public:
                 ) const;
 
                 //- Check edge length
-                bool checkEdgeLength
+                virtual bool checkEdgeLength
                 (
                     const bool report,
                     const scalar minLenSqr,
                     labelHashSet* setPtr = NULL
                 ) const;
 
-                //- Check cell determinant
-                bool checkCellDeterminant
-                (
-                    const bool report = false,
-                    labelHashSet* setPtr = NULL,
-                    const Vector<label>& solutionD = Vector<label>::one
-                ) const;
-
                 //- Check for concave cells by the planes of faces
-                bool checkConcaveCells
+                virtual bool checkConcaveCells
                 (
                     const bool report = false,
                     labelHashSet* setPtr = NULL
@@ -649,22 +730,14 @@ public:
 
             //- Check mesh topology for correctness.
             //  Returns false for no error.
-            bool checkTopology(const bool report = false) const;
+            virtual bool checkTopology(const bool report = false) const;
 
             //- Check mesh geometry (& implicitly topology) for correctness.
             //  Returns false for no error.
-            bool checkGeometry(const bool report = false) const;
+            virtual bool checkGeometry(const bool report = false) const;
 
             //- Check mesh for correctness. Returns false for no error.
-            bool checkMesh(const bool report = false) const;
-
-            //- Check mesh motion for correctness given motion points
-            bool checkMeshMotion
-            (
-                const pointField& newPoints,
-                const bool report = false
-            ) const;
-
+            virtual bool checkMesh(const bool report = false) const;
 
             //- Set the closedness ratio warning threshold
             static scalar setClosedThreshold(const scalar);
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
index 76fa1aa5d1..d805608e44 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
@@ -29,6 +29,7 @@ License
 #include "unitConversion.H"
 #include "SortableList.H"
 #include "EdgeMap.H"
+#include "primitiveMeshTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -41,7 +42,12 @@ Foam::scalar Foam::primitiveMesh::planarCosAngle_   = 1.0e-6;
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-bool Foam::primitiveMesh::checkClosedBoundary(const bool report) const
+bool Foam::primitiveMesh::checkClosedBoundary
+(
+    const vectorField& areas,
+    const bool report,
+    const PackedBoolList& internalOrCoupledFaces
+) const
 {
     if (debug)
     {
@@ -56,12 +62,13 @@ bool Foam::primitiveMesh::checkClosedBoundary(const bool report) const
     vector sumClosed(vector::zero);
     scalar sumMagClosedBoundary = 0;
 
-    const vectorField& areas = faceAreas();
-
     for (label faceI = nInternalFaces(); faceI < areas.size(); faceI++)
     {
-        sumClosed += areas[faceI];
-        sumMagClosedBoundary += mag(areas[faceI]);
+        if (internalOrCoupledFaces.size() && !internalOrCoupledFaces[faceI])
+        {
+            sumClosed += areas[faceI];
+            sumMagClosedBoundary += mag(areas[faceI]);
+        }
     }
 
     reduce(sumClosed, sumOp<vector>());
@@ -95,6 +102,8 @@ bool Foam::primitiveMesh::checkClosedBoundary(const bool report) const
 
 bool Foam::primitiveMesh::checkClosedCells
 (
+    const vectorField& faceAreas,
+    const scalarField& cellVolumes,
     const bool report,
     labelHashSet* setPtr,
     labelHashSet* aspectSetPtr,
@@ -140,66 +149,28 @@ bool Foam::primitiveMesh::checkClosedCells
         return true;
     }
 
-    // Loop through cell faces and sum up the face area vectors for each cell.
-    // This should be zero in all vector components
 
-    vectorField sumClosed(nCells(), vector::zero);
-    vectorField sumMagClosed(nCells(), vector::zero);
-
-    const labelList& own = faceOwner();
-    const labelList& nei = faceNeighbour();
-    const vectorField& areas = faceAreas();
-
-    forAll(own, faceI)
-    {
-        // Add to owner
-        sumClosed[own[faceI]] += areas[faceI];
-        sumMagClosed[own[faceI]] += cmptMag(areas[faceI]);
-    }
-
-    forAll(nei, faceI)
-    {
-        // Subtract from neighbour
-        sumClosed[nei[faceI]] -= areas[faceI];
-        sumMagClosed[nei[faceI]] += cmptMag(areas[faceI]);
-    }
+    scalarField openness;
+    scalarField aspectRatio;
+    primitiveMeshTools::cellClosedness
+    (
+        *this,
+        meshD,
+        faceAreas,
+        cellVolumes,
+        openness,
+        aspectRatio
+    );
 
     label nOpen = 0;
-    scalar maxOpennessCell = 0;
-
+    scalar maxOpennessCell = max(openness);
     label nAspect = 0;
-    scalar maxAspectRatio = 0;
-
-    const scalarField& vols = cellVolumes();
-
-    label nDims = 0;
-    for (direction dir = 0; dir < vector::nComponents; dir++)
-    {
-        if (meshD[dir] == 1)
-        {
-            nDims++;
-        }
-    }
-
+    scalar maxAspectRatio = max(aspectRatio);
 
     // Check the sums
-    forAll(sumClosed, cellI)
+    forAll(openness, cellI)
     {
-        scalar maxOpenness = 0;
-
-        for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
-        {
-            maxOpenness = max
-            (
-                maxOpenness,
-                mag(sumClosed[cellI][cmpt])
-               /(sumMagClosed[cellI][cmpt] + VSMALL)
-            );
-        }
-
-        maxOpennessCell = max(maxOpennessCell, maxOpenness);
-
-        if (maxOpenness > closedThreshold_)
+        if (openness[cellI] > closedThreshold_)
         {
             if (setPtr)
             {
@@ -209,32 +180,7 @@ bool Foam::primitiveMesh::checkClosedCells
             nOpen++;
         }
 
-        // Calculate the aspect ration as the maximum of Cartesian component
-        // aspect ratio to the total area hydraulic area aspect ratio
-        scalar minCmpt = VGREAT;
-        scalar maxCmpt = -VGREAT;
-        for (direction dir = 0; dir < vector::nComponents; dir++)
-        {
-            if (meshD[dir] == 1)
-            {
-                minCmpt = min(minCmpt, sumMagClosed[cellI][dir]);
-                maxCmpt = max(maxCmpt, sumMagClosed[cellI][dir]);
-            }
-        }
-
-        scalar aspectRatio = maxCmpt/(minCmpt + VSMALL);
-        if (nDims == 3)
-        {
-            aspectRatio = max
-            (
-                aspectRatio,
-                1.0/6.0*cmptSum(sumMagClosed[cellI])/pow(vols[cellI], 2.0/3.0)
-            );
-        }
-
-        maxAspectRatio = max(maxAspectRatio, aspectRatio);
-
-        if (aspectRatio > aspectThreshold_)
+        if (aspectRatio[cellI] > aspectThreshold_)
         {
             if (aspectSetPtr)
             {
@@ -290,7 +236,9 @@ bool Foam::primitiveMesh::checkClosedCells
 
 bool Foam::primitiveMesh::checkFaceAreas
 (
+    const vectorField& faceAreas,
     const bool report,
+    const bool detailedReport,
     labelHashSet* setPtr
 ) const
 {
@@ -301,7 +249,7 @@ bool Foam::primitiveMesh::checkFaceAreas
             << "checking face area magnitudes" << endl;
     }
 
-    const scalarField magFaceAreas(mag(faceAreas()));
+    const scalarField magFaceAreas(mag(faceAreas));
 
     scalar minArea = GREAT;
     scalar maxArea = -GREAT;
@@ -314,6 +262,25 @@ bool Foam::primitiveMesh::checkFaceAreas
             {
                 setPtr->insert(faceI);
             }
+            if (detailedReport)
+            {
+                if (isInternalFace(faceI))
+                {
+                    Pout<< "Zero or negative face area detected for "
+                        << "internal face "<< faceI << " between cells "
+                        << faceOwner()[faceI] << " and "
+                        << faceNeighbour()[faceI]
+                        << ".  Face area magnitude = " << magFaceAreas[faceI]
+                        << endl;
+                }
+                else
+                {
+                    Pout<< "Zero or negative face area detected for "
+                        << "boundary face " << faceI << " next to cell "
+                        << faceOwner()[faceI] << ".  Face area magnitude = "
+                        << magFaceAreas[faceI] << endl;
+                }
+            }
         }
 
         minArea = min(minArea, magFaceAreas[faceI]);
@@ -349,7 +316,9 @@ bool Foam::primitiveMesh::checkFaceAreas
 
 bool Foam::primitiveMesh::checkCellVolumes
 (
+    const scalarField& vols,
     const bool report,
+    const bool detailedReport,
     labelHashSet* setPtr
 ) const
 {
@@ -360,8 +329,6 @@ bool Foam::primitiveMesh::checkCellVolumes
             << "checking cell volumes" << endl;
     }
 
-    const scalarField& vols = cellVolumes();
-
     scalar minVolume = GREAT;
     scalar maxVolume = -GREAT;
 
@@ -375,6 +342,11 @@ bool Foam::primitiveMesh::checkCellVolumes
             {
                 setPtr->insert(cellI);
             }
+            if (detailedReport)
+            {
+                Pout<< "Zero or negative cell volume detected for cell "
+                    << cellI << ".  Volume = " << vols[cellI] << endl;
+            }
 
             nNegVolCells++;
         }
@@ -416,6 +388,8 @@ bool Foam::primitiveMesh::checkCellVolumes
 
 bool Foam::primitiveMesh::checkFaceOrthogonality
 (
+    const vectorField& fAreas,
+    const vectorField& cellCtrs,
     const bool report,
     labelHashSet* setPtr
 ) const
@@ -427,35 +401,33 @@ bool Foam::primitiveMesh::checkFaceOrthogonality
             << "checking mesh non-orthogonality" << endl;
     }
 
-    // for all internal faces check that the d dot S product is positive
-    const vectorField& centres = cellCentres();
-    const vectorField& areas = faceAreas();
 
-    const labelList& own = faceOwner();
-    const labelList& nei = faceNeighbour();
+    tmp<scalarField> tortho = primitiveMeshTools::faceOrthogonality
+    (
+        *this,
+        fAreas,
+        cellCtrs
+    );
+    const scalarField& ortho = tortho();
 
     // Severe nonorthogonality threshold
     const scalar severeNonorthogonalityThreshold =
         ::cos(degToRad(nonOrthThreshold_));
 
-    scalar minDDotS = GREAT;
+    scalar minDDotS = min(ortho);
 
-    scalar sumDDotS = 0;
+    scalar sumDDotS = sum(ortho);
 
     label severeNonOrth = 0;
 
     label errorNonOrth = 0;
 
-    forAll(nei, faceI)
+
+    forAll(ortho, faceI)
     {
-        vector d = centres[nei[faceI]] - centres[own[faceI]];
-        const vector& s = areas[faceI];
-
-        scalar dDotS = (d & s)/(mag(d)*mag(s) + VSMALL);
-
-        if (dDotS < severeNonorthogonalityThreshold)
+        if (ortho[faceI] < severeNonorthogonalityThreshold)
         {
-            if (dDotS > SMALL)
+            if (ortho[faceI] > SMALL)
             {
                 if (setPtr)
                 {
@@ -474,13 +446,6 @@ bool Foam::primitiveMesh::checkFaceOrthogonality
                 errorNonOrth++;
             }
         }
-
-        if (dDotS < minDDotS)
-        {
-            minDDotS = dDotS;
-        }
-
-        sumDDotS += dDotS;
     }
 
     reduce(minDDotS, minOp<scalar>());
@@ -490,7 +455,7 @@ bool Foam::primitiveMesh::checkFaceOrthogonality
 
     if (debug || report)
     {
-        label neiSize = nei.size();
+        label neiSize = ortho.size();
         reduce(neiSize, sumOp<label>());
 
         if (neiSize > 0)
@@ -535,7 +500,10 @@ bool Foam::primitiveMesh::checkFaceOrthogonality
 
 bool Foam::primitiveMesh::checkFacePyramids
 (
+    const pointField& points,
+    const vectorField& ctrs,
     const bool report,
+    const bool detailedReport,
     const scalar minPyrVol,
     labelHashSet* setPtr
 ) const
@@ -547,46 +515,63 @@ bool Foam::primitiveMesh::checkFacePyramids
             << "checking face orientation" << endl;
     }
 
-    // check whether face area vector points to the cell with higher label
-    const vectorField& ctrs = cellCentres();
-
     const labelList& own = faceOwner();
     const labelList& nei = faceNeighbour();
-
     const faceList& f = faces();
 
-    const pointField& p = points();
+
+    scalarField ownPyrVol;
+    scalarField neiPyrVol;
+    primitiveMeshTools::facePyramidVolume
+    (
+        *this,
+        points,
+        ctrs,
+        ownPyrVol,
+        neiPyrVol
+    );
+
 
     label nErrorPyrs = 0;
 
-    forAll(f, faceI)
+    forAll(ownPyrVol, faceI)
     {
-        // Create the owner pyramid - it will have negative volume
-        scalar pyrVol = pyramidPointFaceRef(f[faceI], ctrs[own[faceI]]).mag(p);
-
-        if (pyrVol > -minPyrVol)
+        if (ownPyrVol[faceI] < minPyrVol)
         {
             if (setPtr)
             {
                 setPtr->insert(faceI);
             }
+            if (detailedReport)
+            {
+                Pout<< "Negative pyramid volume: " << ownPyrVol[faceI]
+                    << " for face " << faceI << " " << f[faceI]
+                    << "  and owner cell: " << own[faceI] << endl
+                    << "Owner cell vertex labels: "
+                    << cells()[own[faceI]].labels(faces())
+                    << endl;
+            }
 
             nErrorPyrs++;
         }
 
         if (isInternalFace(faceI))
         {
-            // Create the neighbour pyramid - it will have positive volume
-            scalar pyrVol =
-                pyramidPointFaceRef(f[faceI], ctrs[nei[faceI]]).mag(p);
-
-            if (pyrVol < minPyrVol)
+            if (neiPyrVol[faceI] < minPyrVol)
             {
                 if (setPtr)
                 {
                     setPtr->insert(faceI);
                 }
-
+                if (detailedReport)
+                {
+                    Pout<< "Negative pyramid volume: " << neiPyrVol[faceI]
+                        << " for face " << faceI << " " << f[faceI]
+                        << "  and neighbour cell: " << nei[faceI] << nl
+                        << "Neighbour cell vertex labels: "
+                        << cells()[nei[faceI]].labels(faces())
+                        << endl;
+                }
                 nErrorPyrs++;
             }
         }
@@ -619,6 +604,10 @@ bool Foam::primitiveMesh::checkFacePyramids
 
 bool Foam::primitiveMesh::checkFaceSkewness
 (
+    const pointField& points,
+    const vectorField& fCtrs,
+    const vectorField& fAreas,
+    const vectorField& cellCtrs,
     const bool report,
     labelHashSet* setPtr
 ) const
@@ -633,44 +622,24 @@ bool Foam::primitiveMesh::checkFaceSkewness
     // Warn if the skew correction vector is more than skewWarning times
     // larger than the face area vector
 
-    const pointField& p = points();
-    const faceList& fcs = faces();
+    tmp<scalarField> tskewness = primitiveMeshTools::faceSkewness
+    (
+        *this,
+        points,
+        fCtrs,
+        fAreas,
+        cellCtrs
+    );
+    const scalarField& skewness = tskewness();
 
-    const labelList& own = faceOwner();
-    const labelList& nei = faceNeighbour();
-    const vectorField& cellCtrs = cellCentres();
-    const vectorField& faceCtrs = faceCentres();
-    const vectorField& fAreas = faceAreas();
-
-    scalar maxSkew = 0;
+    scalar maxSkew = max(skewness);
     label nWarnSkew = 0;
 
-    forAll(nei, faceI)
+    forAll(skewness, faceI)
     {
-        vector Cpf = faceCtrs[faceI] - cellCtrs[own[faceI]];
-        vector d = cellCtrs[nei[faceI]] - cellCtrs[own[faceI]];
-
-        // Skewness vector
-        vector sv =
-            Cpf - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + SMALL))*d;
-        vector svHat = sv/(mag(sv) + VSMALL);
-
-        // Normalisation distance calculated as the approximate distance
-        // from the face centre to the edge of the face in the direction of
-        // the skewness
-        scalar fd = 0.2*mag(d) + VSMALL;
-        const face& f = fcs[faceI];
-        forAll(f, pi)
-        {
-            fd = max(fd, mag(svHat & (p[f[pi]] - faceCtrs[faceI])));
-        }
-
-        // Normalised skewness
-        scalar skewness = mag(sv)/fd;
-
         // Check if the skewness vector is greater than the PN vector.
         // This does not cause trouble but is a good indication of a poor mesh.
-        if (skewness > skewThreshold_)
+        if (skewness[faceI] > skewThreshold_)
         {
             if (setPtr)
             {
@@ -679,63 +648,8 @@ bool Foam::primitiveMesh::checkFaceSkewness
 
             nWarnSkew++;
         }
-
-        if (skewness > maxSkew)
-        {
-            maxSkew = skewness;
-        }
     }
 
-
-    // Boundary faces: consider them to have only skewness error.
-    // (i.e. treat as if mirror cell on other side)
-
-    for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
-    {
-        vector Cpf = faceCtrs[faceI] - cellCtrs[own[faceI]];
-
-        vector normal = fAreas[faceI];
-        normal /= mag(normal) + VSMALL;
-        vector d = normal*(normal & Cpf);
-
-
-        // Skewness vector
-        vector sv =
-            Cpf - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d;
-        vector svHat = sv/(mag(sv) + VSMALL);
-
-        // Normalisation distance calculated as the approximate distance
-        // from the face centre to the edge of the face in the direction of
-        // the skewness
-        scalar fd = 0.4*mag(d) + VSMALL;
-        const face& f = fcs[faceI];
-        forAll(f, pi)
-        {
-            fd = max(fd, mag(svHat & (p[f[pi]] - faceCtrs[faceI])));
-        }
-
-        // Normalised skewness
-        scalar skewness = mag(sv)/fd;
-
-        // Check if the skewness vector is greater than the PN vector.
-        // This does not cause trouble but is a good indication of a poor mesh.
-        if (skewness > skewThreshold_)
-        {
-            if (setPtr)
-            {
-                setPtr->insert(faceI);
-            }
-
-            nWarnSkew++;
-        }
-
-        if (skewness > maxSkew)
-        {
-            maxSkew = skewness;
-        }
-    }
-
-
     reduce(maxSkew, maxOp<scalar>());
     reduce(nWarnSkew, sumOp<label>());
 
@@ -763,86 +677,14 @@ bool Foam::primitiveMesh::checkFaceSkewness
 }
 
 
-bool Foam::primitiveMesh::checkPoints
-(
-    const bool report,
-    labelHashSet* setPtr
-) const
-{
-    if (debug)
-    {
-        Info<< "bool primitiveMesh::checkPoints"
-            << "(const bool, labelHashSet*) const: "
-            << "checking points" << endl;
-    }
-
-    label nFaceErrors = 0;
-    label nCellErrors = 0;
-
-    const labelListList& pf = pointFaces();
-
-    forAll(pf, pointI)
-    {
-        if (pf[pointI].empty())
-        {
-            if (setPtr)
-            {
-                setPtr->insert(pointI);
-            }
-
-            nFaceErrors++;
-        }
-    }
-
-
-    forAll(pf, pointI)
-    {
-        const labelList& pc = pointCells(pointI);
-
-        if (pc.empty())
-        {
-            if (setPtr)
-            {
-                setPtr->insert(pointI);
-            }
-
-            nCellErrors++;
-        }
-    }
-
-    reduce(nFaceErrors, sumOp<label>());
-    reduce(nCellErrors, sumOp<label>());
-
-    if (nFaceErrors > 0 || nCellErrors > 0)
-    {
-        if (debug || report)
-        {
-            Info<< " ***Unused points found in the mesh, "
-                   "number unused by faces: " << nFaceErrors
-                << " number unused by cells: " << nCellErrors
-                << endl;
-        }
-
-        return true;
-    }
-    else
-    {
-        if (debug || report)
-        {
-            Info<< "    Point usage OK." << endl;
-        }
-
-        return false;
-    }
-}
-
-
 // Check convexity of angles in a face. Allow a slight non-convexity.
 // E.g. maxDeg = 10 allows for angles < 190 (or 10 degrees concavity)
 // (if truly concave and points not visible from face centre the face-pyramid
 //  check in checkMesh will fail)
 bool Foam::primitiveMesh::checkFaceAngles
 (
+    const pointField& points,
+    const vectorField& faceAreas,
     const bool report,
     const scalar maxDeg,
     labelHashSet* setPtr
@@ -867,71 +709,30 @@ bool Foam::primitiveMesh::checkFaceAngles
 
     const scalar maxSin = Foam::sin(degToRad(maxDeg));
 
-    const pointField& p = points();
-    const faceList& fcs = faces();
-    vectorField faceNormals(faceAreas());
-    faceNormals /= mag(faceNormals) + VSMALL;
 
-    scalar maxEdgeSin = 0.0;
+    tmp<scalarField> tfaceAngles = primitiveMeshTools::faceConcavity
+    (
+        maxSin,
+        *this,
+        points,
+        faceAreas
+    );
+    const scalarField& faceAngles = tfaceAngles();
+
+    scalar maxEdgeSin = max(faceAngles);
 
     label nConcave = 0;
 
-    label errorFaceI = -1;
-
-    forAll(fcs, faceI)
+    forAll(faceAngles, faceI)
     {
-        const face& f = fcs[faceI];
-
-        // Get edge from f[0] to f[size-1];
-        vector ePrev(p[f.first()] - p[f.last()]);
-        scalar magEPrev = mag(ePrev);
-        ePrev /= magEPrev + VSMALL;
-
-        forAll(f, fp0)
+        if (faceAngles[faceI] > SMALL)
         {
-            // Get vertex after fp
-            label fp1 = f.fcIndex(fp0);
+            nConcave++;
 
-            // Normalized vector between two consecutive points
-            vector e10(p[f[fp1]] - p[f[fp0]]);
-            scalar magE10 = mag(e10);
-            e10 /= magE10 + VSMALL;
-
-            if (magEPrev > SMALL && magE10 > SMALL)
+            if (setPtr)
             {
-                vector edgeNormal = ePrev ^ e10;
-                scalar magEdgeNormal = mag(edgeNormal);
-
-                if (magEdgeNormal < maxSin)
-                {
-                    // Edges (almost) aligned -> face is ok.
-                }
-                else
-                {
-                    // Check normal
-                    edgeNormal /= magEdgeNormal;
-
-                    if ((edgeNormal & faceNormals[faceI]) < SMALL)
-                    {
-                        if (faceI != errorFaceI)
-                        {
-                            // Count only one error per face.
-                            errorFaceI = faceI;
-                            nConcave++;
-                        }
-
-                        if (setPtr)
-                        {
-                            setPtr->insert(faceI);
-                        }
-
-                        maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
-                    }
-                }
+                setPtr->insert(faceI);
             }
-
-            ePrev = e10;
-            magEPrev = magE10;
         }
     }
 
@@ -965,11 +766,11 @@ bool Foam::primitiveMesh::checkFaceAngles
 }
 
 
-// Check warpage of faces. Is calculated as the difference between areas of
-// individual triangles and the overall area of the face (which ifself is
-// is the average of the areas of the individual triangles).
 bool Foam::primitiveMesh::checkFaceFlatness
 (
+    const pointField& points,
+    const vectorField& faceCentres,
+    const vectorField& faceAreas,
     const bool report,
     const scalar warnFlatness,
     labelHashSet* setPtr
@@ -992,52 +793,34 @@ bool Foam::primitiveMesh::checkFaceFlatness
             << exit(FatalError);
     }
 
-
-    const pointField& p = points();
     const faceList& fcs = faces();
-    const pointField& fctrs = faceCentres();
 
-    // Areas are calculated as the sum of areas. (see
-    // primitiveMeshFaceCentresAndAreas.C)
-    scalarField magAreas(mag(faceAreas()));
+    tmp<scalarField> tfaceFlatness = primitiveMeshTools::faceFlatness
+    (
+        *this,
+        points,
+        faceCentres,
+        faceAreas
+    );
+    const scalarField& faceFlatness = tfaceFlatness();
 
-    label nWarped = 0;
+    scalarField magAreas(mag(faceAreas));
 
     scalar minFlatness = GREAT;
     scalar sumFlatness = 0;
     label nSummed = 0;
+    label nWarped = 0;
 
-    forAll(fcs, faceI)
+    forAll(faceFlatness, faceI)
     {
-        const face& f = fcs[faceI];
-
-        if (f.size() > 3 && magAreas[faceI] > VSMALL)
+        if (fcs[faceI].size() > 3 && magAreas[faceI] > VSMALL)
         {
-            const point& fc = fctrs[faceI];
-
-            // Calculate the sum of magnitude of areas and compare to magnitude
-            // of sum of areas.
-
-            scalar sumA = 0.0;
-
-            forAll(f, fp)
-            {
-                const point& thisPoint = p[f[fp]];
-                const point& nextPoint = p[f.nextLabel(fp)];
-
-                // Triangle around fc.
-                vector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
-                sumA += mag(n);
-            }
-
-            scalar flatness = magAreas[faceI] / (sumA+VSMALL);
-
-            sumFlatness += flatness;
+            sumFlatness += faceFlatness[faceI];
             nSummed++;
 
-            minFlatness = min(minFlatness, flatness);
+            minFlatness = min(minFlatness, faceFlatness[faceI]);
 
-            if (flatness < warnFlatness)
+            if (faceFlatness[faceI] < warnFlatness)
             {
                 nWarped++;
 
@@ -1092,131 +875,103 @@ bool Foam::primitiveMesh::checkFaceFlatness
 }
 
 
-// Check 1D/2Dness of edges. Gets passed the non-empty directions and
-// checks all edges in the mesh whether they:
-// - have no component in a non-empty direction or
-// - are only in a singe non-empty direction.
-// Empty direction info is passed in as a vector of labels (synchronised)
-// which are 1 if the direction is non-empty, 0 if it is.
-bool Foam::primitiveMesh::checkEdgeAlignment
+bool Foam::primitiveMesh::checkConcaveCells
 (
+    const vectorField& fAreas,
+    const pointField& fCentres,
     const bool report,
-    const Vector<label>& directions,
     labelHashSet* setPtr
 ) const
 {
     if (debug)
     {
-        Info<< "bool primitiveMesh::checkEdgeAlignment("
-            << "const bool, const Vector<label>&, labelHashSet*) const: "
-            << "checking edge alignment" << endl;
+        Info<< "bool primitiveMesh::checkConcaveCells(const bool"
+            << ", labelHashSet*) const: "
+            << "checking for concave cells" << endl;
     }
 
-    label nDirs = 0;
-    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+    const cellList& c = cells();
+    const labelList& fOwner = faceOwner();
+
+    label nConcaveCells = 0;
+
+    forAll(c, cellI)
     {
-        if (directions[cmpt] == 1)
+        const cell& cFaces = c[cellI];
+
+        bool concave = false;
+
+        forAll(cFaces, i)
         {
-            nDirs++;
-        }
-        else if (directions[cmpt] != 0)
-        {
-            FatalErrorIn
-            (
-                "primitiveMesh::checkEdgeAlignment"
-                "(const bool, const Vector<label>&, labelHashSet*)"
-            )   << "directions should contain 0 or 1 but is now " << directions
-                << exit(FatalError);
-        }
-    }
-
-    if (nDirs == vector::nComponents)
-    {
-        return false;
-    }
-
-
-
-    const pointField& p = points();
-    const faceList& fcs = faces();
-
-    EdgeMap<label> edgesInError;
-
-    forAll(fcs, faceI)
-    {
-        const face& f = fcs[faceI];
-
-        forAll(f, fp)
-        {
-            label p0 = f[fp];
-            label p1 = f.nextLabel(fp);
-            if (p0 < p1)
+            if (concave)
             {
-                vector d(p[p1]-p[p0]);
-                scalar magD = mag(d);
+                break;
+            }
 
-                if (magD > ROOTVSMALL)
+            label fI = cFaces[i];
+
+            const point& fC = fCentres[fI];
+
+            vector fN = fAreas[fI];
+
+            fN /= max(mag(fN), VSMALL);
+
+            // Flip normal if required so that it is always pointing out of
+            // the cell
+            if (fOwner[fI] != cellI)
+            {
+                fN *= -1;
+            }
+
+            // Is the centre of any other face of the cell on the
+            // wrong side of the plane of this face?
+
+            forAll(cFaces, j)
+            {
+                if (j != i)
                 {
-                    d /= magD;
+                    label fJ = cFaces[j];
 
-                    // Check how many empty directions are used by the edge.
-                    label nEmptyDirs = 0;
-                    label nNonEmptyDirs = 0;
-                    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
-                    {
-                        if (mag(d[cmpt]) > 1e-6)
-                        {
-                            if (directions[cmpt] == 0)
-                            {
-                                nEmptyDirs++;
-                            }
-                            else
-                            {
-                                nNonEmptyDirs++;
-                            }
-                        }
-                    }
+                    const point& pt = fCentres[fJ];
 
-                    if (nEmptyDirs == 0)
+                    // If the cell is concave, the point will be on the
+                    // positive normal side of the plane of f, defined by
+                    // its centre and normal, and the angle between (pt -
+                    // fC) and fN will be less than 90 degrees, so the dot
+                    // product will be positive.
+
+                    vector pC = (pt - fC);
+
+                    pC /= max(mag(pC), VSMALL);
+
+                    if ((pC & fN) > -planarCosAngle_)
                     {
-                        // Purely in ok directions.
-                    }
-                    else if (nEmptyDirs == 1)
-                    {
-                        // Ok if purely in empty directions.
-                        if (nNonEmptyDirs > 0)
+                        // Concave or planar face
+
+                        concave = true;
+
+                        if (setPtr)
                         {
-                            edgesInError.insert(edge(p0, p1), faceI);
+                            setPtr->insert(cellI);
                         }
-                    }
-                    else if (nEmptyDirs > 1)
-                    {
-                        // Always an error
-                        edgesInError.insert(edge(p0, p1), faceI);
+
+                        nConcaveCells++;
+
+                        break;
                     }
                 }
             }
         }
     }
 
-    label nErrorEdges = returnReduce(edgesInError.size(), sumOp<label>());
+    reduce(nConcaveCells, sumOp<label>());
 
-    if (nErrorEdges > 0)
+    if (nConcaveCells > 0)
     {
         if (debug || report)
         {
-            Info<< " ***Number of edges not aligned with or perpendicular to "
-                << "non-empty directions: " << nErrorEdges << endl;
-        }
-
-        if (setPtr)
-        {
-            setPtr->resize(2*edgesInError.size());
-            forAllConstIter(EdgeMap<label>, edgesInError, iter)
-            {
-                setPtr->insert(iter.key()[0]);
-                setPtr->insert(iter.key()[1]);
-            }
+            Info<< " ***Concave cells (using face planes) found,"
+                << " number of cells: " << nConcaveCells << endl;
         }
 
         return true;
@@ -1225,14 +980,18 @@ bool Foam::primitiveMesh::checkEdgeAlignment
     {
         if (debug || report)
         {
-            Info<< "    All edges aligned with or perpendicular to "
-                << "non-empty directions." << endl;
+            Info<< "    Concave cell check OK." << endl;
         }
+
         return false;
     }
+
+    return false;
 }
 
 
+// Topological tests
+
 bool Foam::primitiveMesh::checkUpperTriangular
 (
     const bool report,
@@ -1574,6 +1333,80 @@ bool Foam::primitiveMesh::checkFaceVertices
 }
 
 
+bool Foam::primitiveMesh::checkPoints
+(
+    const bool report,
+    labelHashSet* setPtr
+) const
+{
+    if (debug)
+    {
+        Info<< "bool primitiveMesh::checkPoints"
+            << "(const bool, labelHashSet*) const: "
+            << "checking points" << endl;
+    }
+
+    label nFaceErrors = 0;
+    label nCellErrors = 0;
+
+    const labelListList& pf = pointFaces();
+
+    forAll(pf, pointI)
+    {
+        if (pf[pointI].empty())
+        {
+            if (setPtr)
+            {
+                setPtr->insert(pointI);
+            }
+
+            nFaceErrors++;
+        }
+    }
+
+
+    forAll(pf, pointI)
+    {
+        const labelList& pc = pointCells(pointI);
+
+        if (pc.empty())
+        {
+            if (setPtr)
+            {
+                setPtr->insert(pointI);
+            }
+
+            nCellErrors++;
+        }
+    }
+
+    reduce(nFaceErrors, sumOp<label>());
+    reduce(nCellErrors, sumOp<label>());
+
+    if (nFaceErrors > 0 || nCellErrors > 0)
+    {
+        if (debug || report)
+        {
+            Info<< " ***Unused points found in the mesh, "
+                   "number unused by faces: " << nFaceErrors
+                << " number unused by cells: " << nCellErrors
+                << endl;
+        }
+
+        return true;
+    }
+    else
+    {
+        if (debug || report)
+        {
+            Info<< "    Point usage OK." << endl;
+        }
+
+        return false;
+    }
+}
+
+
 // Check if all points on face are shared between faces.
 bool Foam::primitiveMesh::checkDuplicateFaces
 (
@@ -1891,168 +1724,153 @@ bool Foam::primitiveMesh::checkFaceFaces
 }
 
 
-// Checks cells with 1 or less internal faces. Give numerical problems.
-bool Foam::primitiveMesh::checkCellDeterminant
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::primitiveMesh::checkClosedBoundary(const bool report) const
+{
+    return checkClosedBoundary(faceAreas(), report, PackedBoolList(0));
+}
+
+
+bool Foam::primitiveMesh::checkClosedCells
 (
-    const bool report,    // report,
-    labelHashSet* setPtr, // setPtr
-    const Vector<label>& meshD
+    const bool report,
+    labelHashSet* setPtr,
+    labelHashSet* aspectSetPtr,
+     const Vector<label>& solutionD
 ) const
 {
-    if (debug)
-    {
-        Info<< "bool primitiveMesh::checkCellDeterminant(const bool"
-            << ", labelHashSet*) const: "
-            << "checking for under-determined cells" << endl;
-    }
-
-    // Determine number of dimensions and (for 2D) missing dimension
-    label nDims = 0;
-    label twoD = -1;
-    for (direction dir = 0; dir < vector::nComponents; dir++)
-    {
-        if (meshD[dir] == 1)
-        {
-            nDims++;
-        }
-        else
-        {
-            twoD = dir;
-        }
-    }
+    return checkClosedCells
+    (
+        faceAreas(),
+        cellVolumes(),
+        report,
+        setPtr,
+        aspectSetPtr,
+        solutionD
+    );
+}
 
 
-    const cellList& c = cells();
+bool Foam::primitiveMesh::checkFaceAreas
+(
+    const bool report,
+    labelHashSet* setPtr
+) const
+{
+    return checkFaceAreas
+    (
+        faceAreas(),
+        report,
+        false,  // detailedReport,
+        setPtr
+    );
+}
 
-    label nErrorCells = 0;
 
-    scalar minDet = GREAT;
-    scalar sumDet = 0;
-    label nSummed = 0;
+bool Foam::primitiveMesh::checkCellVolumes
+(
+    const bool report,
+    labelHashSet* setPtr
+) const
+{
+    return checkCellVolumes
+    (
+        cellVolumes(),
+        report,
+        false,  // detailedReport,
+        setPtr
+    );
+}
 
-    if (nDims == 1)
-    {
-        minDet = 1;
-        sumDet = c.size()*minDet;
-        nSummed = c.size();
-    }
-    else
-    {
-        forAll (c, cellI)
-        {
-            const labelList& curFaces = c[cellI];
 
-            // Calculate local normalization factor
-            scalar avgArea = 0;
+bool Foam::primitiveMesh::checkFaceOrthogonality
+(
+    const bool report,
+    labelHashSet* setPtr
+) const
+{
+    return checkFaceOrthogonality
+    (
+        faceAreas(),
+        cellCentres(),
+        report,
+        setPtr
+    );
+}
 
-            label nInternalFaces = 0;
 
-            forAll(curFaces, i)
-            {
-                if (isInternalFace(curFaces[i]))
-                {
-                    avgArea += mag(faceAreas()[curFaces[i]]);
+bool Foam::primitiveMesh::checkFacePyramids
+(
+    const bool report,
+    const scalar minPyrVol,
+    labelHashSet* setPtr
+) const
+{
+    return checkFacePyramids
+    (
+        points(),
+        cellCentres(),
+        report,
+        false,  // detailedReport,
+        minPyrVol,
+        setPtr
+    );
+}
 
-                    nInternalFaces++;
-                }
-            }
 
-            if (nInternalFaces == 0)
-            {
-                if (setPtr)
-                {
-                    setPtr->insert(cellI);
-                }
+bool Foam::primitiveMesh::checkFaceSkewness
+(
+    const bool report,
+    labelHashSet* setPtr
+) const
+{
+    return checkFaceSkewness
+    (
+        points(),
+        faceCentres(),
+        faceAreas(),
+        cellCentres(),
+        report,
+        setPtr
+    );
+}
 
-                nErrorCells++;
-            }
-            else
-            {
-                avgArea /= nInternalFaces;
 
-                symmTensor areaTensor(symmTensor::zero);
+bool Foam::primitiveMesh::checkFaceAngles
+(
+    const bool report,
+    const scalar maxDeg,
+    labelHashSet* setPtr
+) const
+{
+    return checkFaceAngles
+    (
+        points(),
+        faceAreas(),
+        report,
+        maxDeg,
+        setPtr
+    );
+}
 
-                forAll(curFaces, i)
-                {
-                    if (isInternalFace(curFaces[i]))
-                    {
-                        areaTensor += sqr(faceAreas()[curFaces[i]]/avgArea);
-                    }
-                }
 
-                if (nDims == 2)
-                {
-                    // Add the missing eigenvector (such that it does not
-                    // affect the determinant)
-                    if (twoD == 0)
-                    {
-                        areaTensor.xx() = 1;
-                    }
-                    else if (twoD == 1)
-                    {
-                        areaTensor.yy() = 1;
-                    }
-                    else
-                    {
-                        areaTensor.zz() = 1;
-                    }
-                }
-
-                scalar determinant = mag(det(areaTensor));
-
-                minDet = min(determinant, minDet);
-                sumDet += determinant;
-                nSummed++;
-
-                if (determinant < 1e-3)
-                {
-                    if (setPtr)
-                    {
-                        setPtr->insert(cellI);
-                    }
-
-                    nErrorCells++;
-                }
-            }
-        }
-    }
-
-    reduce(nErrorCells, sumOp<label>());
-    reduce(minDet, minOp<scalar>());
-    reduce(sumDet, sumOp<scalar>());
-    reduce(nSummed, sumOp<label>());
-
-    if (debug || report)
-    {
-        if (nSummed > 0)
-        {
-            Info<< "    Cell determinant (wellposedness) : minimum: " << minDet
-                << " average: " << sumDet/nSummed
-                << endl;
-        }
-    }
-
-    if (nErrorCells > 0)
-    {
-        if (debug || report)
-        {
-            Info<< " ***Cells with small determinant found, number of cells: "
-                << nErrorCells << endl;
-        }
-
-        return true;
-    }
-    else
-    {
-        if (debug || report)
-        {
-            Info<< "    Cell determinant check OK." << endl;
-        }
-
-        return false;
-    }
-
-    return false;
+bool Foam::primitiveMesh::checkFaceFlatness
+(
+    const bool report,
+    const scalar warnFlatness,
+    labelHashSet* setPtr
+) const
+{
+    return checkFaceFlatness
+    (
+        points(),
+        faceCentres(),
+        faceAreas(),
+        report,
+        warnFlatness,
+        setPtr
+    );
 }
 
 
@@ -2062,117 +1880,16 @@ bool Foam::primitiveMesh::checkConcaveCells
     labelHashSet* setPtr
 ) const
 {
-    if (debug)
-    {
-        Info<< "bool primitiveMesh::checkConcaveCells(const bool"
-            << ", labelHashSet*) const: "
-            << "checking for concave cells" << endl;
-    }
-
-    const cellList& c = cells();
-    const labelList& fOwner = faceOwner();
-    const vectorField& fAreas = faceAreas();
-    const pointField& fCentres = faceCentres();
-
-    label nConcaveCells = 0;
-
-    forAll(c, cellI)
-    {
-        const cell& cFaces = c[cellI];
-
-        bool concave = false;
-
-        forAll(cFaces, i)
-        {
-            if (concave)
-            {
-                break;
-            }
-
-            label fI = cFaces[i];
-
-            const point& fC = fCentres[fI];
-
-            vector fN = fAreas[fI];
-
-            fN /= max(mag(fN), VSMALL);
-
-            // Flip normal if required so that it is always pointing out of
-            // the cell
-            if (fOwner[fI] != cellI)
-            {
-                fN *= -1;
-            }
-
-            // Is the centre of any other face of the cell on the
-            // wrong side of the plane of this face?
-
-            forAll(cFaces, j)
-            {
-                if (j != i)
-                {
-                    label fJ = cFaces[j];
-
-                    const point& pt = fCentres[fJ];
-
-                    // If the cell is concave, the point will be on the
-                    // positive normal side of the plane of f, defined by
-                    // its centre and normal, and the angle between (pt -
-                    // fC) and fN will be less than 90 degrees, so the dot
-                    // product will be positive.
-
-                    vector pC = (pt - fC);
-
-                    pC /= max(mag(pC), VSMALL);
-
-                    if ((pC & fN) > -planarCosAngle_)
-                    {
-                        // Concave or planar face
-
-                        concave = true;
-
-                        if (setPtr)
-                        {
-                            setPtr->insert(cellI);
-                        }
-
-                        nConcaveCells++;
-
-                        break;
-                    }
-                }
-            }
-        }
-    }
-
-    reduce(nConcaveCells, sumOp<label>());
-
-    if (nConcaveCells > 0)
-    {
-        if (debug || report)
-        {
-            Info<< " ***Concave cells (using face planes) found,"
-                << " number of cells: " << nConcaveCells << endl;
-        }
-
-        return true;
-    }
-    else
-    {
-        if (debug || report)
-        {
-            Info<< "    Concave cell check OK." << endl;
-        }
-
-        return false;
-    }
-
-    return false;
+    return checkConcaveCells
+    (
+        faceAreas(),
+        faceCentres(),
+        report,
+        setPtr
+    );
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
 bool Foam::primitiveMesh::checkTopology(const bool report) const
 {
     label noFailedChecks = 0;
@@ -2182,7 +1899,6 @@ bool Foam::primitiveMesh::checkTopology(const bool report) const
     if (checkCellsZipUp(report)) noFailedChecks++;
     if (checkFaceVertices(report)) noFailedChecks++;
     if (checkFaceFaces(report)) noFailedChecks++;
-    //if (checkCellDeterminant(report)) noFailedChecks++;
 
     if (noFailedChecks == 0)
     {
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheckMotion.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheckMotion.C
deleted file mode 100644
index 4b0d2da63a..0000000000
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheckMotion.C
+++ /dev/null
@@ -1,266 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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/>.
-
-Description
-    Given a set of points, find out if the mesh resulting from point motion will
-    be valid without actually changing the mesh.
-
-\*---------------------------------------------------------------------------*/
-
-#include "primitiveMesh.H"
-#include "pyramidPointFaceRef.H"
-#include "cell.H"
-#include "unitConversion.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-bool Foam::primitiveMesh::checkMeshMotion
-(
-    const pointField& newPoints,
-    const bool report
-) const
-{
-    if (debug || report)
-    {
-        Pout<< "bool primitiveMesh::checkMeshMotion("
-            << "const pointField& newPoints, const bool report) const: "
-            << "checking mesh motion" << endl;
-    }
-
-    bool error = false;
-
-    const faceList& f = faces();
-
-    const labelList& own = faceOwner();
-    const labelList& nei = faceNeighbour();
-
-    vectorField fCtrs(nFaces());
-    vectorField fAreas(nFaces());
-
-    makeFaceCentresAndAreas(newPoints, fCtrs, fAreas);
-
-    // Check cell volumes and calculate new cell centres
-    vectorField cellCtrs(nCells());
-    scalarField cellVols(nCells());
-
-    makeCellCentresAndVols(fCtrs, fAreas, cellCtrs, cellVols);
-
-    scalar minVolume = GREAT;
-    label nNegVols = 0;
-
-    forAll(cellVols, cellI)
-    {
-        if (cellVols[cellI] < VSMALL)
-        {
-            if (debug || report)
-            {
-                Pout<< "Zero or negative cell volume detected for cell "
-                    << cellI << ".  Volume = " << cellVols[cellI] << endl;
-            }
-
-            nNegVols++;
-        }
-
-        minVolume = min(minVolume, cellVols[cellI]);
-    }
-
-    if (nNegVols > 0)
-    {
-        error = true;
-
-        Pout<< "Zero or negative cell volume in mesh motion in " << nNegVols
-            << " cells.  Min volume: " << minVolume << endl;
-    }
-    else
-    {
-        if (debug || report)
-        {
-            Pout<< "Min volume = " << minVolume
-                << ".  Total volume = " << sum(cellVols)
-                << ".  Cell volumes OK." << endl;
-        }
-    }
-
-    // Check face areas
-
-    scalar minArea = GREAT;
-    label nNegAreas = 0;
-    label nPyrErrors = 0;
-    label nDotProductErrors = 0;
-
-    forAll(f, faceI)
-    {
-        const scalar a = Foam::mag(fAreas[faceI]);
-
-        if (a < VSMALL)
-        {
-            if (debug || report)
-            {
-                if (isInternalFace(faceI))
-                {
-                    Pout<< "Zero or negative face area detected for "
-                        << "internal face "<< faceI << " between cells "
-                        << own[faceI] << " and " << nei[faceI]
-                        << ".  Face area magnitude = " << a << endl;
-                }
-                else
-                {
-                    Pout<< "Zero or negative face area detected for "
-                        << "boundary face " << faceI << " next to cell "
-                        << own[faceI] << ".  Face area magnitude = "
-                        << a << endl;
-                }
-            }
-
-            nNegAreas++;
-        }
-
-        minArea = min(minArea, a);
-
-        // Create the owner pyramid - it will have negative volume
-        scalar pyrVol =
-            pyramidPointFaceRef(f[faceI], cellCtrs[own[faceI]]).mag(newPoints);
-
-        if (pyrVol > SMALL)
-        {
-            if (debug || report)
-            {
-                Pout<< "Negative pyramid volume: " << -pyrVol
-                    << " for face " << faceI << " " << f[faceI]
-                    << "  and owner cell: " << own[faceI] << endl
-                    << "Owner cell vertex labels: "
-                    << cells()[own[faceI]].labels(f)
-                    << endl;
-            }
-
-            nPyrErrors++;
-        }
-
-        if (isInternalFace(faceI))
-        {
-            // Create the neighbour pyramid - it will have positive volume
-            scalar pyrVol =
-                pyramidPointFaceRef
-                (
-                    f[faceI],
-                    cellCtrs[nei[faceI]]
-                ).mag(newPoints);
-
-            if (pyrVol < -SMALL)
-            {
-                if (debug || report)
-                {
-                    Pout<< "Negative pyramid volume: " << pyrVol
-                        << " for face " << faceI << " " << f[faceI]
-                        << "  and neighbour cell: " << nei[faceI] << nl
-                        << "Neighbour cell vertex labels: "
-                        << cells()[nei[faceI]].labels(f)
-                        << endl;
-                }
-
-                nPyrErrors++;
-            }
-
-            const vector d = cellCtrs[nei[faceI]] - cellCtrs[own[faceI]];
-            const vector& s = fAreas[faceI];
-            scalar dDotS = (d & s)/(mag(d)*mag(s) + VSMALL);
-
-            // Only write full message the first time
-            if (dDotS < SMALL && nDotProductErrors == 0)
-            {
-                // Non-orthogonality greater than 90 deg
-                WarningIn
-                (
-                    "primitiveMesh::checkMeshMotion"
-                    "(const pointField& newPoints, const bool report) const"
-                )   << "Severe non-orthogonality in mesh motion for face "
-                    << faceI
-                    << " between cells " << own[faceI] << " and " << nei[faceI]
-                    << ": Angle = " << radToDeg(::acos(dDotS))
-                    << " deg." << endl;
-
-                nDotProductErrors++;
-            }
-        }
-    }
-
-    if (nNegAreas > 0)
-    {
-        error = true;
-
-        WarningIn
-        (
-            "primitiveMesh::checkMeshMotion"
-            "(const pointField& newPoints, const bool report) const"
-        )   << "Zero or negative face area in mesh motion in " << nNegAreas
-            << " faces.  Min area: " << minArea << endl;
-    }
-    else
-    {
-        if (debug || report)
-        {
-            Pout<< "Min area = " << minArea
-                << ".  Face areas OK." << endl;
-        }
-    }
-
-    if (nPyrErrors > 0)
-    {
-        Pout<< "Detected " << nPyrErrors
-            << " negative pyramid volume in mesh motion" << endl;
-
-        error = true;
-    }
-    else
-    {
-        if (debug || report)
-        {
-            Pout<< "Pyramid volumes OK." << endl;
-        }
-    }
-
-    if (nDotProductErrors > 0)
-    {
-        Pout<< "Detected " << nDotProductErrors
-            << " in non-orthogonality in mesh motion." << endl;
-
-        error = true;
-    }
-    else
-    {
-        if (debug || report)
-        {
-            Pout<< "Non-orthogonality check OK." << endl;
-        }
-    }
-
-    if (!error && (debug || report))
-    {
-        Pout<< "Mesh motion check OK." << endl;
-    }
-
-    return error;
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C
new file mode 100644
index 0000000000..e94c8edbf7
--- /dev/null
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C
@@ -0,0 +1,638 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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 "primitiveMeshTools.H"
+#include "syncTools.H"
+#include "pyramidPointFaceRef.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::primitiveMeshTools::faceSkewness
+(
+    const primitiveMesh& mesh,
+    const pointField& p,
+    const vectorField& fCtrs,
+    const vectorField& fAreas,
+
+    const label faceI,
+    const point& ownCc,
+    const point& neiCc
+)
+{
+    vector Cpf = fCtrs[faceI] - ownCc;
+    vector d = neiCc - ownCc;
+
+    // Skewness vector
+    vector sv =
+        Cpf
+      - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + SMALL))*d;
+    vector svHat = sv/(mag(sv) + VSMALL);
+
+    // Normalisation distance calculated as the approximate distance
+    // from the face centre to the edge of the face in the direction
+    // of the skewness
+    scalar fd = 0.2*mag(d) + VSMALL;
+    const face& f = mesh.faces()[faceI];
+    forAll(f, pi)
+    {
+        fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI])));
+    }
+
+    // Normalised skewness
+    return mag(sv)/fd;
+}
+
+
+Foam::scalar Foam::primitiveMeshTools::faceOrthogonality
+(
+    const point& ownCc,
+    const point& neiCc,
+    const vector& s
+)
+{
+    vector d = neiCc - ownCc;
+
+    return (d & s)/(mag(d)*mag(s) + VSMALL);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceOrthogonality
+(
+    const primitiveMesh& mesh,
+    const vectorField& areas,
+    const vectorField& cc
+)
+{
+    const labelList& own = mesh.faceOwner();
+    const labelList& nei = mesh.faceNeighbour();
+
+    tmp<scalarField> tortho(new scalarField(mesh.nInternalFaces()));
+    scalarField& ortho = tortho();
+
+    // Internal faces
+    forAll(nei, faceI)
+    {
+        ortho[faceI] = faceOrthogonality
+        (
+            cc[own[faceI]],
+            cc[nei[faceI]],
+            areas[faceI]
+        );
+    }
+
+    return tortho;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceSkewness
+(
+    const primitiveMesh& mesh,
+    const pointField& p,
+    const vectorField& fCtrs,
+    const vectorField& fAreas,
+    const vectorField& cellCtrs
+)
+{
+    const labelList& own = mesh.faceOwner();
+    const labelList& nei = mesh.faceNeighbour();
+    const faceList& fcs = mesh.faces();
+
+    tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
+    scalarField& skew = tskew();
+
+    forAll(nei, faceI)
+    {
+        skew[faceI] = faceSkewness
+        (
+            mesh,
+            p,
+            fCtrs,
+            fAreas,
+
+            faceI,
+            cellCtrs[own[faceI]],
+            cellCtrs[nei[faceI]]
+        );
+    }
+
+
+    // Boundary faces: consider them to have only skewness error.
+    // (i.e. treat as if mirror cell on other side)
+
+    for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
+    {
+        vector Cpf = fCtrs[faceI] - cellCtrs[own[faceI]];
+
+        vector normal = fAreas[faceI];
+        normal /= mag(normal) + VSMALL;
+        vector d = normal*(normal & Cpf);
+
+
+        // Skewness vector
+        vector sv =
+            Cpf
+          - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d;
+        vector svHat = sv/(mag(sv) + VSMALL);
+
+        // Normalisation distance calculated as the approximate distance
+        // from the face centre to the edge of the face in the direction
+        // of the skewness
+        scalar fd = 0.4*mag(d) + VSMALL;
+        const face& f = fcs[faceI];
+        forAll(f, pi)
+        {
+            fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI])));
+        }
+
+        // Normalised skewness
+        skew[faceI] = mag(sv)/fd;
+    }
+
+    return tskew;
+}
+
+void Foam::primitiveMeshTools::facePyramidVolume
+(
+    const primitiveMesh& mesh,
+    const pointField& points,
+    const vectorField& ctrs,
+
+    scalarField& ownPyrVol,
+    scalarField& neiPyrVol
+)
+{
+    const labelList& own = mesh.faceOwner();
+    const labelList& nei = mesh.faceNeighbour();
+    const faceList& f = mesh.faces();
+
+    ownPyrVol.setSize(mesh.nFaces());
+    neiPyrVol.setSize(mesh.nInternalFaces());
+
+    forAll(f, faceI)
+    {
+        // Create the owner pyramid
+        ownPyrVol[faceI] = -pyramidPointFaceRef
+        (
+            f[faceI],
+            ctrs[own[faceI]]
+        ).mag(points);
+
+        if (mesh.isInternalFace(faceI))
+        {
+            // Create the neighbour pyramid - it will have positive volume
+            neiPyrVol[faceI] = pyramidPointFaceRef
+            (
+                f[faceI],
+                ctrs[nei[faceI]]
+            ).mag(points);
+        }
+    }
+}
+
+
+void Foam::primitiveMeshTools::cellClosedness
+(
+    const primitiveMesh& mesh,
+    const Vector<label>& meshD,
+    const vectorField& areas,
+    const scalarField& vols,
+
+    scalarField& openness,
+    scalarField& aratio
+)
+{
+    const labelList& own = mesh.faceOwner();
+    const labelList& nei = mesh.faceNeighbour();
+
+    // Loop through cell faces and sum up the face area vectors for each cell.
+    // This should be zero in all vector components
+
+    vectorField sumClosed(mesh.nCells(), vector::zero);
+    vectorField sumMagClosed(mesh.nCells(), vector::zero);
+
+    forAll(own, faceI)
+    {
+        // Add to owner
+        sumClosed[own[faceI]] += areas[faceI];
+        sumMagClosed[own[faceI]] += cmptMag(areas[faceI]);
+    }
+
+    forAll(nei, faceI)
+    {
+        // Subtract from neighbour
+        sumClosed[nei[faceI]] -= areas[faceI];
+        sumMagClosed[nei[faceI]] += cmptMag(areas[faceI]);
+    }
+
+
+    label nDims = 0;
+    for (direction dir = 0; dir < vector::nComponents; dir++)
+    {
+        if (meshD[dir] == 1)
+        {
+            nDims++;
+        }
+    }
+
+
+    // Check the sums
+    openness.setSize(mesh.nCells());
+    aratio.setSize(mesh.nCells());
+
+    forAll(sumClosed, cellI)
+    {
+        scalar maxOpenness = 0;
+
+        for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+        {
+            maxOpenness = max
+            (
+                maxOpenness,
+                mag(sumClosed[cellI][cmpt])
+               /(sumMagClosed[cellI][cmpt] + VSMALL)
+            );
+        }
+        openness[cellI] = maxOpenness;
+
+        // Calculate the aspect ration as the maximum of Cartesian component
+        // aspect ratio to the total area hydraulic area aspect ratio
+        scalar minCmpt = VGREAT;
+        scalar maxCmpt = -VGREAT;
+        for (direction dir = 0; dir < vector::nComponents; dir++)
+        {
+            if (meshD[dir] == 1)
+            {
+                minCmpt = min(minCmpt, sumMagClosed[cellI][dir]);
+                maxCmpt = max(maxCmpt, sumMagClosed[cellI][dir]);
+            }
+        }
+
+        scalar aspectRatio = maxCmpt/(minCmpt + VSMALL);
+        if (nDims == 3)
+        {
+            aspectRatio = max
+            (
+                aspectRatio,
+                1.0/6.0*cmptSum(sumMagClosed[cellI])/pow(vols[cellI], 2.0/3.0)
+            );
+        }
+
+        aratio[cellI] = aspectRatio;
+    }
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceConcavity
+(
+    const scalar maxSin,
+    const primitiveMesh& mesh,
+    const pointField& p,
+    const vectorField& faceAreas
+)
+{
+    const faceList& fcs = mesh.faces();
+
+    vectorField faceNormals(faceAreas);
+    faceNormals /= mag(faceNormals) + VSMALL;
+
+    tmp<scalarField> tfaceAngles(new scalarField(mesh.nFaces()));
+    scalarField& faceAngles = tfaceAngles();
+
+
+    forAll(fcs, faceI)
+    {
+        const face& f = fcs[faceI];
+
+        // Get edge from f[0] to f[size-1];
+        vector ePrev(p[f.first()] - p[f.last()]);
+        scalar magEPrev = mag(ePrev);
+        ePrev /= magEPrev + VSMALL;
+
+        scalar maxEdgeSin = 0.0;
+
+        forAll(f, fp0)
+        {
+            // Get vertex after fp
+            label fp1 = f.fcIndex(fp0);
+
+            // Normalized vector between two consecutive points
+            vector e10(p[f[fp1]] - p[f[fp0]]);
+            scalar magE10 = mag(e10);
+            e10 /= magE10 + VSMALL;
+
+            if (magEPrev > SMALL && magE10 > SMALL)
+            {
+                vector edgeNormal = ePrev ^ e10;
+                scalar magEdgeNormal = mag(edgeNormal);
+
+                if (magEdgeNormal < maxSin)
+                {
+                    // Edges (almost) aligned -> face is ok.
+                }
+                else
+                {
+                    // Check normal
+                    edgeNormal /= magEdgeNormal;
+
+                    if ((edgeNormal & faceNormals[faceI]) < SMALL)
+                    {
+                        maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
+                    }
+                }
+            }
+
+            ePrev = e10;
+            magEPrev = magE10;
+        }
+
+        faceAngles[faceI] = maxEdgeSin;
+    }
+
+    return tfaceAngles;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceFlatness
+(
+    const primitiveMesh& mesh,
+    const pointField& p,
+    const vectorField& fCtrs,
+    const vectorField& faceAreas
+)
+{
+    const faceList& fcs = mesh.faces();
+
+    // Areas are calculated as the sum of areas. (see
+    // primitiveMeshFaceCentresAndAreas.C)
+    scalarField magAreas(mag(faceAreas));
+
+    tmp<scalarField> tfaceFlatness(new scalarField(mesh.nFaces(), 1.0));
+    scalarField& faceFlatness = tfaceFlatness();
+
+
+    forAll(fcs, faceI)
+    {
+        const face& f = fcs[faceI];
+
+        if (f.size() > 3 && magAreas[faceI] > VSMALL)
+        {
+            const point& fc = fCtrs[faceI];
+
+            // Calculate the sum of magnitude of areas and compare to magnitude
+            // of sum of areas.
+
+            scalar sumA = 0.0;
+
+            forAll(f, fp)
+            {
+                const point& thisPoint = p[f[fp]];
+                const point& nextPoint = p[f.nextLabel(fp)];
+
+                // Triangle around fc.
+                vector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
+                sumA += mag(n);
+            }
+
+            faceFlatness[faceI] = magAreas[faceI] / (sumA+VSMALL);
+        }
+    }
+
+    return tfaceFlatness;
+}
+
+
+//Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::edgeAlignment
+//(
+//    const primitiveMesh& mesh,
+//    const Vector<label>& meshD,
+//    const pointField& p
+//)
+//{
+//    label nDirs = 0;
+//    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+//    {
+//        if (meshD[cmpt] == 1)
+//        {
+//            nDirs++;
+//        }
+//        else if (meshD[cmpt] != 0)
+//        {
+//            FatalErrorIn
+//            (
+//                "primitiveMeshTools::edgeAlignment"
+//                "(const primitiveMesh&, const Vector<label>&"
+//                ", const pointField&)"
+//            )   << "directions should contain 0 or 1 but is now " << meshD
+//                << exit(FatalError);
+//        }
+//    }
+//
+//    tmp<scalarField> tedgeAlignment(new scalarField(mesh.nFaces(), 0.0));
+//    scalarField& edgeAlignment = tedgeAlignment();
+//
+//    if (nDirs != vector::nComponents)
+//    {
+//        const faceList& fcs = mesh.faces();
+//
+//        forAll(fcs, faceI)
+//        {
+//            const face& f = fcs[faceI];
+//
+//            forAll(f, fp)
+//            {
+//                label p0 = f[fp];
+//                label p1 = f.nextLabel(fp);
+//                if (p0 < p1)
+//                {
+//                    vector d(p[p1]-p[p0]);
+//                    scalar magD = mag(d);
+//
+//                    if (magD > ROOTVSMALL)
+//                    {
+//                        d /= magD;
+//
+//                        // Check how many empty directions are used by the
+//                        // edge.
+//                        label nEmptyDirs = 0;
+//                        label nNonEmptyDirs = 0;
+//                        for
+//                        (
+//                              direction cmpt=0;
+//                              cmpt<vector::nComponents;
+//                              cmpt++
+//                        )
+//                        {
+//                            if (mag(d[cmpt]) > 1e-6)
+//                            {
+//                                if (meshD[cmpt] == 0)
+//                                {
+//                                    nEmptyDirs++;
+//                                }
+//                                else
+//                                {
+//                                    nNonEmptyDirs++;
+//                                }
+//                            }
+//                        }
+//
+//                        if (nEmptyDirs == 0)
+//                        {
+//                            // Purely in ok directions.
+//                        }
+//                        else if (nEmptyDirs == 1)
+//                        {
+//                            // Ok if purely in empty directions.
+//                            if (nNonEmptyDirs > 0)
+//                            {
+//                                edgeAlignment[faceI] = max
+//                                (
+//                                    edgeAlignment[faceI],
+//                                    1
+//                                );
+//                            }
+//                        }
+//                        else if (nEmptyDirs > 1)
+//                        {
+//                            // Always an error
+//                            edgeAlignment[faceI] = max
+//                            (
+//                                edgeAlignment[faceI],
+//                                2
+//                            );
+//                        }
+//                    }
+//                }
+//            }
+//        }
+//    }
+//
+//    return tedgeAlignment;
+//}
+
+
+// Checks cells with 1 or less internal faces. Give numerical problems.
+Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::cellDeterminant
+(
+    const primitiveMesh& mesh,
+    const Vector<label>& meshD,
+    const vectorField& faceAreas,
+    const PackedBoolList& internalOrCoupledFace
+)
+{
+    // Determine number of dimensions and (for 2D) missing dimension
+    label nDims = 0;
+    label twoD = -1;
+    for (direction dir = 0; dir < vector::nComponents; dir++)
+    {
+        if (meshD[dir] == 1)
+        {
+            nDims++;
+        }
+        else
+        {
+            twoD = dir;
+        }
+    }
+
+
+
+    tmp<scalarField> tcellDeterminant(new scalarField(mesh.nCells()));
+    scalarField& cellDeterminant = tcellDeterminant();
+
+    const cellList& c = mesh.cells();
+
+    if (nDims == 1)
+    {
+        cellDeterminant = 1.0;
+    }
+    else
+    {
+        forAll (c, cellI)
+        {
+            const labelList& curFaces = c[cellI];
+
+            // Calculate local normalization factor
+            scalar avgArea = 0;
+
+            label nInternalFaces = 0;
+
+            forAll(curFaces, i)
+            {
+                if (internalOrCoupledFace[curFaces[i]])
+                {
+                    avgArea += mag(faceAreas[curFaces[i]]);
+
+                    nInternalFaces++;
+                }
+            }
+
+            if (nInternalFaces == 0)
+            {
+                cellDeterminant[cellI] = 0;
+            }
+            else
+            {
+                avgArea /= nInternalFaces;
+
+                symmTensor areaTensor(symmTensor::zero);
+
+                forAll(curFaces, i)
+                {
+                    if (internalOrCoupledFace[curFaces[i]])
+                    {
+                        areaTensor += sqr(faceAreas[curFaces[i]]/avgArea);
+                    }
+                }
+
+                if (nDims == 2)
+                {
+                    // Add the missing eigenvector (such that it does not
+                    // affect the determinant)
+                    if (twoD == 0)
+                    {
+                        areaTensor.xx() = 1;
+                    }
+                    else if (twoD == 1)
+                    {
+                        areaTensor.yy() = 1;
+                    }
+                    else
+                    {
+                        areaTensor.zz() = 1;
+                    }
+                }
+
+                cellDeterminant[cellI] = mag(det(areaTensor));
+            }
+        }
+    }
+
+    return tcellDeterminant;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H
new file mode 100644
index 0000000000..5e3b2c07fd
--- /dev/null
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H
@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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/>.
+
+Namespace
+    Foam::primitiveMeshTools
+
+Description
+    Collection of static functions operating on primitiveMesh (mainly checks).
+
+SourceFiles
+    primitiveMeshTools.C
+
+\*---------------------------------------------------------------------------*/
+#ifndef primitiveMeshTools_H
+#define primitiveMeshTools_H
+
+#include "primitiveMesh.H"
+#include "PackedBoolList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Namespace primitiveMeshTools Declaration
+\*---------------------------------------------------------------------------*/
+
+class primitiveMeshTools
+{
+
+protected:
+
+   // Protected Member Functions
+
+        //- Skewness of single face
+        static scalar faceSkewness
+        (
+            const primitiveMesh& mesh,
+            const pointField& p,
+            const vectorField& fCtrs,
+            const vectorField& fAreas,
+
+            const label faceI,
+            const point& ownCc,
+            const point& neiCc
+        );
+
+        //- Orthogonality of single face
+        static scalar faceOrthogonality
+        (
+            const point& ownCc,
+            const point& neiCc,
+            const vector& s
+        );
+
+public:
+
+    //- Generate non-orthogonality field (internal faces only)
+    static tmp<scalarField> faceOrthogonality
+    (
+        const primitiveMesh& mesh,
+        const vectorField& fAreas,
+        const vectorField& cellCtrs
+    );
+
+    //- Generate face pyramid volume fields
+    static void facePyramidVolume
+    (
+        const primitiveMesh& mesh,
+        const pointField& points,
+        const vectorField& cellCtrs,
+        scalarField& ownPyrVol,
+        scalarField& neiPyrVol
+    );
+
+    //- Generate skewness field
+    static tmp<scalarField> faceSkewness
+    (
+        const primitiveMesh& mesh,
+        const pointField& points,
+        const vectorField& fCtrs,
+        const vectorField& fAreas,
+        const vectorField& cellCtrs
+    );
+
+    //- Generate cell openness and cell ascpect ratio field
+    static void cellClosedness
+    (
+        const primitiveMesh& mesh,
+        const Vector<label>& meshD,
+        const vectorField& areas,
+        const scalarField& vols,
+        scalarField& openness,
+        scalarField& aratio
+    );
+
+    //- Generate face concavity field. Returns per face the (sin of the)
+    //  most concave angle between two consecutive edges
+    static tmp<scalarField> faceConcavity
+    (
+        const scalar maxSin,
+        const primitiveMesh& mesh,
+        const pointField& p,
+        const vectorField& faceAreas
+    );
+
+    //- Generate face flatness field. Compares the individual triangles'
+    //  normals against the face average normal. Between 0 (fully warped)
+    //  and 1 (fully flat)
+    static tmp<scalarField> faceFlatness
+    (
+        const primitiveMesh& mesh,
+        const pointField& p,
+        const vectorField& fCtrs,
+        const vectorField& faceAreas
+    );
+
+    //- Generate edge alignment field. Is per face the minimum aligned edge
+    //  (does not use edge addressing)
+    static tmp<scalarField> edgeAlignment
+    (
+        const primitiveMesh& mesh,
+        const Vector<label>& directions,
+        const pointField& p
+    );
+
+    //- Generate cell determinant field
+    static tmp<scalarField> cellDeterminant
+    (
+        const primitiveMesh& mesh,
+        const Vector<label>& directions,
+        const vectorField& faceAreas,
+        const PackedBoolList& internalOrCoupledFace
+    );
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //