diff --git a/applications/utilities/mesh/advanced/combinePatchFaces/Make/options b/applications/utilities/mesh/advanced/combinePatchFaces/Make/options
index 7d4e1304f0..4f26f666de 100644
--- a/applications/utilities/mesh/advanced/combinePatchFaces/Make/options
+++ b/applications/utilities/mesh/advanced/combinePatchFaces/Make/options
@@ -5,4 +5,5 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lmeshTools \
+    -lsampling \
     -ldynamicMesh
diff --git a/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C b/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C
index 877579fa5f..8c8be624f0 100644
--- a/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C
+++ b/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.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
@@ -53,105 +53,19 @@ Description
 #include "polyMesh.H"
 #include "mapPolyMesh.H"
 #include "unitConversion.H"
+#include "motionSmoother.H"
 
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-
-// Same check as snapMesh
-void checkSnapMesh
-(
-    const Time& runTime,
-    const polyMesh& mesh,
-    labelHashSet& wrongFaces
-)
-{
-    IOdictionary snapDict
-    (
-        IOobject
-        (
-            "snapMeshDict",
-            runTime.system(),
-            mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    );
-
-    // Max nonorthogonality allowed
-    scalar maxNonOrtho(readScalar(snapDict.lookup("maxNonOrtho")));
-    // Max concaveness allowed.
-    scalar maxConcave(readScalar(snapDict.lookup("maxConcave")));
-    // Min volume allowed (factor of minimum cellVolume)
-    scalar relMinVol(readScalar(snapDict.lookup("minVol")));
-    const scalar minCellVol = min(mesh.cellVolumes());
-    const scalar minPyrVol = relMinVol*minCellVol;
-    // Min area
-    scalar minArea(readScalar(snapDict.lookup("minArea")));
-
-    if (maxNonOrtho < 180.0-SMALL)
-    {
-        Pout<< "Checking non orthogonality" << endl;
-
-        label nOldSize = wrongFaces.size();
-        mesh.setNonOrthThreshold(maxNonOrtho);
-        mesh.checkFaceOrthogonality(false, &wrongFaces);
-
-        Pout<< "Detected " << wrongFaces.size() - nOldSize
-            << " faces with non-orthogonality > " << maxNonOrtho << " degrees"
-            << endl;
-    }
-
-    if (minPyrVol > -GREAT)
-    {
-        Pout<< "Checking face pyramids" << endl;
-
-        label nOldSize = wrongFaces.size();
-        mesh.checkFacePyramids(false, minPyrVol, &wrongFaces);
-        Pout<< "Detected additional " << wrongFaces.size() - nOldSize
-            << " faces with illegal face pyramids" << endl;
-    }
-
-    if (maxConcave < 180.0-SMALL)
-    {
-        Pout<< "Checking face angles" << endl;
-
-        label nOldSize = wrongFaces.size();
-        mesh.checkFaceAngles(false, maxConcave, &wrongFaces);
-        Pout<< "Detected additional " << wrongFaces.size() - nOldSize
-            << " faces with concavity > " << maxConcave << " degrees"
-            << endl;
-    }
-
-    if (minArea > -SMALL)
-    {
-        Pout<< "Checking face areas" << endl;
-
-        label nOldSize = wrongFaces.size();
-
-        const scalarField magFaceAreas(mag(mesh.faceAreas()));
-
-        forAll(magFaceAreas, faceI)
-        {
-            if (magFaceAreas[faceI] < minArea)
-            {
-                wrongFaces.insert(faceI);
-            }
-        }
-        Pout<< "Detected additional " << wrongFaces.size() - nOldSize
-            << " faces with area < " << minArea << " m^2" << endl;
-    }
-}
-
-
 // Merge faces on the same patch (usually from exposing refinement)
 // Can undo merges if these cause problems.
 label mergePatchFaces
 (
     const scalar minCos,
     const scalar concaveSin,
-    const bool snapMeshDict,
+    const autoPtr<IOdictionary>& qualDictPtr,
     const Time& runTime,
     polyMesh& mesh
 )
@@ -212,9 +126,9 @@ label mergePatchFaces
         // Faces in error.
         labelHashSet errorFaces;
 
-        if (snapMeshDict)
+        if (qualDictPtr.valid())
         {
-            checkSnapMesh(runTime, mesh, errorFaces);
+            motionSmoother::checkMesh(false, mesh, qualDictPtr(), errorFaces);
         }
         else
         {
@@ -437,8 +351,8 @@ int main(int argc, char *argv[])
     );
     argList::addBoolOption
     (
-        "snapMesh",
-        "use system/snapMeshDict"
+        "meshQuality",
+        "read user-defined mesh quality criterions from system/meshQualityDict"
     );
 
 #   include "setRootCase.H"
@@ -455,8 +369,8 @@ int main(int argc, char *argv[])
     scalar concaveAngle = args.optionLookupOrDefault("concaveAngle", 30.0);
     scalar concaveSin = Foam::sin(degToRad(concaveAngle));
 
-    const bool snapMeshDict = args.optionFound("snapMesh");
     const bool overwrite = args.optionFound("overwrite");
+    const bool meshQuality = args.optionFound("meshQuality");
 
     Info<< "Merging all faces of a cell" << nl
         << "    - which are on the same patch" << nl
@@ -468,23 +382,47 @@ int main(int argc, char *argv[])
         << "      (sin:" << concaveSin << ')' << nl
         << endl;
 
+    autoPtr<IOdictionary> qualDict;
+    if (meshQuality)
+    {
+        Info<< "Enabling user-defined geometry checks." << nl << endl;
+
+        qualDict.reset
+        (
+            new IOdictionary
+            (
+                IOobject
+                (
+                    "meshQualityDict",
+                    mesh.time().system(),
+                    mesh,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                )
+           )
+        );
+    }
+
+
     if (!overwrite)
     {
         runTime++;
     }
 
+
+
     // Merge faces on same patch
     label nChanged = mergePatchFaces
     (
         minCos,
         concaveSin,
-        snapMeshDict,
+        qualDict,
         runTime,
         mesh
     );
 
     // Merge points on straight edges and remove unused points
-    if (snapMeshDict)
+    if (qualDict.valid())
     {
         Info<< "Merging all 'loose' points on surface edges, "
             << "regardless of the angle they make." << endl;
diff --git a/applications/utilities/mesh/generation/cvMesh/Make/options b/applications/utilities/mesh/generation/cvMesh/Make/options
index 2255a9a271..2228fc2b12 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 97e045f305..62649cf232 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 ecd7f102ba..f15ebea29a 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 c7c073ab17..35c90b1f48 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 b7fd6d3bd2..94ff17ee99 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 2372b144ec..6b63430918 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/applications/utilities/mesh/manipulation/setsToZones/Make/options b/applications/utilities/mesh/manipulation/setsToZones/Make/options
index 54c035b8f5..ec38e1cbdb 100644
--- a/applications/utilities/mesh/manipulation/setsToZones/Make/options
+++ b/applications/utilities/mesh/manipulation/setsToZones/Make/options
@@ -2,4 +2,5 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
-    -lmeshTools
+    -lmeshTools \
+    -lsampling
diff --git a/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C b/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C
index 2f56bb5693..c3d3d6297a 100644
--- a/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C
+++ b/applications/utilities/mesh/manipulation/setsToZones/setsToZones.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
@@ -48,6 +48,7 @@ Description
 #include "IFstream.H"
 #include "IOobjectList.H"
 #include "SortableList.H"
+#include "timeSelector.H"
 
 using namespace Foam;
 
@@ -57,6 +58,7 @@ using namespace Foam;
 
 int main(int argc, char *argv[])
 {
+    timeSelector::addOptions(true, false);
     argList::addNote
     (
         "add point/face/cell Zones from similar named point/face/cell Sets"
@@ -76,15 +78,7 @@ int main(int argc, char *argv[])
     const bool noFlipMap = args.optionFound("noFlipMap");
 
     // Get times list
-    instantList Times = runTime.times();
-
-    label startTime = Times.size()-1;
-    label endTime = Times.size();
-
-    // check -time and -latestTime options
-    #include "checkTimeOption.H"
-
-    runTime.setTime(Times[startTime], startTime);
+    (void)timeSelector::selectIfPresent(runTime, args);
 
     #include "createNamedPolyMesh.H"
 
diff --git a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H
index 397dc7ecc2..363f09d7f3 100644
--- a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H
+++ b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H
@@ -446,13 +446,6 @@ class vtkPV3Foam
             template<class PatchType>
             vtkPolyData* patchVTKMesh(const word& name, const PatchType&);
 
-            //- Add face zone mesh
-            vtkPolyData* faceZoneVTKMesh
-            (
-                const fvMesh&,
-                const labelList& faceLabels
-            );
-
             //- Add point zone
             vtkPolyData* pointZoneVTKMesh
             (
diff --git a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMesh.C b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMesh.C
index b229da9479..b80acb4949 100644
--- a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMesh.C
+++ b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMesh.C
@@ -450,7 +450,8 @@ void Foam::vtkPV3Foam::convertMeshFaceZones
                 << zoneName << endl;
         }
 
-        vtkPolyData* vtkmesh = faceZoneVTKMesh(mesh, zMesh[zoneId]);
+        vtkPolyData* vtkmesh = patchVTKMesh(zoneName, zMesh[zoneId]());
+
         if (vtkmesh)
         {
             AddToBlock(output, vtkmesh, range, datasetNo, zoneName);
diff --git a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshZone.C b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshZone.C
index 28a1fd4c5f..bfa3e1acc9 100644
--- a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshZone.C
+++ b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshZone.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
@@ -35,78 +35,6 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-vtkPolyData* Foam::vtkPV3Foam::faceZoneVTKMesh
-(
-    const fvMesh& mesh,
-    const labelList& faceLabels
-)
-{
-    vtkPolyData* vtkmesh = vtkPolyData::New();
-
-    if (debug)
-    {
-        Info<< "<beg> Foam::vtkPV3Foam::faceZoneVTKMesh" << endl;
-        printMemory();
-    }
-
-    // Construct primitivePatch of faces in faceZone
-
-    const faceList& meshFaces = mesh.faces();
-    faceList patchFaces(faceLabels.size());
-    forAll(faceLabels, faceI)
-    {
-        patchFaces[faceI] = meshFaces[faceLabels[faceI]];
-    }
-    primitiveFacePatch p(patchFaces, mesh.points());
-
-
-    // The balance of this routine should be identical to patchVTKMesh
-
-    // Convert OpenFOAM mesh vertices to VTK
-    const pointField& points = p.localPoints();
-
-    vtkPoints* vtkpoints = vtkPoints::New();
-    vtkpoints->Allocate(points.size());
-    forAll(points, i)
-    {
-        vtkInsertNextOpenFOAMPoint(vtkpoints, points[i]);
-    }
-
-    vtkmesh->SetPoints(vtkpoints);
-    vtkpoints->Delete();
-
-
-    // Add faces as polygons
-    const faceList& faces = p.localFaces();
-
-    vtkCellArray* vtkcells = vtkCellArray::New();
-    vtkcells->Allocate(faces.size());
-
-    forAll(faces, faceI)
-    {
-        const face& f = faces[faceI];
-        vtkIdType nodeIds[f.size()];
-
-        forAll(f, fp)
-        {
-            nodeIds[fp] = f[fp];
-        }
-        vtkcells->InsertNextCell(f.size(), nodeIds);
-    }
-
-    vtkmesh->SetPolys(vtkcells);
-    vtkcells->Delete();
-
-    if (debug)
-    {
-        Info<< "<end> Foam::vtkPV3Foam::faceZoneVTKMesh" << endl;
-        printMemory();
-    }
-
-    return vtkmesh;
-}
-
-
 vtkPolyData* Foam::vtkPV3Foam::pointZoneVTKMesh
 (
     const fvMesh& mesh,
diff --git a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamPointFields.H b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamPointFields.H
index 78582fe7e3..e74937cd2f 100644
--- a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamPointFields.H
+++ b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamPointFields.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
@@ -129,6 +129,42 @@ void Foam::vtkPV3Foam::convertPointFields
                 datasetNo
             );
         }
+
+        //
+        // Convert faceZones - if activated
+        //
+        for
+        (
+            int partId = arrayRangeFaceZones_.start();
+            partId < arrayRangeFaceZones_.end();
+            ++partId
+        )
+        {
+            const word zoneName = getPartName(partId);
+            const label datasetNo = partDataset_[partId];
+            const label zoneId = mesh.faceZones().findZoneID(zoneName);
+
+            if (!partStatus_[partId] || datasetNo < 0 || zoneId < 0)
+            {
+                continue;
+            }
+
+            // Extract the field on the zone
+            Field<Type> fld
+            (
+                ptf.internalField(),
+                mesh.faceZones()[zoneId]().meshPoints()
+            );
+
+            convertPatchPointField
+            (
+                fieldName,
+                fld,
+                output,
+                arrayRangeFaceZones_,
+                datasetNo
+            );
+        }
     }
 }
 
diff --git a/bin/tools/CleanFunctions b/bin/tools/CleanFunctions
index 758627a29a..8a6301a545 100644
--- a/bin/tools/CleanFunctions
+++ b/bin/tools/CleanFunctions
@@ -101,7 +101,7 @@ cleanCase()
     rm -rf constant/tetDualMesh > /dev/null 2>&1
 
     rm -rf VTK > /dev/null 2>&1
-    rm -f 0/cellLevel 0/pointLevel
+    rm -f 0/cellLevel 0/pointLevel 0/cellDist constant/cellDecomposition
 
     if [ -e constant/polyMesh/blockMeshDict.m4 ]
     then
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index 2fb7c7dfe2..d1061458e8 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -526,7 +526,7 @@ $(Fields)/symmTensorField/symmTensorField.C
 $(Fields)/tensorField/tensorField.C
 $(Fields)/complexFields/complexFields.C
 
-$(Fields)/labelField/labelIOField.
+$(Fields)/labelField/labelIOField.C
 $(Fields)/labelField/labelFieldIOField.C
 $(Fields)/scalarField/scalarIOField.C
 $(Fields)/scalarField/scalarFieldIOField.C
diff --git a/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C b/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C
index 43d14a5740..1ea9ae6e11 100644
--- a/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.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
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "SlicedGeometricField.H"
+#include "processorFvPatch.H"
 
 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
 
@@ -40,7 +41,8 @@ slicedBoundaryField
 (
     const Mesh& mesh,
     const Field<Type>& completeField,
-    const bool preserveCouples
+    const bool preserveCouples,
+    const bool preserveProcessorOnly
 )
 {
     tmp<FieldField<PatchField, Type> > tbf
@@ -52,7 +54,15 @@ slicedBoundaryField
 
     forAll(mesh.boundary(), patchi)
     {
-        if (preserveCouples && mesh.boundary()[patchi].coupled())
+        if
+        (
+            preserveCouples
+         && mesh.boundary()[patchi].coupled()
+         && (
+               !preserveProcessorOnly
+            || isA<processorFvPatch>(mesh.boundary()[patchi])
+            )
+        )
         {
             // For coupled patched construct the correct patch field type
             bf.set
@@ -243,7 +253,8 @@ SlicedGeometricField
     const dimensionSet& ds,
     const Field<Type>& completeIField,
     const Field<Type>& completeBField,
-    const bool preserveCouples
+    const bool preserveCouples,
+    const bool preserveProcessorOnly
 )
 :
     GeometricField<Type, PatchField, GeoMesh>
@@ -252,7 +263,13 @@ SlicedGeometricField
         mesh,
         ds,
         Field<Type>(),
-        slicedBoundaryField(mesh, completeBField, preserveCouples)
+        slicedBoundaryField
+        (
+            mesh,
+            completeBField,
+            preserveCouples,
+            preserveProcessorOnly
+        )
     )
 {
     // Set the internalField to the slice of the complete field
diff --git a/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.H b/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.H
index e04690e042..dd23712e67 100644
--- a/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.H
+++ b/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.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
@@ -83,7 +83,8 @@ private:
         (
             const Mesh& mesh,
             const Field<Type>& completeField,
-            const bool preserveCouples
+            const bool preserveCouples,
+            const bool preserveProcessorOnly = false
         );
 
         //- Slice the given field and a create a PtrList of SlicedPatchField
@@ -133,7 +134,8 @@ public:
             const dimensionSet&,
             const Field<Type>& completeIField,
             const Field<Type>& completeBField,
-            const bool preserveCouples=true
+            const bool preserveCouples=true,
+            const bool preserveProcessorOnly = false
         );
 
         //- Construct from GeometricField. Reuses full internal and
diff --git a/src/conversion/meshReader/createPolyBoundary.C b/src/conversion/meshReader/createPolyBoundary.C
index e0c898156a..47fb2008f7 100644
--- a/src/conversion/meshReader/createPolyBoundary.C
+++ b/src/conversion/meshReader/createPolyBoundary.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
@@ -303,14 +303,12 @@ void Foam::meshReader::createPolyBoundary()
 
     Info<< "Added " << nMissingFaces << " unmatched faces" << endl;
 
+    // Add missing faces to last patch ('Default_Empty' etc.)
     if (nMissingFaces > 0)
     {
         patchSizes_.last() = nMissingFaces;
     }
-    else
-    {
-        patchStarts_.setSize(patchStarts_.size() - 1);
-    }
+
 
     // reset the size of the face list
     meshFaces_.setSize(nCreatedFaces);
diff --git a/src/finiteVolume/fvMesh/fvMeshGeometry.C b/src/finiteVolume/fvMesh/fvMeshGeometry.C
index 0a06ef7364..5abdb45d11 100644
--- a/src/finiteVolume/fvMesh/fvMeshGeometry.C
+++ b/src/finiteVolume/fvMesh/fvMeshGeometry.C
@@ -133,6 +133,8 @@ void fvMesh::makeC() const
             << abort(FatalError);
     }
 
+    // Construct as slices. Only preserve processor (not e.g. cyclic)
+
     CPtr_ = new slicedVolVectorField
     (
         IOobject
@@ -148,33 +150,10 @@ void fvMesh::makeC() const
         *this,
         dimLength,
         cellCentres(),
-        faceCentres()
+        faceCentres(),
+        true,               //preserveCouples
+        true                //preserveProcOnly
     );
-
-
-    // Need to correct for cyclics transformation since absolute quantity.
-    // Ok on processor patches since hold opposite cell centre (no
-    // transformation)
-    slicedVolVectorField& C = *CPtr_;
-
-    forAll(C.boundaryField(), patchi)
-    {
-        if
-        (
-            isA<cyclicFvPatchVectorField>(C.boundaryField()[patchi])
-         || isA<cyclicAMIFvPatchVectorField>(C.boundaryField()[patchi])
-        )
-        {
-            // Note: cyclic is not slice but proper field
-            C.boundaryField()[patchi] == static_cast<const vectorField&>
-            (
-                static_cast<const List<vector>&>
-                (
-                    boundary_[patchi].patchSlice(faceCentres())
-                )
-            );
-        }
-    }
 }
 
 
diff --git a/src/fvMotionSolver/Make/options b/src/fvMotionSolver/Make/options
index 7c440dd78f..fa13513b50 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 ca8cb9e2e4..0ee4f07bb0 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/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index 7cf1e79f61..e4cea95be2 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -167,7 +167,7 @@ class autoSnapDriver
                 void correctAttraction
                 (
                     const DynamicList<point>& surfacePoints,
-                    const DynamicList<label>& surfaceCount,
+                    const DynamicList<label>& surfaceCounts,
                     const point& edgePt,
                     const vector& edgeNormal,   // normalised normal
                     const point& pt,
@@ -213,7 +213,7 @@ class autoSnapDriver
 
                     DynamicList<point>& surfacePoints,
                     DynamicList<vector>& surfaceNormals,
-                    DynamicList<label>& surfaceCount
+                    DynamicList<label>& surfaceCounts
                 ) const;
                 void binFeatureFaces
                 (
@@ -224,13 +224,13 @@ class autoSnapDriver
                     const scalarField& snapDist,
                     const label pointI,
 
-                    const List<List<point> >& pointFaceNormals,
+                    const List<List<point> >& pointFaceSurfNormals,
                     const List<List<point> >& pointFaceDisp,
                     const List<List<point> >& pointFaceCentres,
 
                     DynamicList<point>& surfacePoints,
                     DynamicList<vector>& surfaceNormals,
-                    DynamicList<label>& surfaceCount
+                    DynamicList<label>& surfaceCounts
                 ) const;
 
 
@@ -259,7 +259,7 @@ class autoSnapDriver
                     const indirectPrimitivePatch& pp,
                     const scalarField& snapDist,
 
-                    const List<List<point> >& pointFaceNormals,
+                    const List<List<point> >& pointFaceSurfNormals,
                     const List<List<point> >& pointFaceDisp,
                     const List<List<point> >& pointFaceCentres,
                     const labelListList& pointFacePatchID,
@@ -277,7 +277,7 @@ class autoSnapDriver
                     const indirectPrimitivePatch&,
                     const scalarField&,
 
-                    const List<List<point> >& pointFaceNormals,
+                    const List<List<point> >& pointFaceSurfNormals,
                     const List<List<point> >& pointFaceDisp,
                     const List<List<point> >& pointFaceCentres,
                     const labelListList& pointFacePatchID,
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index f32d9ce4bd..b9990cb078 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -690,7 +690,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
 void Foam::autoSnapDriver::correctAttraction
 (
     const DynamicList<point>& surfacePoints,
-    const DynamicList<label>& surfaceCount,
+    const DynamicList<label>& surfaceCounts,
     const point& edgePt,
     const vector& edgeNormal,       // normalised normal
     const point& pt,
@@ -702,7 +702,7 @@ void Foam::autoSnapDriver::correctAttraction
     scalar tang = ((pt-edgePt)&edgeNormal);
 
     labelList order;
-    Foam::sortedOrder(surfaceCount, order);
+    Foam::sortedOrder(surfaceCounts, order);
 
     if (order[0] < order[1])
     {
@@ -763,7 +763,7 @@ void Foam::autoSnapDriver::binFeatureFace
 
     DynamicList<point>& surfacePoints,
     DynamicList<vector>& surfaceNormals,
-    DynamicList<label>& surfaceCount
+    DynamicList<label>& surfaceCounts
 ) const
 {
     // What to do with very far attraction? For now just ignore the face
@@ -783,7 +783,7 @@ void Foam::autoSnapDriver::binFeatureFace
             )
             {
                 same = true;
-                surfaceCount[j]++;
+                surfaceCounts[j]++;
                 break;
             }
         }
@@ -796,7 +796,7 @@ void Foam::autoSnapDriver::binFeatureFace
             {
                 surfacePoints.append(pt);
                 surfaceNormals.append(faceSurfaceNormal);
-                surfaceCount.append(1);
+                surfaceCounts.append(1);
             }
             else if (surfacePoints.size() == 2)
             {
@@ -810,7 +810,7 @@ void Foam::autoSnapDriver::binFeatureFace
                     // Definitely makes a feature point
                     surfacePoints.append(pt);
                     surfaceNormals.append(faceSurfaceNormal);
-                    surfaceCount.append(1);
+                    surfaceCounts.append(1);
                 }
             }
             else if (surfacePoints.size() == 3)
@@ -834,7 +834,7 @@ void Foam::autoSnapDriver::binFeatureFace
                         // Different feature point
                         surfacePoints.append(pt);
                         surfaceNormals.append(faceSurfaceNormal);
-                        surfaceCount.append(1);
+                        surfaceCounts.append(1);
                     }
                 }
             }
@@ -860,15 +860,15 @@ void Foam::autoSnapDriver::binFeatureFaces
 
     DynamicList<point>& surfacePoints,
     DynamicList<vector>& surfaceNormals,
-    DynamicList<label>& surfaceCount
+    DynamicList<label>& surfaceCounts
 ) const
 {
-    const List<point>& pfNormals = pointFaceSurfNormals[pointI];
+    const List<point>& pfSurfNormals = pointFaceSurfNormals[pointI];
     const List<point>& pfDisp = pointFaceDisp[pointI];
     const List<point>& pfCentres = pointFaceCentres[pointI];
 
     // Collect all different directions
-    forAll(pfNormals, i)
+    forAll(pfSurfNormals, i)
     {
         binFeatureFace
         (
@@ -879,12 +879,12 @@ void Foam::autoSnapDriver::binFeatureFaces
             snapDist[pointI],
 
             pfCentres[i],
-            pfNormals[i],
+            pfSurfNormals[i],
             pfDisp[i],
 
             surfacePoints,
             surfaceNormals,
-            surfaceCount
+            surfaceCounts
         );
     }
 }
@@ -914,7 +914,7 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
     // Collect all different directions
     DynamicList<point> surfacePoints(4);
     DynamicList<vector> surfaceNormals(4);
-    DynamicList<label> surfaceCount(4);
+    DynamicList<label> surfaceCounts(4);
 
     binFeatureFaces
     (
@@ -931,7 +931,7 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
 
         surfacePoints,
         surfaceNormals,
-        surfaceCount
+        surfaceCounts
     );
 
     const point& pt = pp.localPoints()[pointI];
@@ -966,11 +966,13 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
         vector d = r.refPoint()-pt;
         d -= (d&n)*n;
 
+        //- This does not help much but distorts a perfectly aligned mesh
+        //  so disabled for now.
         //// Correct for attraction to non-dominant face
         //correctAttraction
         //(
         //    surfacePoints,
-        //    surfaceCount,
+        //    surfaceCounts,
         //    r.refPoint(),
         //    n,                  // normalised normal
         //    pt,
@@ -1813,8 +1815,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     }
 
     // Reverse: from pp point to nearest feature
-    vectorField allPatchAttraction(pp.nPoints(), vector::zero);
-    List<pointConstraint> allPatchConstraints(pp.nPoints());
+    vectorField rawPatchAttraction(pp.nPoints(), vector::zero);
+    List<pointConstraint> rawPatchConstraints(pp.nPoints());
 
     determineFeatures
     (
@@ -1837,15 +1839,15 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
         edgeAttractors,
         edgeConstraints,
         // pp point to nearest feature
-        allPatchAttraction,
-        allPatchConstraints
+        rawPatchAttraction,
+        rawPatchConstraints
     );
 
 
 
     // Baffle handling
     // ~~~~~~~~~~~~~~~
-    // Override pointAttractor, edgeAttractor, allPatchAttration etc. to
+    // Override pointAttractor, edgeAttractor, rawPatchAttration etc. to
     // implement 'baffle' handling.
     // Baffle: the mesh pp point originates from a loose standing
     // baffle.
@@ -1983,8 +1985,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                     featI,
                     edgeAttractors,
                     edgeConstraints,
-                    allPatchAttraction,
-                    allPatchConstraints
+                    rawPatchAttraction,
+                    rawPatchConstraints
                 );
 
                 if (!nearInfo.hit())
@@ -2033,8 +2035,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                             vector::zero;
 
                         // Store for later use
-                        allPatchAttraction[pointI] = featPt-pt;
-                        allPatchConstraints[pointI] =
+                        rawPatchAttraction[pointI] = featPt-pt;
+                        rawPatchConstraints[pointI] =
                             pointConstraints[featI][featPointI];
 
                         if (oldPointI != -1)
@@ -2053,8 +2055,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                                 edgeFeatI,
                                 edgeAttractors,
                                 edgeConstraints,
-                                allPatchAttraction,
-                                allPatchConstraints
+                                rawPatchAttraction,
+                                rawPatchConstraints
                             );
                         }
                     }
@@ -2084,8 +2086,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                         featI,
                         edgeAttractors,
                         edgeConstraints,
-                        allPatchAttraction,
-                        allPatchConstraints
+                        rawPatchAttraction,
+                        rawPatchConstraints
                     );
                 }
             }
@@ -2296,11 +2298,11 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                 if
                 (
                     patchConstraints[pointI].first() <= 1
-                 && allPatchConstraints[pointI].first() > 1
+                 && rawPatchConstraints[pointI].first() > 1
                 )
                 {
-                    patchAttraction[pointI] = allPatchAttraction[pointI];
-                    patchConstraints[pointI] = allPatchConstraints[pointI];
+                    patchAttraction[pointI] = rawPatchAttraction[pointI];
+                    patchConstraints[pointI] = rawPatchConstraints[pointI];
 
                     if (multiPatchStr.valid())
                     {
@@ -2434,7 +2436,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     // Snap edges to feature edges
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Walk existing edges and snap remaining ones (that are marked as
-    // feature edges in allPatchConstraints)
+    // feature edges in rawPatchConstraints)
 
     stringFeatureEdges
     (
@@ -2444,8 +2446,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
         pp,
         snapDist,
 
-        allPatchAttraction,
-        allPatchConstraints,
+        rawPatchAttraction,
+        rawPatchConstraints,
 
         patchAttraction,
         patchConstraints
diff --git a/src/meshTools/Make/options b/src/meshTools/Make/options
index 45765759c1..ac717c7516 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 c42133c29f..653d708e49 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 cc319679d4..107e343ef3 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.
diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.C b/src/meshTools/searchableSurface/triSurfaceMesh.C
index 86fcd5889b..2508cdb919 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.C
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.C
@@ -625,6 +625,12 @@ Foam::triSurfaceMesh::edgeTree() const
 }
 
 
+Foam::scalar Foam::triSurfaceMesh::tolerance() const
+{
+    return tolerance_;
+}
+
+
 const Foam::wordList& Foam::triSurfaceMesh::regions() const
 {
     if (regions_.empty())
diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.H b/src/meshTools/searchableSurface/triSurfaceMesh.H
index a1d0378487..aa462e2554 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.H
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.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
@@ -182,6 +182,7 @@ public:
         //- Demand driven contruction of octree for boundary edges
         const indexedOctree<treeDataEdge>& edgeTree() const;
 
+        scalar tolerance() const;
 
         // searchableSurface implementation
 
diff --git a/src/surfMesh/surfaceFormats/obj/OBJstream.C b/src/surfMesh/surfaceFormats/obj/OBJstream.C
index ca24b6b6aa..fcd446eaab 100644
--- a/src/surfMesh/surfaceFormats/obj/OBJstream.C
+++ b/src/surfMesh/surfaceFormats/obj/OBJstream.C
@@ -238,6 +238,38 @@ Foam::Ostream& Foam::OBJstream::write(const linePointRef& ln)
 }
 
 
+Foam::Ostream& Foam::OBJstream::write
+(
+    const triPointRef& f,
+    const bool lines
+)
+{
+    label start = nVertices_;
+    write(f.a());
+    write(f.b());
+    write(f.c());
+    if (lines)
+    {
+        write('l');
+        for (int i = 0; i < 3; i++)
+        {
+            write(' ') << start+1+i;
+        }
+        write(' ') << start+1 << '\n';
+    }
+    else
+    {
+        write('f');
+        for (int i = 0; i < 3; i++)
+        {
+            write(' ') << start+1+i;
+        }
+        write('\n');
+    }
+    return *this;
+}
+
+
 Foam::Ostream& Foam::OBJstream::write
 (
     const face& f,
diff --git a/src/surfMesh/surfaceFormats/obj/OBJstream.H b/src/surfMesh/surfaceFormats/obj/OBJstream.H
index 9e62cbdbeb..ec7e93d24e 100644
--- a/src/surfMesh/surfaceFormats/obj/OBJstream.H
+++ b/src/surfMesh/surfaceFormats/obj/OBJstream.H
@@ -39,13 +39,14 @@ SourceFiles
 #include "point.H"
 #include "edge.H"
 #include "face.H"
+#include "triPointRef.H"
+#include "linePointRef.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-
 /*---------------------------------------------------------------------------*\
                           Class OBJstream Declaration
 \*---------------------------------------------------------------------------*/
@@ -134,6 +135,9 @@ public:
             //- Write line
             Ostream& write(const linePointRef&);
 
+            //- Write triangle as points with lines or filled polygon
+            Ostream& write(const triPointRef&, const bool lines = true);
+
             //- Write face as points with lines or filled polygon
             Ostream& write
             (
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/constant/polyMesh/blockMeshDict b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/constant/polyMesh/blockMeshDict
index 707a53ea74..57de714b2f 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/constant/polyMesh/blockMeshDict
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/constant/polyMesh/blockMeshDict
@@ -31,7 +31,7 @@ vertices
 
 blocks
 (
-    hex (0 1 2 3 4 5 6 7) (20 20 20) simpleGrading (1 1 1)
+    hex (0 1 2 3 4 5 6 7) (20 20 10) simpleGrading (1 1 1)
 );
 
 edges
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict
index 8297dae5af..788a479390 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict
@@ -42,38 +42,6 @@ geometry
         max (1 1 1);
     }
 
-    fridgeFreezer
-    {
-        type searchableSurfaceCollection;
-
-        mergeSubRegions true;
-
-        freezer
-        {
-            surface box1;
-            scale (1 1 1);
-            transform
-            {
-                type    cartesian;
-                origin  (0 0 0);
-                e1      (1 0 0);
-                e3      (0 0 1);
-            }
-        }
-        fridge
-        {
-            surface box1;
-            scale (1 1 1.1);
-            transform
-            {
-                type    cartesian;
-                origin  (0 0 1);
-                e1      (1 0 0);
-                e3      (0 0 1);
-            }
-        }
-    }
-
     twoFridgeFreezers
     {
         type searchableSurfaceCollection;
@@ -150,7 +118,7 @@ castellatedMeshControls
     (
         {
             file "fridgeA.eMesh";
-            level 3;
+            levels ((0.01 3));
         }
     );
 
diff --git a/tutorials/incompressible/pimpleFoam/channel395/Allrun b/tutorials/incompressible/pimpleFoam/channel395/Allrun
index 61c1b54530..7f75af44c2 100755
--- a/tutorials/incompressible/pimpleFoam/channel395/Allrun
+++ b/tutorials/incompressible/pimpleFoam/channel395/Allrun
@@ -8,7 +8,15 @@ cd ${0%/*} || exit 1    # run from this directory
 application=`getApplication`
 
 runApplication blockMesh
-runApplication $application
+
+#- Run serial
+#runApplication $application
+
+#- Run parallel
+runApplication decomposePar -cellDist
+runParallel $application 5
+runApplication reconstructPar
+
 runApplication postChannel
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/incompressible/pimpleFoam/channel395/constant/polyMesh/boundary b/tutorials/incompressible/pimpleFoam/channel395/constant/polyMesh/boundary
index 655c5c0b54..12dc93fdc2 100644
--- a/tutorials/incompressible/pimpleFoam/channel395/constant/polyMesh/boundary
+++ b/tutorials/incompressible/pimpleFoam/channel395/constant/polyMesh/boundary
@@ -32,6 +32,7 @@ FoamFile
     sides1_half0
     {
         type            cyclic;
+        inGroups        1(cyclic);
         nFaces          1000;
         startFace       177700;
         matchTolerance  0.0001;
@@ -40,6 +41,7 @@ FoamFile
     sides1_half1
     {
         type            cyclic;
+        inGroups        1(cyclic);
         nFaces          1000;
         startFace       178700;
         matchTolerance  0.0001;
@@ -48,6 +50,7 @@ FoamFile
     sides2_half0
     {
         type            cyclic;
+        inGroups        1(cyclic);
         nFaces          1000;
         startFace       179700;
         matchTolerance  0.0001;
@@ -56,6 +59,7 @@ FoamFile
     sides2_half1
     {
         type            cyclic;
+        inGroups        1(cyclic);
         nFaces          1000;
         startFace       180700;
         matchTolerance  0.0001;
@@ -64,6 +68,7 @@ FoamFile
     inout1_half0
     {
         type            cyclic;
+        inGroups        1(cyclic);
         nFaces          750;
         startFace       181700;
         matchTolerance  0.0001;
@@ -72,6 +77,7 @@ FoamFile
     inout1_half1
     {
         type            cyclic;
+        inGroups        1(cyclic);
         nFaces          750;
         startFace       182450;
         matchTolerance  0.0001;
@@ -80,6 +86,7 @@ FoamFile
     inout2_half0
     {
         type            cyclic;
+        inGroups        1(cyclic);
         nFaces          750;
         startFace       183200;
         matchTolerance  0.0001;
@@ -88,6 +95,7 @@ FoamFile
     inout2_half1
     {
         type            cyclic;
+        inGroups        1(cyclic);
         nFaces          750;
         startFace       183950;
         matchTolerance  0.0001;
diff --git a/tutorials/incompressible/pimpleFoam/channel395/system/decomposeParDict b/tutorials/incompressible/pimpleFoam/channel395/system/decomposeParDict
new file mode 100644
index 0000000000..9f0365f2de
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/channel395/system/decomposeParDict
@@ -0,0 +1,135 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "mesh decomposition control dictionary";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  5;
+
+//- Keep owner and neighbour on same processor for faces in zones:
+// preserveFaceZones (heater solid1 solid3);
+
+//- Keep owner and neighbour on same processor for faces in patches:
+//  (makes sense only for cyclic patches)
+//preservePatches (cyclic_half0 cyclic_half1);
+
+//- Keep all of faceSet on a single processor. This puts all cells
+//  connected with a point, edge or face on the same processor.
+//  (just having face connected cells might not guarantee a balanced
+//  decomposition)
+// The processor can be -1 (the decompositionMethod chooses the processor
+// for a good load balance) or explicitly provided (upsets balance).
+//singleProcessorFaceSets ((f0 -1));
+
+
+//- Use the volScalarField named here as a weight for each cell in the
+//  decomposition.  For example, use a particle population field to decompose
+//  for a balanced number of particles in a lagrangian simulation.
+// weightField dsmcRhoNMean;
+
+method          scotch;
+//method          hierarchical;
+// method          simple;
+// method          metis;
+// method          manual;
+// method          multiLevel;
+// method          structured;  // does 2D decomposition of structured mesh
+
+multiLevelCoeffs
+{
+    // Decomposition methods to apply in turn. This is like hierarchical but
+    // fully general - every method can be used at every level.
+
+    level0
+    {
+        numberOfSubdomains  64;
+        //method simple;
+        //simpleCoeffs
+        //{
+        //    n           (2 1 1);
+        //    delta       0.001;
+        //}
+        method scotch;
+    }
+    level1
+    {
+        numberOfSubdomains  4;
+        method scotch;
+    }
+}
+
+// Desired output
+
+simpleCoeffs
+{
+    n           (2 1 1);
+    delta       0.001;
+}
+
+hierarchicalCoeffs
+{
+    n           (1 2 1);
+    delta       0.001;
+    order       xyz;
+}
+
+metisCoeffs
+{
+ /*
+    processorWeights
+    (
+        1
+        1
+        1
+        1
+    );
+  */
+}
+
+scotchCoeffs
+{
+    //processorWeights
+    //(
+    //    1
+    //    1
+    //    1
+    //    1
+    //);
+    //writeGraph  true;
+    //strategy "b";
+}
+
+manualCoeffs
+{
+    dataFile    "decompositionData";
+}
+
+structuredCoeffs
+{
+    // Patches to do 2D decomposition on. Structured mesh only; cells have
+    // to be in 'columns' on top of patches.
+    patches     (bottomPatch);
+}
+
+//// Is the case distributed? Note: command-line argument -roots takes
+//// precedence
+//distributed     yes;
+//// Per slave (so nProcs-1 entries) the directory above the case.
+//roots
+//(
+//    "/tmp"
+//    "/tmp"
+//);
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/motorBike/Allrun b/tutorials/incompressible/simpleFoam/motorBike/Allrun
index 714dba48dc..d3466dcfab 100755
--- a/tutorials/incompressible/simpleFoam/motorBike/Allrun
+++ b/tutorials/incompressible/simpleFoam/motorBike/Allrun
@@ -6,6 +6,7 @@ cd ${0%/*} || exit 1    # run from this directory
 
 # copy motorbike surface from resources folder
 cp $FOAM_TUTORIALS/resources/geometry/motorBike.obj.gz constant/triSurface/
+runApplication surfaceFeatureExtract
 
 runApplication blockMesh
 
@@ -19,6 +20,7 @@ runParallel snappyHexMesh 6 -overwrite
 ls -d processor* | xargs -i rm -rf ./{}/0 $1
 ls -d processor* | xargs -i cp -r 0.org ./{}/0 $1
 
+runParallel patchSummary 6
 runParallel potentialFoam 6 -noFunctionObjects -writep
 runParallel `getApplication` 6
 
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict b/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
index 68819f01c4..517ef75b9c 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
@@ -90,10 +90,10 @@ castellatedMeshControls
     // This is a featureEdgeMesh, read from constant/triSurface for now.
     features
     (
-        //{
-        //    file "someLine.eMesh";
-        //    level 2;
-        //}
+        {
+            file "motorBike.eMesh";
+            level 0;
+        }
     );
 
 
@@ -189,10 +189,21 @@ snapControls
     //  before upon reaching a correct mesh.
     nRelaxIter 5;
 
-    //- Highly experimental and wip: number of feature edge snapping
-    //  iterations. Leave out altogether to disable.
-    //  Do not use here since mesh resolution too low and baffles present
-    //nFeatureSnapIter 10;
+    // Feature snapping
+
+        //- Number of feature edge snapping iterations.
+        //  Leave out altogether to disable.
+        nFeatureSnapIter 10;
+
+        //- Detect (geometric only) features by sampling the surface
+        //  (default=false).
+        implicitFeatureSnap false;
+
+        //- Use castellatedMeshControls::features (default = true)
+        explicitFeatureSnap true;
+
+        //- Detect points on multiple surfaces (only for explicitFeatureSnap)
+        multiRegionFeatureSnap true;
 }
 
 
@@ -239,7 +250,11 @@ addLayersControls
 
     //- When not to extrude surface. 0 is flat surface, 90 is when two faces
     //  make straight angle.
-    featureAngle 30;
+    featureAngle 60;
+
+    //- At non-patched sides allow mesh to slip if extrusion direction makes
+    //  angle larger than slipFeatureAngle.
+    slipFeatureAngle 30;
 
     //- Maximum number of snapping relaxation iterations. Should stop
     //  before upon reaching a correct mesh.
diff --git a/tutorials/lagrangian/reactingParcelFoam/filter/Allrun b/tutorials/lagrangian/reactingParcelFoam/filter/Allrun
index b7428e064f..b82b7aad36 100755
--- a/tutorials/lagrangian/reactingParcelFoam/filter/Allrun
+++ b/tutorials/lagrangian/reactingParcelFoam/filter/Allrun
@@ -13,9 +13,6 @@ runApplication blockMesh
 #setSet -batch system/sets.setSet > log.setSet1 2>&1
 runApplication topoSet
 
-# convert sets to zones
-setsToZones -noFlipMap > log.setsToZones 2>&1
-
 # create the first cyclic - lhs of porous zone
 # Note that we don't know what value to give these patches-out-of-nothing so
 # - use binary writing to avoid 'nan'
diff --git a/tutorials/lagrangian/reactingParcelFoam/filter/system/topoSetDict b/tutorials/lagrangian/reactingParcelFoam/filter/system/topoSetDict
index fd7ee0e747..e77e55346f 100644
--- a/tutorials/lagrangian/reactingParcelFoam/filter/system/topoSetDict
+++ b/tutorials/lagrangian/reactingParcelFoam/filter/system/topoSetDict
@@ -28,6 +28,17 @@ actions
             box (1.5 -10 -10) (2 10 10);
         }
     }
+    {
+        name    filter;
+        type    cellZoneSet;
+        action  new;
+        source  setToCellZone;
+        sourceInfo
+        {
+            set filter;
+        }
+    }
+
 
     {
         name    leftFluid;
@@ -39,6 +50,16 @@ actions
             box (-10 -10 -10) (1.5 10 10);
         }
     }
+    {
+        name    leftFluid;
+        type    cellZoneSet;
+        action  new;
+        source  setToCellZone;
+        sourceInfo
+        {
+            set leftFluid;
+        }
+    }
     {
         name    rightFluid;
         type    cellSet;
@@ -49,6 +70,16 @@ actions
             box (2 -1 -1) (10 10 10);
         }
     }
+    {
+        name    rightFluid;
+        type    cellZoneSet;
+        action  new;
+        source  setToCellZone;
+        sourceInfo
+        {
+            set rightFluid;
+        }
+    }
 
 
     // cycLeft
@@ -74,6 +105,17 @@ actions
             option  all;
         }
     }
+    // Create faceZone from cycLeft
+    {
+        name    cycLeft;
+        type    faceZoneSet;
+        action  new;
+        source  setToFaceZone;
+        sourceInfo
+        {
+            faceSet cycLeft;    // name of faceSet
+        }
+    }
 
     // cycRight
     {
@@ -98,6 +140,18 @@ actions
             option  all;
         }
     }
+    // Create faceZone from cycRight
+    {
+        name    cycRight;
+        type    faceZoneSet;
+        action  new;
+        source  setToFaceZone;
+        sourceInfo
+        {
+            faceSet cycRight;   // name of faceSet
+        }
+    }
+
 );
 
 // ************************************************************************* //