diff --git a/applications/utilities/mesh/generation/cvMesh/Make/options b/applications/utilities/mesh/generation/cvMesh/Make/options
index 2255a9a27164813355db2e5b730673a3441dc4c9..2228fc2b121966cd0063aa3a5afeadd2f4a31202 100644
--- a/applications/utilities/mesh/generation/cvMesh/Make/options
+++ b/applications/utilities/mesh/generation/cvMesh/Make/options
@@ -17,6 +17,7 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
     -I$(LIB_SRC)/edgeMesh/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude
 
@@ -29,5 +30,6 @@ EXE_LIBS = \
     -ldecompositionMethods \
     -L$(FOAM_LIBBIN)/dummy -lptscotchDecomp \
     -ledgeMesh \
+    -lsampling \
     -ltriSurface \
     -ldynamicMesh
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options
index 97e045f30502c9dc417b7d110e9492f03d29b289..62649cf23282eb467b02bbc51b116c928b76e4c8 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options
@@ -17,6 +17,7 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
     -I$(LIB_SRC)/edgeMesh/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I../vectorTools
diff --git a/applications/utilities/mesh/generation/cvMesh/cvMesh.C b/applications/utilities/mesh/generation/cvMesh/cvMesh.C
index ecd7f102baa7ab22e1161473c7defc665ef58a3d..f15ebea29abb5ff77b1e8566a28d39a5d6050942 100644
--- a/applications/utilities/mesh/generation/cvMesh/cvMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/cvMesh.C
@@ -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
@@ -36,6 +36,7 @@ Description
 
 #include "argList.H"
 #include "conformalVoronoiMesh.H"
+#include "vtkSetWriter.H"
 
 using namespace Foam;
 
@@ -48,6 +49,11 @@ int main(int argc, char *argv[])
         "noFilter",
         "Do not filter the mesh"
     );
+    Foam::argList::addBoolOption
+    (
+        "checkGeometry",
+        "check all surface geometry for quality"
+    );
 
     #include "setRootCase.H"
     #include "createTime.H"
@@ -55,6 +61,7 @@ int main(int argc, char *argv[])
     runTime.functionObjects().off();
 
     const bool noFilter = !args.optionFound("noFilter");
+    const bool checkGeometry = args.optionFound("checkGeometry");
 
     Info<< "Mesh filtering is " << (noFilter ? "on" : "off") << endl;
 
@@ -74,6 +81,29 @@ int main(int argc, char *argv[])
 
     conformalVoronoiMesh mesh(runTime, cvMeshDict);
 
+
+    if (checkGeometry)
+    {
+        const searchableSurfaces& allGeometry = mesh.allGeometry();
+
+        // Write some stats
+        allGeometry.writeStats(List<wordList>(0), Info);
+        // Check topology
+        allGeometry.checkTopology(true);
+        // Check geometry
+        allGeometry.checkGeometry
+        (
+            100.0,      // max size ratio
+            1e-9,       // intersection tolerance
+            autoPtr<writer<scalar> >(new vtkSetWriter<scalar>()),
+            0.01,       // min triangle quality
+            true
+        );
+
+        return 0;
+    }
+
+
     while (runTime.loop())
     {
         Info<< nl << "Time = " << runTime.timeName() << endl;
diff --git a/applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/Make/options b/applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/Make/options
index c7c073ab1764bc3672e63a10d8aa3ec580b88156..35c90b1f48a814892bce172cb3a1a4d31152b163 100644
--- a/applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/Make/options
+++ b/applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/Make/options
@@ -9,6 +9,7 @@ EXE_INC = \
     -I$(FASTDUALOCTREE_SRC_PATH) \
     -I../conformalVoronoiMesh/lnInclude \
     -I$(LIB_SRC)/edgeMesh/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
@@ -21,6 +22,7 @@ EXE_LIBS = \
     -lconformalVoronoiMesh \
     -ldecompositionMethods -L$(FOAM_LIBBIN)/dummy -lscotchDecomp \
     -ledgeMesh \
+    -lsampling \
     -ltriSurface \
     -lmeshTools \
     -ldynamicMesh
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/Make/options b/applications/utilities/mesh/generation/snappyHexMesh/Make/options
index b7fd6d3bd227e7ed7d2bf9f5b6c79121d34a1446..94ff17ee9979e9c63d98f3513d7d24a564277430 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/Make/options
+++ b/applications/utilities/mesh/generation/snappyHexMesh/Make/options
@@ -3,6 +3,7 @@ EXE_INC = \
     -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
     -I$(LIB_SRC)/mesh/autoMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/edgeMesh/lnInclude \
@@ -13,5 +14,6 @@ EXE_LIBS = \
     -ldecompositionMethods \
     -L$(FOAM_LIBBIN)/dummy -lptscotchDecomp \
     -lmeshTools \
+    -lsampling \
     -ldynamicMesh \
     -lautoMesh
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index 2372b144ec9af2d0a7e5920b3e22f9ff8281e407..6b6343091872e0e1dc15997225e917f1a7146516 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -46,7 +46,7 @@ Description
 #include "refinementParameters.H"
 #include "snapParameters.H"
 #include "layerParameters.H"
-
+#include "vtkSetWriter.H"
 
 using namespace Foam;
 
@@ -122,6 +122,12 @@ void writeMesh
 int main(int argc, char *argv[])
 {
 #   include "addOverwriteOption.H"
+    Foam::argList::addBoolOption
+    (
+        "checkGeometry",
+        "check all surface geometry for quality"
+    );
+
 #   include "setRootCase.H"
 #   include "createTime.H"
     runTime.functionObjects().off();
@@ -131,6 +137,7 @@ int main(int argc, char *argv[])
         << runTime.cpuTimeIncrement() << " s" << endl;
 
     const bool overwrite = args.optionFound("overwrite");
+    const bool checkGeometry = args.optionFound("checkGeometry");
 
     // Check patches and faceZones are synchronised
     mesh.boundaryMesh().checkParallelSync(true);
@@ -244,6 +251,56 @@ int main(int argc, char *argv[])
         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
 
 
+    // Checking only?
+
+    if (checkGeometry)
+    {
+        // Extract patchInfo
+        List<wordList> patchTypes(allGeometry.size());
+
+        const PtrList<dictionary>& patchInfo = surfaces.patchInfo();
+        const labelList& surfaceGeometry = surfaces.surfaces();
+        forAll(surfaceGeometry, surfI)
+        {
+            label geomI = surfaceGeometry[surfI];
+            const wordList& regNames = allGeometry.regionNames()[geomI];
+
+            patchTypes[geomI].setSize(regNames.size());
+            forAll(regNames, regionI)
+            {
+                label globalRegionI = surfaces.globalRegion(surfI, regionI);
+
+                if (patchInfo.set(globalRegionI))
+                {
+                    patchTypes[geomI][regionI] =
+                        word(patchInfo[globalRegionI].lookup("type"));
+                }
+                else
+                {
+                    patchTypes[geomI][regionI] = wallPolyPatch::typeName;
+                }
+            }
+        }
+
+        // Write some stats
+        allGeometry.writeStats(patchTypes, Info);
+        // Check topology
+        allGeometry.checkTopology(true);
+        // Check geometry
+        allGeometry.checkGeometry
+        (
+            100.0,      // max size ratio
+            1e-9,       // intersection tolerance
+            autoPtr<writer<scalar> >(new vtkSetWriter<scalar>()),
+            0.01,       // min triangle quality
+            true
+        );
+
+        return 0;
+    }
+
+
+
     // Read refinement shells
     // ~~~~~~~~~~~~~~~~~~~~~~
 
diff --git a/src/fvMotionSolver/Make/options b/src/fvMotionSolver/Make/options
index 7c440dd78fcfee40befe27a4b976b0a5d0d1943e..fa13513b50fbf34e09087c06eb56bc413d357a5e 100644
--- a/src/fvMotionSolver/Make/options
+++ b/src/fvMotionSolver/Make/options
@@ -3,6 +3,7 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/postProcessing/functionObjects/forces/lnInclude \
 
 LIB_LIBS = \
diff --git a/src/mesh/autoMesh/Make/options b/src/mesh/autoMesh/Make/options
index ca8cb9e2e44d51cc79d846a27a48e1731738e775..0ee4f07bb039ccdf6b15f6182b77b60c47d89ec0 100644
--- a/src/mesh/autoMesh/Make/options
+++ b/src/mesh/autoMesh/Make/options
@@ -4,6 +4,7 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/edgeMesh/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude
diff --git a/src/meshTools/Make/options b/src/meshTools/Make/options
index 45765759c16acadd0d907685d44e79783df5cc54..ac717c75165175856a5b5d27b056357166ed0956 100644
--- a/src/meshTools/Make/options
+++ b/src/meshTools/Make/options
@@ -1,5 +1,6 @@
 EXE_INC = \
-    -I$(LIB_SRC)/triSurface/lnInclude
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
 
 LIB_LIBS = \
     -ltriSurface
diff --git a/src/meshTools/searchableSurface/searchableSurfaces.C b/src/meshTools/searchableSurface/searchableSurfaces.C
index c42133c29f490d586fa1a3b9bad59b97d93052ef..653d708e49df678f72e8b0546ad5df7ccf3174c1 100644
--- a/src/meshTools/searchableSurface/searchableSurfaces.C
+++ b/src/meshTools/searchableSurface/searchableSurfaces.C
@@ -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
@@ -27,12 +27,42 @@ License
 #include "searchableSurfacesQueries.H"
 #include "ListOps.H"
 #include "Time.H"
+//#include "vtkSetWriter.H"
+#include "DynamicField.H"
+//#include "OBJstream.H"
+#include "PatchTools.H"
+#include "triSurfaceMesh.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 defineTypeNameAndDebug(Foam::searchableSurfaces, 0);
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+//- Is edge connected to triangle
+bool Foam::searchableSurfaces::connected
+(
+    const triSurface& s,
+    const label edgeI,
+    const pointIndexHit& hit
+)
+{
+    const triFace& localFace = s.localFaces()[hit.index()];
+    const edge& e = s.edges()[edgeI];
+
+    forAll(localFace, i)
+    {
+        if (e.otherVertex(localFace[i]) != -1)
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct with length.
@@ -378,6 +408,468 @@ Foam::pointIndexHit Foam::searchableSurfaces::facesIntersection
     );
 }
 
+
+bool Foam::searchableSurfaces::checkClosed(const bool report) const
+{
+    if (report)
+    {
+        Info<< "Checking for closedness." << endl;
+    }
+
+    bool hasError = false;
+
+    forAll(*this, surfI)
+    {
+        if (!operator[](surfI).hasVolumeType())
+        {
+            hasError = true;
+
+            if (report)
+            {
+                Info<< "    " << names()[surfI]
+                    << " : not closed" << endl;
+            }
+
+            if (isA<triSurface>(operator[](surfI)))
+            {
+                const triSurface& s = dynamic_cast<const triSurface&>
+                (
+                    operator[](surfI)
+                );
+                const labelListList& edgeFaces = s.edgeFaces();
+
+                label nSingleEdges = 0;
+                forAll(edgeFaces, edgeI)
+                {
+                    if (edgeFaces[edgeI].size() == 1)
+                    {
+                        nSingleEdges++;
+                    }
+                }
+
+                label nMultEdges = 0;
+                forAll(edgeFaces, edgeI)
+                {
+                    if (edgeFaces[edgeI].size() > 2)
+                    {
+                        nMultEdges++;
+                    }
+                }
+
+                if (report && (nSingleEdges != 0 || nMultEdges != 0))
+                {
+                    Info<< "        connected to one face : "
+                        << nSingleEdges << nl
+                        << "        connected to >2 faces : "
+                        << nMultEdges << endl;
+                }
+            }
+        }
+    }
+
+    if (report)
+    {
+        Info<< endl;
+    }
+
+    return returnReduce(hasError, orOp<bool>());
+}
+
+
+bool Foam::searchableSurfaces::checkNormalOrientation(const bool report) const
+{
+    if (report)
+    {
+        Info<< "Checking for normal orientation." << endl;
+    }
+
+    bool hasError = false;
+
+    forAll(*this, surfI)
+    {
+        if (isA<triSurface>(operator[](surfI)))
+        {
+            const triSurface& s = dynamic_cast<const triSurface&>
+            (
+                operator[](surfI)
+            );
+
+            labelHashSet borderEdge(s.size()/1000);
+            PatchTools::checkOrientation(s, false, &borderEdge);
+            // Colour all faces into zones using borderEdge
+            labelList normalZone;
+            label nZones = PatchTools::markZones(s, borderEdge, normalZone);
+            if (nZones > 1)
+            {
+                hasError = true;
+
+                if (report)
+                {
+                    Info<< "    " << names()[surfI]
+                        << " : has multiple orientation zones ("
+                        << nZones << ")" << endl;
+                }
+            }
+        }
+    }
+
+    if (report)
+    {
+        Info<< endl;
+    }
+
+    return returnReduce(hasError, orOp<bool>());
+}
+
+
+bool Foam::searchableSurfaces::checkSizes
+(
+    const scalar maxRatio,
+    const bool report
+) const
+{
+    if (report)
+    {
+        Info<< "Checking for size." << endl;
+    }
+
+    bool hasError = false;
+
+    forAll(*this, i)
+    {
+        const boundBox& bb = operator[](i).bounds();
+
+        for (label j = i+1; j < size(); j++)
+        {
+            scalar ratio = bb.mag()/operator[](j).bounds().mag();
+
+            if (ratio > maxRatio || ratio < 1.0/maxRatio)
+            {
+                hasError = true;
+
+                if (report)
+                {
+                    Info<< "    " << names()[i]
+                        << " bounds differ from " << names()[j]
+                        << " by more than a factor 100:" << nl
+                        << "        bounding box : " << bb << nl
+                        << "        bounding box : " << operator[](j).bounds()
+                        << endl;
+                }
+                break;
+            }
+        }
+    }
+
+    if (report)
+    {
+        Info<< endl;
+    }
+
+    return returnReduce(hasError, orOp<bool>());
+}
+
+
+bool Foam::searchableSurfaces::checkIntersection
+(
+    const scalar tolerance,
+    const autoPtr<writer<scalar> >& setWriter,
+    const bool report
+) const
+{
+    if (report)
+    {
+        Info<< "Checking for intersection." << endl;
+    }
+
+    //cpuTime timer;
+
+    bool hasError = false;
+
+    forAll(*this, i)
+    {
+        if (isA<triSurfaceMesh>(operator[](i)))
+        {
+            const triSurfaceMesh& s0 = dynamic_cast<const triSurfaceMesh&>
+            (
+                operator[](i)
+            );
+            const edgeList& edges0 = s0.edges();
+            const pointField& localPoints0 = s0.localPoints();
+
+            // Collect all the edge vectors
+            pointField start(edges0.size());
+            pointField end(edges0.size());
+            forAll(edges0, edgeI)
+            {
+                const edge& e = edges0[edgeI];
+                start[edgeI] = localPoints0[e[0]];
+                end[edgeI] = localPoints0[e[1]];
+            }
+
+            // Test all other surfaces for intersection
+            forAll(*this, j)
+            {
+                List<pointIndexHit> hits;
+
+                if (i == j)
+                {
+                    // Slightly shorten the edges to prevent finding lots of
+                    // intersections. The fast triangle intersection routine
+                    // has problems with rays passing through a point of the
+                    // triangle.
+                    vectorField d
+                    (
+                        max(tolerance, 10*s0.tolerance())
+                       *(end-start)
+                    );
+                    start += d;
+                    end -= d;
+                }
+
+                operator[](j).findLineAny(start, end, hits);
+
+                label nHits = 0;
+                DynamicField<point> intersections(edges0.size()/100);
+                DynamicField<scalar> intersectionEdge(intersections.capacity());
+
+                forAll(hits, edgeI)
+                {
+                    if
+                    (
+                        hits[edgeI].hit()
+                     && (i != j || !connected(s0, edgeI, hits[edgeI]))
+                    )
+                    {
+                        intersections.append(hits[edgeI].hitPoint());
+                        intersectionEdge.append(1.0*edgeI);
+                        nHits++;
+                    }
+                }
+
+                if (nHits > 0)
+                {
+                    if (report)
+                    {
+                        Info<< "    " << names()[i]
+                            << " intersects " << names()[j]
+                            << " at " << nHits
+                            << " locations."
+                            << endl;
+
+                        //vtkSetWriter<scalar> setWriter;
+                        if (setWriter.valid())
+                        {
+                            scalarField dist(mag(intersections));
+                            coordSet track
+                            (
+                                names()[i] + '_' + names()[j],
+                                "xyz",
+                                intersections.xfer(),
+                                dist
+                            );
+                            wordList valueSetNames(1, "edgeIndex");
+                            List<const scalarField*> valueSets
+                            (
+                                1,
+                                &intersectionEdge
+                            );
+
+                            fileName fName
+                            (
+                                setWriter().getFileName(track, valueSetNames)
+                            );
+                            Info<< "    Writing intersection locations to "
+                                << fName << endl;
+                            OFstream os
+                            (
+                                s0.searchableSurface::time().path()
+                               /fName
+                            );
+                            setWriter().write
+                            (
+                                track,
+                                valueSetNames,
+                                valueSets,
+                                os
+                            );
+                        }
+                    }
+
+                    hasError = true;
+                    break;
+                }
+            }
+        }
+    }
+
+    if (report)
+    {
+        Info<< endl;
+    }
+
+    return returnReduce(hasError, orOp<bool>());
+}
+
+
+bool Foam::searchableSurfaces::checkQuality
+(
+    const scalar minQuality,
+    const bool report
+) const
+{
+    if (report)
+    {
+        Info<< "Checking for triangle quality." << endl;
+    }
+
+    bool hasError = false;
+
+    forAll(*this, surfI)
+    {
+        if (isA<triSurface>(operator[](surfI)))
+        {
+            const triSurface& s = dynamic_cast<const triSurface&>
+            (
+                operator[](surfI)
+            );
+
+            label nBadTris = 0;
+            forAll(s, faceI)
+            {
+                const labelledTri& f = s[faceI];
+
+                scalar q = triPointRef
+                (
+                    s.points()[f[0]],
+                    s.points()[f[1]],
+                    s.points()[f[2]]
+                ).quality();
+
+                if (q < minQuality)
+                {
+                    nBadTris++;
+                }
+            }
+
+            if (nBadTris > 0)
+            {
+                hasError = true;
+
+                if (report)
+                {
+                    Info<< "    " << names()[surfI]
+                        << " : has " << nBadTris << " bad quality triangles "
+                        << " (quality < " << minQuality << ")" << endl;
+                }
+            }
+        }
+    }
+
+    if (report)
+    {
+        Info<< endl;
+    }
+
+    return returnReduce(hasError, orOp<bool>());
+
+}
+
+
+// Check all surfaces. Return number of failures.
+Foam::label Foam::searchableSurfaces::checkTopology
+(
+    const bool report
+) const
+{
+    label noFailedChecks = 0;
+
+    if (checkClosed(report))
+    {
+        noFailedChecks++;
+    }
+
+    if (checkNormalOrientation(report))
+    {
+        noFailedChecks++;
+    }
+    return noFailedChecks;
+}
+
+
+Foam::label Foam::searchableSurfaces::checkGeometry
+(
+    const scalar maxRatio,
+    const scalar tol,
+    const autoPtr<writer<scalar> >& setWriter,
+    const scalar minQuality,
+    const bool report
+) const
+{
+    label noFailedChecks = 0;
+
+    if (maxRatio > 0 && checkSizes(maxRatio, report))
+    {
+        noFailedChecks++;
+    }
+
+    if (checkIntersection(tol, setWriter, report))
+    {
+        noFailedChecks++;
+    }
+
+    if (checkQuality(minQuality, report))
+    {
+        noFailedChecks++;
+    }
+
+    return noFailedChecks;
+}
+
+
+void Foam::searchableSurfaces::writeStats
+(
+    const List<wordList>& patchTypes,
+    Ostream& os
+) const
+{
+    Info<< "Statistics." << endl;
+    forAll(*this, surfI)
+    {
+        Info<< "    " << names()[surfI] << ':' << endl;
+
+        const searchableSurface& s = operator[](surfI);
+
+        Info<< "        type      : " << s.type() << nl
+            << "        size      : " << s.globalSize() << nl;
+        if (isA<triSurfaceMesh>(s))
+        {
+            const triSurfaceMesh& ts = dynamic_cast<const triSurfaceMesh&>(s);
+            Info<< "        edges     : " << ts.nEdges() << nl
+                << "        points    : " << ts.points().size() << nl;
+        }
+        Info<< "        bounds    : " << s.bounds() << nl
+            << "        closed    : " << Switch(s.hasVolumeType()) << endl;
+
+        if (patchTypes.size() && patchTypes[surfI].size() >= 1)
+        {
+            wordList unique(HashSet<word>(patchTypes[surfI]).sortedToc());
+            Info<< "        patches   : ";
+            forAll(unique, i)
+            {
+                Info<< unique[i];
+                if (i < unique.size()-1)
+                {
+                    Info<< ',';
+                }
+            }
+            Info<< endl;
+        }
+    }
+    Info<< endl;
+}
+
+
 // * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
 
 const Foam::searchableSurface& Foam::searchableSurfaces::operator[]
diff --git a/src/meshTools/searchableSurface/searchableSurfaces.H b/src/meshTools/searchableSurface/searchableSurfaces.H
index cc319679d4fbc85fc912c8f8e7f1bbd58653c705..107e343ef3d262f077a105dfce0441f757b435a5 100644
--- a/src/meshTools/searchableSurface/searchableSurfaces.H
+++ b/src/meshTools/searchableSurface/searchableSurfaces.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
@@ -37,6 +37,7 @@ SourceFiles
 
 #include "searchableSurface.H"
 #include "labelPair.H"
+#include "writer.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +45,7 @@ namespace Foam
 {
 
 // Forward declaration of classes
+class triSurface;
 
 /*---------------------------------------------------------------------------*\
                            Class searchableSurfaces Declaration
@@ -70,6 +72,15 @@ class searchableSurfaces
 
     // Private Member Functions
 
+        //- Is edge on face
+        static bool connected
+        (
+            const triSurface& s,
+            const label edgeI,
+            const pointIndexHit& hit
+        );
+
+
         //- Disallow default bitwise copy construct
         searchableSurfaces(const searchableSurfaces&);
 
@@ -189,6 +200,49 @@ public:
             ) const;
 
 
+        // Checking
+
+            //- Are all surfaces closed and manifold
+            bool checkClosed(const bool report) const;
+
+            //- Are all (triangulated) surfaces consistent normal orientation
+            bool checkNormalOrientation(const bool report) const;
+
+            //- Are all bounding boxes of similar size
+            bool checkSizes(const scalar maxRatio, const bool report) const;
+
+            //- Do surfaces self-intersect or intersect others
+            bool checkIntersection
+            (
+                const scalar tol,
+                const autoPtr<writer<scalar> >&,
+                const bool report
+            ) const;
+
+            //- Check triangle quality
+            bool checkQuality
+            (
+                const scalar minQuality,
+                const bool report
+            ) const;
+
+            //- All topological checks. Return number of failed checks
+            label checkTopology(const bool report) const;
+
+            //- All geometric checks. Return number of failed checks
+            label checkGeometry
+            (
+                const scalar maxRatio,
+                const scalar tolerance,
+                const autoPtr<writer<scalar> >& setWriter,
+                const scalar minQuality,
+                const bool report
+            ) const;
+
+            //- Write some stats
+            void writeStats(const List<wordList>&, Ostream&) const;
+
+
     // Member Operators
 
         //- Return const and non-const reference to searchableSurface by index.