diff --git a/.gitignore b/.gitignore
index 96c3093b71..150288aae2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -55,7 +55,7 @@ doc/[Dd]oxygen/man
 # ignore .timeStamp in the main directory
 /.timeStamp
 
-# ignore .ebrowse in the main directory
-/.ebrowse
+# ignore .tags in the main directory
+/.tags
 
 # end-of-file
diff --git a/applications/solvers/incompressible/simpleFoam/Make/options b/applications/solvers/incompressible/simpleFoam/Make/options
index 61091811bf..1223bdd06f 100644
--- a/applications/solvers/incompressible/simpleFoam/Make/options
+++ b/applications/solvers/incompressible/simpleFoam/Make/options
@@ -8,5 +8,4 @@ EXE_INC = \
 EXE_LIBS = \
     -lincompressibleRASModels \
     -lincompressibleTransportModels \
-    -lfiniteVolume \
-    -lmeshTools
+    -lfiniteVolume
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index 43d550f406..58f324160e 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -316,6 +316,11 @@ meshQualityControls
     //must be >0 for Fluent compatibility
     minTriangleTwist -1;
 
+    //- if >0 : preserve single cells with all points on the surface if the
+    //  resulting volume after snapping is larger than minVolFraction times old
+    //  volume. If <0 : delete always.
+    minVolFraction 0.1;
+
 
     // Advanced
 
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
index 52f23658a8..0adcd30927 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
@@ -23,16 +23,13 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
     scalar minDistSqr = magSqr(1e-6 * globalBb.span());
 
     // Non-empty directions
-    const Vector<label> validDirs = (mesh.directions() + Vector<label>::one)/2;
+    const Vector<label> validDirs = (mesh.geometricD() + Vector<label>::one)/2;
+    Info<< "    Mesh (non-empty, non-wedge) directions " << validDirs << endl;
 
-    Info<< "    Mesh (non-empty) directions " << validDirs << endl;
+    const Vector<label> solDirs = (mesh.solutionD() + Vector<label>::one)/2;
+    Info<< "    Mesh (non-empty) directions " << solDirs << endl;
 
-    scalar nGeomDims = mesh.nGeometricD();
-
-    Info<< "    Mesh (non-empty, non-wedge) dimensions "
-        << nGeomDims << endl;
-
-    if (nGeomDims < 3)
+    if (mesh.nGeometricD() < 3)
     {
         pointSet nonAlignedPoints(mesh, "nonAlignedEdges", mesh.nPoints()/100);
 
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/PV3FoamReader_SM.xml b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/PV3FoamReader_SM.xml
index 644b794a44..ef25ac0378 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/PV3FoamReader_SM.xml
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/PV3FoamReader_SM.xml
@@ -30,16 +30,16 @@
 
     <!-- Global settings -->
 
-    <!-- Extrapolate Walls check-box -->
+    <!-- Extrapolate Patches check-box -->
     <IntVectorProperty
-      name="ExtrapolateWalls"
-      command="SetExtrapolateWalls"
+      name="ExtrapolatePatches"
+      command="SetExtrapolatePatches"
       number_of_elements="1"
       default_values="0"
       animateable="0">
       <BooleanDomain name="bool"/>
       <Documentation>
-        Extrapolate internalField to wall and empty patches
+        Extrapolate internalField to non-constraint patches
       </Documentation>
     </IntVectorProperty>
 
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.cxx b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.cxx
index 8fc518d745..c16a43deae 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.cxx
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.cxx
@@ -64,7 +64,7 @@ vtkPV3FoamReader::vtkPV3FoamReader()
 
     CacheMesh = 1;
 
-    ExtrapolateWalls = 0;
+    ExtrapolatePatches = 0;
     IncludeSets = 0;
     IncludeZones = 0;
     ShowPatchNames = 0;
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.h b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.h
index fde87527de..bc21cd8ce8 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.h
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.h
@@ -65,9 +65,9 @@ public:
     vtkGetMacro(CacheMesh, int);
 
     // Description:
-    // FOAM extrapolate internal values onto the walls
-    vtkSetMacro(ExtrapolateWalls, int);
-    vtkGetMacro(ExtrapolateWalls, int);
+    // FOAM extrapolate internal values onto the patches
+    vtkSetMacro(ExtrapolatePatches, int);
+    vtkGetMacro(ExtrapolatePatches, int);
 
     // FOAM read sets control
     vtkSetMacro(IncludeSets, int);
@@ -183,7 +183,7 @@ private:
     int TimeStepRange[2];
     int CacheMesh;
 
-    int ExtrapolateWalls;
+    int ExtrapolatePatches;
     int IncludeSets;
     int IncludeZones;
     int ShowPatchNames;
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
index 7a8fd89e9e..0ed3d53185 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
@@ -659,29 +659,55 @@ void Foam::vtkPV3Foam::addPatchNames(vtkRenderer* renderer)
         }
     }
 
+    // Count number of zones we're actually going to display. This is truncated
+    // to a max per patch
+
+    const label MAXPATCHZONES = 20;
+
+    label displayZoneI = 0;
+
+    forAll(pbMesh, patchI)
+    {
+        displayZoneI += min(MAXPATCHZONES, nZones[patchI]);
+    }
+
+
     zoneCentre.shrink();
 
     if (debug)
     {
         Info<< "patch zone centres = " << zoneCentre << nl
+            << "displayed zone centres = " << displayZoneI << nl
             << "zones per patch = " << nZones << endl;
     }
 
     // Set the size of the patch labels to max number of zones
-    patchTextActorsPtrs_.setSize(zoneCentre.size());
+    patchTextActorsPtrs_.setSize(displayZoneI);
 
     if (debug)
     {
         Info<< "constructing patch labels" << endl;
     }
 
+    // Actor index
+    displayZoneI = 0;
+
+    // Index in zone centres
     label globalZoneI = 0;
+
     forAll(pbMesh, patchI)
     {
         const polyPatch& pp = pbMesh[patchI];
 
         // Only selected patches will have a non-zero number of zones
-        for (label i=0; i<nZones[patchI]; i++)
+        label nDisplayZones = min(MAXPATCHZONES, nZones[patchI]);
+        label increment = 1;
+        if (nZones[patchI] >= MAXPATCHZONES)
+        {
+            increment = nZones[patchI]/MAXPATCHZONES;
+        }
+
+        for (label i = 0; i < nDisplayZones; i++)
         {
             if (debug)
             {
@@ -719,14 +745,15 @@ void Foam::vtkPV3Foam::addPatchNames(vtkRenderer* renderer)
 
             // Maintain a list of text labels added so that they can be
             // removed later
-            patchTextActorsPtrs_[globalZoneI] = txt;
+            patchTextActorsPtrs_[displayZoneI] = txt;
 
-            globalZoneI++;
+            globalZoneI += increment;
+            displayZoneI++;
         }
     }
 
     // Resize the patch names list to the actual number of patch names added
-    patchTextActorsPtrs_.setSize(globalZoneI);
+    patchTextActorsPtrs_.setSize(displayZoneI);
 
     if (debug)
     {
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H
index f217d1a889..6c8d32f470 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H
@@ -132,8 +132,8 @@ void Foam::vtkPV3Foam::convertVolFields
                 isType<emptyFvPatchField<Type> >(ptf)
              ||
                 (
-                    typeid(patches[patchId]) == typeid(wallPolyPatch)
-                 && reader_->GetExtrapolateWalls()
+                    reader_->GetExtrapolatePatches()
+                && !polyPatch::constraintType(patches[patchId].type())
                 )
             )
             {
diff --git a/bin/foamEbrowse b/bin/foamEbrowse
index e303a77264..66fafc1e2b 100755
--- a/bin/foamEbrowse
+++ b/bin/foamEbrowse
@@ -27,26 +27,27 @@
 #     foamEbrowse
 #
 # Description
-#     Build the Ebrowse database for all the .C and .H files
+#     Build the Ebrowse database for all the .H and .C files
 #
 #------------------------------------------------------------------------------
-headersFile=${TMPDIR:-/tmp}/headersFile.$$
 sourcesFile=${TMPDIR:-/tmp}/sourcesFile.$$
 
 if [ $# -ne 0 ]; then
    echo "Usage : ${0##*/}"
    echo ""
-   echo "Build the Ebrowse dadbase for all the .C and .H files"
+   echo "Build the Ebrowse dadbase for all the .H and .C files"
    echo ""
    exit 1
 fi
 
 # Clean up on termination and on Ctrl-C
-trap 'rm -f $headersFile $sourcesFile 2>/dev/null; exit 0' EXIT TERM INT
+trap 'rm -f $sourcesFile 2>/dev/null; exit 0' EXIT TERM INT
 
 cd $WM_PROJECT_DIR
-find -H . -name "*.H" | fgrep -v lnInclude > $headersFile
-find -H . -name "*.C" | fgrep -v lnInclude > $sourcesFile
-ebrowse --files=$headersFile --files=$sourcesFile --output-file=.ebrowse
+mkdir .tags 2>/dev/null
+cd .tags
+
+find -H .. \( -name "*.[HC]" -not -name "lnInclude" -not -name "Doxygen" \) -print > $sourcesFile
+ebrowse --files=$sourcesFile --output-file=ebrowse
 
 #------------------------------------------------------------------------------
diff --git a/bin/foamTags b/bin/foamTags
new file mode 100755
index 0000000000..1cfd186f52
--- /dev/null
+++ b/bin/foamTags
@@ -0,0 +1,56 @@
+#!/bin/sh
+#------------------------------------------------------------------------------
+# =========                 |
+# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+#  \\    /   O peration     |
+#   \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+#    \\/     M anipulation  |
+#-------------------------------------------------------------------------------
+# License
+#     This file is part of OpenFOAM.
+#
+#     OpenFOAM is free software; you can redistribute it and/or modify it
+#     under the terms of the GNU General Public License as published by the
+#     Free Software Foundation; either version 2 of the License, or (at your
+#     option) any later version.
+#
+#     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+#     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+#     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+#     for more details.
+#
+#     You should have received a copy of the GNU General Public License
+#     along with OpenFOAM; if not, write to the Free Software Foundation,
+#     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+#
+# Script
+#     foamTags
+#
+# Description
+#     Build the tags files for all the .C and .H files
+#
+#------------------------------------------------------------------------------
+
+if [ $# -ne 0 ]; then
+   echo "Usage : ${0##*/}"
+   echo ""
+   echo "Build the tags files for all the .C and .H files"
+   echo ""
+   exit 1
+fi
+
+cd $WM_PROJECT_DIR
+mkdir .tags 2>/dev/null
+
+find -H . \( -name "*.[HC]" -not -name "lnInclude" -not -name "Doxygen" \) | \
+    etags --declarations -l c++ -o .tags/etags -
+find -H . \( -name "*.[HC]" -not -name "lnInclude" -not -name "Doxygen" \) | \
+    etags -l c++ -o .tags/etagsDef -
+find -H . \( -name "*.H" -not -name "lnInclude" -not -name "Doxygen" \) | \
+    etags --declarations -l c++ -o .tags/etagsDec -
+
+gtags -i --gtagsconf bin/tools/gtagsrc .tags
+
+foamEbrowse
+
+#------------------------------------------------------------------------------
diff --git a/bin/tools/gtagsrc b/bin/tools/gtagsrc
new file mode 100644
index 0000000000..ce66e0c084
--- /dev/null
+++ b/bin/tools/gtagsrc
@@ -0,0 +1,63 @@
+#------------------------------------------------------------------------------
+# =========                 |
+# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+#  \\    /   O peration     |
+#   \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+#    \\/     M anipulation  |
+#-------------------------------------------------------------------------------
+# License
+#     This file is part of OpenFOAM.
+#
+#     OpenFOAM is free software; you can redistribute it and/or modify it
+#     under the terms of the GNU General Public License as published by the
+#     Free Software Foundation; either version 2 of the License, or (at your
+#     option) any later version.
+#
+#     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+#     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+#     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+#     for more details.
+#
+#     You should have received a copy of the GNU General Public License
+#     along with OpenFOAM; if not, write to the Free Software Foundation,
+#     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+#
+# Canfiguration file
+#     gtagsrc
+#
+# Description
+#     Configuration file for gtags(1).
+#
+#------------------------------------------------------------------------------
+default:\
+	:tc=gtags:tc=htags:
+#------------------------------------------------------------------------------
+# Configuration for gtags(1)
+# See gtags(1).
+#------------------------------------------------------------------------------
+common:\
+	:skip=GPATH,GTAGS,GRTAGS,GSYMS,HTML/,HTML.pub/,html/,tags,TAGS,ID,.ebrowse,.etags,.etagsDef,.etagsDec,y.tab.c,y.tab.h,.notfunction,cscope.out,cscope.po.out,cscope.in.out,.gdbinit,SCCS/,RCS/,CVS/,CVSROOT/,{arch}/,.svn/,.git/,.cvsrc,.cvsignore,.gitignore,.cvspass,.cvswrappers,.deps/,autom4te.cache/,.snprj/:\
+	:langmap=c\:.c.h,yacc\:.y,asm\:.s.S,java\:.java,cpp\:.c++.cc.cpp.cxx.hxx.hpp.C.H,php\:.php.php3.phtml:
+gtags:\
+	:tc=common:\
+	:GTAGS=gtags-parser %s:\
+	:GRTAGS=gtags-parser -r %s:\
+	:GSYMS=gtags-parser -s %s:\
+	:skip=lnInclude/,tutorials/,wmake/,doc/,lib/,etc/:
+#
+#------------------------------------------------------------------------------
+# Configuration for htags(1)
+#------------------------------------------------------------------------------
+htags:\
+	:body_begin=<body text='#191970' bgcolor='#f5f5dc' vlink='gray'>:body_end=</body>:\
+	:table_begin=<table>:table_end=</table>:\
+	:title_begin=<h1><font color='#cc0000'>:title_end=</font></h1>:\
+	:comment_begin=<i><font color='green'>:comment_end=</font></i>:\
+	:sharp_begin=<font color='darkred'>:sharp_end=</font>:\
+	:brace_begin=<font color='red'>:brace_end=</font>:\
+	:warned_line_begin=<span style='background-color\:yellow'>:warned_line_end=</span>:\
+	:reserved_begin=<b>:reserved_end=</b>:script_alias=/cgi-bin/:\
+	:ncol#4:tabs#8:normal_suffix=html:gzipped_suffix=ghtml:\
+	:definition_header=no:
+
+#------------------------------------------------------------------------------
\ No newline at end of file
diff --git a/src/Allwmake b/src/Allwmake
index f22e0729db..c3ce648ba0 100755
--- a/src/Allwmake
+++ b/src/Allwmake
@@ -5,8 +5,8 @@ set -x
 # update OpenFOAM version strings if required
 wmakePrintBuild -check || /bin/rm -f OpenFOAM/Make/$WM_OPTIONS/global.? 2>/dev/null
 
-wmakeLnInclude -f OpenFOAM
-wmakeLnInclude -f OSspecific/$WM_OS
+wmakeLnInclude -f OpenFOAM -sf
+wmakeLnInclude -f OSspecific/$WM_OS -sf
 Pstream/Allwmake
 
 wmake libo  OSspecific/$WM_OS
diff --git a/src/OpenFOAM/db/Time/findInstance.C b/src/OpenFOAM/db/Time/findInstance.C
index f21a524c80..a19d3bda63 100644
--- a/src/OpenFOAM/db/Time/findInstance.C
+++ b/src/OpenFOAM/db/Time/findInstance.C
@@ -58,8 +58,9 @@ Foam::word Foam::Time::findInstance
     {
         if (debug)
         {
-            Info<< "Time::findInstance(const fileName&, const word&) : "
-                << "found \"" << name
+            Info<< "Time::findInstance"
+                "(const fileName&, const word&, const IOobject::readOption)"
+                << " : found \"" << name
                 << "\" in " << timeName()/dir
                 << endl;
         }
@@ -98,8 +99,8 @@ Foam::word Foam::Time::findInstance
             if (debug)
             {
                 Info<< "Time::findInstance"
-                    "(const fileName&,const word&) : "
-                    << "found \"" << name
+                    "(const fileName&, const word&, const IOobject::readOption)"
+                    << " : found \"" << name
                     << "\" in " << ts[instanceI].name()/dir
                     << endl;
             }
@@ -129,8 +130,8 @@ Foam::word Foam::Time::findInstance
         if (debug)
         {
             Info<< "Time::findInstance"
-                   "(const fileName&,const word&) : "
-                << "found \"" << name
+                "(const fileName&, const word&, const IOobject::readOption)"
+                << " : found \"" << name
                 << "\" in " << constant()/dir
                 << endl;
         }
@@ -141,10 +142,10 @@ Foam::word Foam::Time::findInstance
     if (rOpt == IOobject::MUST_READ)
     {
         FatalErrorIn
-            (
-                "Time::findInstance(const fileName&,const word&)"
-            )
-            << "Cannot find file \"" << name << "\" in directory "
+        (
+            "Time::findInstance"
+            "(const fileName&, const word&, const IOobject::readOption)"
+        )   << "Cannot find file \"" << name << "\" in directory "
             << constant()/dir
             << exit(FatalError);
     }
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C b/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C
index 70f738fcfb..fa6a2a58d0 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C
@@ -154,7 +154,7 @@ tmp<GeometricField<scalar, PatchField, GeoMesh> > pow
             pow
             (
                 gsf1.dimensions(),
-                dimensionedScalar("1", 1.0, gsf2.dimensions())
+                dimensionedScalar("1", gsf2.dimensions(), 1.0)
             )
         )
     );
@@ -183,7 +183,7 @@ tmp<GeometricField<scalar, PatchField, GeoMesh> > pow
             pow
             (
                 gsf1.dimensions(),
-                dimensionedScalar("1", 1.0, gsf2.dimensions())
+                dimensionedScalar("1", gsf2.dimensions(), 1.0)
             )
         )
     );
@@ -214,7 +214,7 @@ tmp<GeometricField<scalar, PatchField, GeoMesh> > pow
             pow
             (
                 gsf1.dimensions(),
-                dimensionedScalar("1", 1.0, gsf2.dimensions())
+                dimensionedScalar("1", gsf2.dimensions(), 1.0)
             )
         )
     );
@@ -247,7 +247,7 @@ tmp<GeometricField<scalar, PatchField, GeoMesh> > pow
             pow
             (
                 gsf1.dimensions(),
-                dimensionedScalar("1", 1.0, gsf2.dimensions())
+                dimensionedScalar("1", gsf2.dimensions(), 1.0)
             )
         )
     );
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C
index 666d84a5c1..05958edf68 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C
+++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C
@@ -60,7 +60,7 @@ void Foam::lduMatrix::Amul
         interfaceBouCoeffs,
         interfaces,
         psi,
-        Apsi, 
+        Apsi,
         cmpt
     );
 
@@ -79,7 +79,7 @@ void Foam::lduMatrix::Amul
 
     register const label nFaces = upper().size();
     #ifdef ICC_IA64_PREFETCH
-    #pragma swp  
+    #pragma swp
     #endif
     for (register label face=0; face<nFaces; face++)
     {
@@ -227,7 +227,7 @@ void Foam::lduMatrix::sumA
     }
 
     #ifdef ICC_IA64_PREFETCH
-    #pragma swp  
+    #pragma swp
     #endif
 
     for (register label face=0; face<nFaces; face++)
@@ -316,7 +316,7 @@ void Foam::lduMatrix::residual
         mBouCoeffs,
         interfaces,
         psi,
-        rA, 
+        rA,
         cmpt
     );
 
@@ -336,7 +336,7 @@ void Foam::lduMatrix::residual
 
     register const label nFaces = upper().size();
     #ifdef ICC_IA64_PREFETCH
-    #pragma swp  
+    #pragma swp
     #endif
     for (register label face=0; face<nFaces; face++)
     {
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
index 0428b252d7..e0b4c0536c 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
@@ -54,40 +54,79 @@ void Foam::polyMesh::calcDirections() const
 {
     for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
     {
-        directions_[cmpt] = 1;
+        solutionD_[cmpt] = 1;
     }
 
-    label nEmptyPatches = 0;
+    // Knock out empty and wedge directions. Note:they will be present on all
+    // domains.
 
-    vector dirVec = vector::zero;
+    label nEmptyPatches = 0;
+    label nWedgePatches = 0;
+
+    vector emptyDirVec = vector::zero;
+    vector wedgeDirVec = vector::zero;
 
     forAll(boundaryMesh(), patchi)
     {
-        if (isA<emptyPolyPatch>(boundaryMesh()[patchi]))
+        if (boundaryMesh()[patchi].size())
         {
-            if (boundaryMesh()[patchi].size())
+            if (isA<emptyPolyPatch>(boundaryMesh()[patchi]))
             {
                 nEmptyPatches++;
-                dirVec += sum(cmptMag(boundaryMesh()[patchi].faceAreas()));
+                emptyDirVec += sum(cmptMag(boundaryMesh()[patchi].faceAreas()));
+            }
+            else if (isA<wedgePolyPatch>(boundaryMesh()[patchi]))
+            {
+                const wedgePolyPatch& wpp = refCast<const wedgePolyPatch>
+                (
+                    boundaryMesh()[patchi]
+                );
+
+                nWedgePatches++;
+                wedgeDirVec += cmptMag(wpp.centreNormal());
             }
         }
     }
 
     if (nEmptyPatches)
     {
-        reduce(dirVec, sumOp<vector>());
+        reduce(emptyDirVec, sumOp<vector>());
 
-        dirVec /= mag(dirVec);
+        emptyDirVec /= mag(emptyDirVec);
 
         for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
         {
-            if (dirVec[cmpt] > 1e-6)
+            if (emptyDirVec[cmpt] > 1e-6)
             {
-                directions_[cmpt] = -1;
+                solutionD_[cmpt] = -1;
             }
             else
             {
-                directions_[cmpt] = 1;
+                solutionD_[cmpt] = 1;
+            }
+        }
+    }
+
+
+    // Knock out wedge directions
+
+    geometricD_ = solutionD_;
+
+    if (nWedgePatches)
+    {
+        reduce(wedgeDirVec, sumOp<vector>());
+
+        wedgeDirVec /= mag(wedgeDirVec);
+
+        for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+        {
+            if (wedgeDirVec[cmpt] > 1e-6)
+            {
+                geometricD_[cmpt] = -1;
+            }
+            else
+            {
+                geometricD_[cmpt] = 1;
             }
         }
     }
@@ -163,7 +202,8 @@ Foam::polyMesh::polyMesh(const IOobject& io)
         *this
     ),
     bounds_(points_),
-    directions_(Vector<label>::zero),
+    geometricD_(Vector<label>::zero),
+    solutionD_(Vector<label>::zero),
     pointZones_
     (
         IOobject
@@ -350,7 +390,8 @@ Foam::polyMesh::polyMesh
         0
     ),
     bounds_(points_, syncPar),
-    directions_(Vector<label>::zero),
+    geometricD_(Vector<label>::zero),
+    solutionD_(Vector<label>::zero),
     pointZones_
     (
         IOobject
@@ -505,7 +546,8 @@ Foam::polyMesh::polyMesh
         0
     ),
     bounds_(points_, syncPar),
-    directions_(Vector<label>::zero),
+    geometricD_(Vector<label>::zero),
+    solutionD_(Vector<label>::zero),
     pointZones_
     (
         IOobject
@@ -766,44 +808,37 @@ const Foam::fileName& Foam::polyMesh::facesInstance() const
 }
 
 
-const Foam::Vector<Foam::label>& Foam::polyMesh::directions() const
+const Foam::Vector<Foam::label>& Foam::polyMesh::geometricD() const
 {
-    if (directions_.x() == 0)
+    if (geometricD_.x() == 0)
     {
         calcDirections();
     }
 
-    return directions_;
+    return geometricD_;
 }
 
 
 Foam::label Foam::polyMesh::nGeometricD() const
 {
-    label nWedges = 0;
+    return cmptSum(geometricD() + Vector<label>::one)/2;
+}
 
-    forAll(boundary_, patchi)
+
+const Foam::Vector<Foam::label>& Foam::polyMesh::solutionD() const
+{
+    if (solutionD_.x() == 0)
     {
-        if (isA<wedgePolyPatch>(boundary_[patchi]))
-        {
-            nWedges++;
-        }
+        calcDirections();
     }
 
-    if (nWedges != 0 && nWedges != 2 && nWedges != 4)
-    {
-        FatalErrorIn("label polyMesh::nGeometricD() const")
-            << "Number of wedge patches " << nWedges << " is incorrect, "
-               "should be 0, 2 or 4"
-            << exit(FatalError);
-    }
-
-    return nSolutionD() - nWedges/2;
+    return solutionD_;
 }
 
 
 Foam::label Foam::polyMesh::nSolutionD() const
 {
-    return cmptSum(directions() + Vector<label>::one)/2;
+    return cmptSum(solutionD() + Vector<label>::one)/2;
 }
 
 
@@ -823,6 +858,10 @@ void Foam::polyMesh::addPatches
             << abort(FatalError);
     }
 
+    // Reset valid directions
+    geometricD_ = Vector<label>::zero;
+    solutionD_ = Vector<label>::zero;
+
     boundary_.setSize(p.size());
 
     // Copy the patch pointers
@@ -1037,6 +1076,10 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
     faceZones_.movePoints(points_);
     cellZones_.movePoints(points_);
 
+    // Reset valid directions (could change with rotation)
+    geometricD_ = Vector<label>::zero;
+    solutionD_ = Vector<label>::zero;
+
 
     // Hack until proper callbacks. Below are all the polyMeh MeshObjects with a
     // movePoints function.
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.H b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
index 588394d1c2..ecb805ad7f 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
@@ -120,9 +120,13 @@ private:
             //  Created from points on construction, updated when the mesh moves
             boundBox bounds_;
 
+            //- vector of non-constrained directions in mesh
+            //  defined according to the presence of empty and wedge patches
+            mutable Vector<label> geometricD_;
+
             //- vector of valid directions in mesh
             //  defined according to the presence of empty patches
-            mutable Vector<label> directions_;
+            mutable Vector<label> solutionD_;
 
 
         // Zoning information
@@ -309,17 +313,22 @@ public:
                 return bounds_;
             }
 
-            //- Return the vector of valid directions in mesh.
-            //  Defined according to the presence of empty patches.
-            //  1 indicates valid direction and -1 an invalid direction.
-            const Vector<label>& directions() const;
+            //- Return the vector of geometric directions in mesh.
+            //  Defined according to the presence of empty and wedge patches.
+            //  1 indicates unconstrained direction and -1 a constrained
+            //  direction.
+            const Vector<label>& geometricD() const;
 
             //- Return the number of valid geometric dimensions in the mesh
             label nGeometricD() const;
 
-            //- Return the number of valid solution dimensions in the mesh.
-            //  For wedge cases this includes the circumferential direction
-            //  in case of swirl.
+            //- Return the vector of solved-for directions in mesh.
+            //  Differs from geometricD in that it includes for wedge cases
+            //  the circumferential direction in case of swirl.
+            //  1 indicates valid direction and -1 an invalid direction.
+            const Vector<label>& solutionD() const;
+
+            //- Return the number of valid solved-for dimensions in the mesh
             label nSolutionD() const;
 
             //- Return point zone mesh
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C b/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
index 52f4446444..7183120984 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
@@ -68,6 +68,10 @@ void Foam::polyMesh::clearGeom()
         boundary_[patchI].clearGeom();
     }
 
+    // Reset valid directions (could change with rotation)
+    geometricD_ = Vector<label>::zero;
+    solutionD_ = Vector<label>::zero;
+
     pointMesh::Delete(*this);
 }
 
@@ -87,6 +91,10 @@ void Foam::polyMesh::clearAddressing()
     // recalculation
     deleteDemandDrivenData(globalMeshDataPtr_);
 
+    // Reset valid directions
+    geometricD_ = Vector<label>::zero;
+    solutionD_ = Vector<label>::zero;
+
     pointMesh::Delete(*this);
 }
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C
index d12dffff62..f1cc51cb9e 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C
@@ -217,7 +217,8 @@ Foam::polyMesh::polyMesh
         boundaryFaces.size() + 1    // add room for a default patch
     ),
     bounds_(points_, syncPar),
-    directions_(Vector<label>::zero),
+    geometricD_(Vector<label>::zero),
+    solutionD_(Vector<label>::zero),
     pointZones_
     (
         IOobject
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C b/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C
index 24d9ebbb13..7dcd3ba6b4 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C
@@ -68,6 +68,11 @@ void Foam::polyMesh::updateMesh(const mapPolyMesh& mpm)
         newMotionPoints.map(oldMotionPoints, mpm.pointMap());
     }
 
+    // Reset valid directions (could change by faces put into empty patches)
+    geometricD_ = Vector<label>::zero;
+    solutionD_ = Vector<label>::zero;
+
+
     // Hack until proper callbacks. Below are all the polyMesh-MeshObjects.
 
     // pointMesh
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
index ecfff62ac3..ddf162dd0f 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
@@ -29,6 +29,7 @@ License
 #include "mathematicalConstants.H"
 #include "refinementSurfaces.H"
 #include "searchableSurfaces.H"
+#include "regExp.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -300,10 +301,44 @@ Foam::layerParameters::layerParameters
             //    readScalar(layerDict.lookup("minThickness"));
         }
     }
+
+
+    // Check whether layer specification matches any patches
+    const List<keyType> wildCards = layersDict.keys(true);
+
+    forAll(wildCards, i)
+    {
+        regExp re(wildCards[i]);
+
+        bool hasMatch = false;
+        forAll(boundaryMesh, patchI)
+        {
+            if (re.match(boundaryMesh[patchI].name()))
+            {
+                hasMatch = true;
+                break;
+            }
+        }
+        if (!hasMatch)
+        {
+            IOWarningIn("layerParameters::layerParameters(..)", layersDict)
+                << "Wildcard layer specification for " << wildCards[i]
+                << " does not match any patch." << endl;
+        }
+    }
+
+    const List<keyType> nonWildCards = layersDict.keys(false);
+
+    forAll(nonWildCards, i)
+    {
+        if (boundaryMesh.findPatchID(nonWildCards[i]) == -1)
+        {
+            IOWarningIn("layerParameters::layerParameters(..)", layersDict)
+                << "Layer specification for " << nonWildCards[i]
+                << " does not match any patch." << endl;
+        }
+    }
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //
diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index a626103cd1..6bbb1c8d39 100644
--- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -265,101 +265,112 @@ Foam::refinementSurfaces::refinementSurfaces
     zoneInside_(surfacesDict.size()),
     regionOffset_(surfacesDict.size())
 {
-    labelList globalMinLevel(surfacesDict.size(), 0);
-    labelList globalMaxLevel(surfacesDict.size(), 0);
-    scalarField globalAngle(surfacesDict.size(), -GREAT);
-    List<Map<label> > regionMinLevel(surfacesDict.size());
-    List<Map<label> > regionMaxLevel(surfacesDict.size());
-    List<Map<scalar> > regionAngle(surfacesDict.size());
+    // Wilcard specification : loop over all surface, all regions
+    // and try to find a match.
 
+    // Count number of surfaces.
     label surfI = 0;
-    forAllConstIter(dictionary, surfacesDict, iter)
+    forAll(allGeometry.names(), geomI)
     {
-        names_[surfI] = iter().keyword();
+        const word& geomName = allGeometry_.names()[geomI];
 
-        surfaces_[surfI] = allGeometry_.findSurfaceID(names_[surfI]);
-
-        if (surfaces_[surfI] == -1)
+        if (surfacesDict.found(geomName))
         {
-            FatalErrorIn
-            (
-                "refinementSurfaces::refinementSurfaces"
-                "(const searchableSurfaces&, const dictionary>&"
-            )   << "No surface called " << iter().keyword() << endl
-                << "Valid surfaces are " << allGeometry_.names()
-                << exit(FatalError);
+            surfI++;
         }
-        const dictionary& dict = surfacesDict.subDict(iter().keyword());
+    }
 
-        const labelPair refLevel(dict.lookup("level"));
-        globalMinLevel[surfI] = refLevel[0];
-        globalMaxLevel[surfI] = refLevel[1];
+    // Size lists
+    surfaces_.setSize(surfI);
+    names_.setSize(surfI);
+    faceZoneNames_.setSize(surfI);
+    cellZoneNames_.setSize(surfI);
+    zoneInside_.setSize(surfI);
+    regionOffset_.setSize(surfI);
 
-        // Global zone names per surface
-        if (dict.found("faceZone"))
+    labelList globalMinLevel(surfI, 0);
+    labelList globalMaxLevel(surfI, 0);
+    scalarField globalAngle(surfI, -GREAT);
+    List<Map<label> > regionMinLevel(surfI);
+    List<Map<label> > regionMaxLevel(surfI);
+    List<Map<scalar> > regionAngle(surfI);
+
+    surfI = 0;
+    forAll(allGeometry.names(), geomI)
+    {
+        const word& geomName = allGeometry_.names()[geomI];
+
+        if (surfacesDict.found(geomName))
         {
-            dict.lookup("faceZone") >> faceZoneNames_[surfI];
-            dict.lookup("cellZone") >> cellZoneNames_[surfI];
-            dict.lookup("zoneInside") >> zoneInside_[surfI];
-        }
+            const dictionary& dict = surfacesDict.subDict(geomName);
 
-        // Global perpendicular angle
-        if (dict.found("perpendicularAngle"))
-        {
-            globalAngle[surfI] = readScalar(dict.lookup("perpendicularAngle"));
-        }
+            names_[surfI] = geomName;
+            surfaces_[surfI] = geomI;
 
-        if (dict.found("regions"))
-        {
-            const dictionary& regionsDict = dict.subDict("regions");
-            const wordList& regionNames =
-                allGeometry_[surfaces_[surfI]].regions();
+            const labelPair refLevel(dict.lookup("level"));
+            globalMinLevel[surfI] = refLevel[0];
+            globalMaxLevel[surfI] = refLevel[1];
 
-            forAllConstIter(dictionary, regionsDict, iter)
+            // Global zone names per surface
+            if (dict.found("faceZone"))
             {
-                const word& key = iter().keyword();
+                dict.lookup("faceZone") >> faceZoneNames_[surfI];
+                dict.lookup("cellZone") >> cellZoneNames_[surfI];
+                dict.lookup("zoneInside") >> zoneInside_[surfI];
+            }
 
-                if (regionsDict.isDict(key))
+            // Global perpendicular angle
+            if (dict.found("perpendicularAngle"))
+            {
+                globalAngle[surfI] = readScalar
+                (
+                    dict.lookup("perpendicularAngle")
+                );
+            }
+
+            if (dict.found("regions"))
+            {
+                const dictionary& regionsDict = dict.subDict("regions");
+                const wordList& regionNames =
+                    allGeometry_[surfaces_[surfI]].regions();
+
+                forAll(regionNames, regionI)
                 {
-                    // Get the dictionary for region iter.keyword()
-                    const dictionary& regionDict = regionsDict.subDict(key);
-
-                    label regionI = findIndex(regionNames, key);
-                    if (regionI == -1)
+                    if (regionsDict.found(regionNames[regionI]))
                     {
-                        FatalErrorIn
+                        // Get the dictionary for region 
+                        const dictionary& regionDict = regionsDict.subDict
                         (
-                            "refinementSurfaces::refinementSurfaces"
-                            "(const searchableSurfaces&, const dictionary>&"
-                        )   << "No region called " << key << " on surface "
-                            << allGeometry_[surfaces_[surfI]].name() << endl
-                            << "Valid regions are " << regionNames
-                            << exit(FatalError);
-                    }
-
-                    const labelPair refLevel(regionDict.lookup("level"));
-
-                    regionMinLevel[surfI].insert(regionI, refLevel[0]);
-                    regionMaxLevel[surfI].insert(regionI, refLevel[1]);
-
-                    if (regionDict.found("perpendicularAngle"))
-                    {
-                        regionAngle[surfI].insert
-                        (
-                            regionI,
-                            readScalar(regionDict.lookup("perpendicularAngle"))
+                            regionNames[regionI]
                         );
+
+                        const labelPair refLevel(regionDict.lookup("level"));
+
+                        regionMinLevel[surfI].insert(regionI, refLevel[0]);
+                        regionMaxLevel[surfI].insert(regionI, refLevel[1]);
+
+                        if (regionDict.found("perpendicularAngle"))
+                        {
+                            regionAngle[surfI].insert
+                            (
+                                regionI,
+                                readScalar
+                                (
+                                    regionDict.lookup("perpendicularAngle")
+                                )
+                            );
+                        }
                     }
                 }
             }
+            surfI++;
         }
-        surfI++;
     }
 
     // Calculate local to global region offset
     label nRegions = 0;
 
-    forAll(surfacesDict, surfI)
+    forAll(surfaces_, surfI)
     {
         regionOffset_[surfI] = nRegions;
         nRegions += allGeometry_[surfaces_[surfI]].regions().size();
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
index 8bc4d3b938..ab9709e993 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
@@ -29,6 +29,10 @@ License
 #include "volFields.H"
 #include "surfaceFields.H"
 #include "fvMatrices.H"
+#include "PackedList.H"
+#include "syncTools.H"
+
+#include "faceSet.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -39,12 +43,15 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
     dict_(is),
     cellZoneID_(mesh_.cellZones().findZoneID(name_)),
     faceZoneID_(mesh_.faceZones().findZoneID(name_)),
+    outsideFaces_(0),
     patchNames_(dict_.lookup("patches")),
     origin_(dict_.lookup("origin")),
     axis_(dict_.lookup("axis")),
     omega_(dict_.lookup("omega")),
     Omega_("Omega", omega_*axis_)
 {
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
     axis_ = axis_/mag(axis_);
     Omega_ = omega_*axis_;
 
@@ -65,18 +72,71 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
 
     if (!faceZoneFound)
     {
-        FatalErrorIn
-        (
-            "Foam::MRFZone::MRFZone(const fvMesh& , const dictionary&)"
-        )   << "cannot find MRF faceZone " << name_
-            << exit(FatalError);
+        // Determine faces in cell zone
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // (does not construct cells)
+
+        const labelList& own = mesh_.faceOwner();
+        const labelList& nei = mesh_.faceNeighbour();
+
+        // Cells in zone
+        PackedBoolList zoneCell(mesh_.nCells());
+
+        if (cellZoneID_ != -1)
+        {
+            const labelList& cellLabels = mesh_.cellZones()[cellZoneID_];
+            forAll(cellLabels, i)
+            {
+                zoneCell[cellLabels[i]] = 1u;
+            }
+        }
+
+
+        // Faces in zone
+        PackedBoolList zoneFacesSet(mesh_.nFaces());
+
+        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+        {
+            if
+            (
+                zoneCell.get(own[faceI]) == 1u
+             || zoneCell.get(nei[faceI]) == 1u
+            )
+            {
+                zoneFacesSet[faceI] = 1u;
+            }
+        }
+        syncTools::syncFaceList(mesh_, zoneFacesSet, orEqOp<unsigned int>());
+
+
+        // Transfer to labelList
+        label n = zoneFacesSet.count();
+        outsideFaces_.setSize(n);
+        n = 0;
+        forAll(zoneFacesSet, faceI)
+        {
+            if (zoneFacesSet.get(faceI) == 1u)
+            {
+                outsideFaces_[n++] = faceI;
+            }
+        }
+
+        Info<< nl
+            << "MRFZone " << name_ << " : did not find a faceZone; using "
+            << returnReduce(outsideFaces_.size(), sumOp<label>())
+            << " faces internal to cellZone instead." << endl;
+
+
+        // Flag use of outsideFaces
+        faceZoneID_ = -2;
     }
 
+
     patchLabels_.setSize(patchNames_.size());
 
     forAll(patchNames_, i)
     {
-        patchLabels_[i] = mesh_.boundaryMesh().findPatchID(patchNames_[i]);
+        patchLabels_[i] = patches.findPatchID(patchNames_[i]);
 
         if (patchLabels_[i] == -1)
         {
@@ -125,7 +185,13 @@ void Foam::MRFZone::relativeFlux(surfaceScalarField& phi) const
     const vector& origin = origin_.value();
     const vector& Omega = Omega_.value();
 
-    const labelList& faces = mesh_.faceZones()[faceZoneID_];
+    // Use either zone faces or outsideFaces_
+    const labelList& faces =
+    (
+        faceZoneID_ == -2
+      ? outsideFaces_
+      : mesh_.faceZones()[faceZoneID_]
+    );
 
     forAll(faces, i)
     {
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H b/src/finiteVolume/cfdTools/general/MRF/MRFZone.H
index 3ffb4e461a..7c489ed25b 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.H
@@ -26,7 +26,7 @@ Class
     Foam::MRFZone
 
 Description
-    MRF zone definition based on both cell and face zones and parameters
+    MRF zone definition based on cell zone and optional face zone and parameters
     obtained from a control dictionary constructed from the given stream.
 
     The rotation of the MRF region is defined by an origin and axis of
@@ -68,18 +68,26 @@ class MRFZone
 
         const fvMesh& mesh_;
 
-        word name_;
+        const word name_;
 
-        dictionary dict_;
+        const dictionary dict_;
 
         label cellZoneID_;
+
+        //- label of face zone with faces on outside of cell zone.
+        //  If -2 determines outside faces itself
         label faceZoneID_;
-        wordList patchNames_;
+
+        //- outside faces (only if faceZoneID = -2)
+        labelList outsideFaces_;
+
+
+        const wordList patchNames_;
         labelList patchLabels_;
 
-        dimensionedVector origin_;
+        const dimensionedVector origin_;
         dimensionedVector axis_;
-        dimensionedScalar omega_;
+        const dimensionedScalar omega_;
         dimensionedVector Omega_;
 
 
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C
index 133f1e7933..6bef8de896 100644
--- a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C
@@ -140,11 +140,11 @@ void Foam::quadraticFitSnGradData::findFaceDirs
     #ifndef SPHERICAL_GEOMETRY
         if (mesh.nGeometricD() <= 2) // find the normal direcion
         {
-            if (mesh.directions()[0] == -1)
+            if (mesh.geometricD()[0] == -1)
             {
                 kdir = vector(1, 0, 0);
             }
-            else if (mesh.directions()[1] == -1)
+            else if (mesh.geometricD()[1] == -1)
             {
                 kdir = vector(0, 1, 0);
             }
@@ -153,7 +153,7 @@ void Foam::quadraticFitSnGradData::findFaceDirs
                 kdir = vector(0, 0, 1);
             }
         }
-        else // 3D so find a direction in the place of the face
+        else // 3D so find a direction in the plane of the face
         {
             const face& f = mesh.faces()[faci];
             kdir = mesh.points()[f[0]] - mesh.points()[f[1]];
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index e1ef6135a6..2134420ead 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -713,7 +713,7 @@ Foam::fvMatrix<Type>::H() const
     (
         pow
         (
-            psi_.mesh().directions(),
+            psi_.mesh().solutionD(),
             pTraits<typename powProduct<Vector<label>, Type::rank>::type>::zero
         )
     );
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
index e8f4e6de49..21f3738a02 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
@@ -82,7 +82,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
     (
         pow
         (
-            psi_.mesh().directions(),
+            psi_.mesh().solutionD(),
             pTraits<typename powProduct<Vector<label>, Type::rank>::type>::zero
         )
     );
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C
index 0872cb81a9..28a3c41e64 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C
@@ -83,11 +83,11 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::findFaceDirs
 #   ifndef SPHERICAL_GEOMETRY
     if (mesh.nGeometricD() <= 2) // find the normal direction
     {
-        if (mesh.directions()[0] == -1)
+        if (mesh.geometricD()[0] == -1)
         {
             kdir = vector(1, 0, 0);
         }
-        else if (mesh.directions()[1] == -1)
+        else if (mesh.geometricD()[1] == -1)
         {
             kdir = vector(0, 1, 0);
         }
@@ -115,7 +115,7 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::findFaceDirs
 
         if (magk < SMALL)
         {
-            FatalErrorIn("findFaceDirs") << " calculated kdir = zero"
+            FatalErrorIn("findFaceDirs(..)") << " calculated kdir = zero"
                 << exit(FatalError);
         }
         else
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index a93cd58400..9f849c5894 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -60,6 +60,7 @@ $(searchableSurface)/searchablePlane.C
 $(searchableSurface)/searchablePlate.C
 $(searchableSurface)/searchableSphere.C
 $(searchableSurface)/searchableSurface.C
+$(searchableSurface)/searchableSurfaceCollection.C
 $(searchableSurface)/searchableSurfaces.C
 $(searchableSurface)/searchableSurfacesQueries.C
 $(searchableSurface)/searchableSurfaceWithGaps.C
diff --git a/src/meshTools/cellDist/wallPoint/wallPoint.H b/src/meshTools/cellDist/wallPoint/wallPoint.H
index 0544b6ce36..82d6b385af 100644
--- a/src/meshTools/cellDist/wallPoint/wallPoint.H
+++ b/src/meshTools/cellDist/wallPoint/wallPoint.H
@@ -51,6 +51,12 @@ namespace Foam
 // Forward declaration of classes
 class polyPatch;
 class polyMesh;
+class wallPoint;
+
+// Forward declaration of friend functions and operators
+Ostream& operator<<(Ostream&, const wallPoint&);
+Istream& operator>>(Istream&, wallPoint&);
+
 
 /*---------------------------------------------------------------------------*\
                            Class wallPoint Declaration
@@ -78,12 +84,15 @@ class wallPoint
             const scalar tol
         );
 
+
 public:
+
     // Static data members
 
         //- initial point far away.
         static point greatPoint;
 
+
     // Constructors
 
         //- Construct null
@@ -102,6 +111,7 @@ public:
             const wallPoint&
         );
 
+
     // Member Functions
 
         // Access
@@ -184,11 +194,11 @@ public:
                 const scalar tol
             );
 
+
     // Member Operators
 
         // Needed for List IO
         inline bool operator==(const wallPoint&) const;
-
         inline bool operator!=(const wallPoint&) const;
 
 
diff --git a/src/meshTools/meshTools/meshTools.C b/src/meshTools/meshTools/meshTools.C
index a4f57a8de8..f826675b30 100644
--- a/src/meshTools/meshTools/meshTools.C
+++ b/src/meshTools/meshTools/meshTools.C
@@ -25,7 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "meshTools.H"
-#include "primitiveMesh.H"
+#include "polyMesh.H"
 #include "hexMatcher.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -635,6 +635,103 @@ Foam::label Foam::meshTools::walkFace
 }
 
 
+void Foam::meshTools::constrainToMeshCentre(const polyMesh& mesh, point& pt)
+{
+    const Vector<label>& dirs = mesh.geometricD();
+    const point& min = mesh.bounds().min();
+    const point& max = mesh.bounds().max();
+
+    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+    {
+        if (dirs[cmpt] == -1)
+        {
+            pt[cmpt] = 0.5*(min[cmpt]+max[cmpt]);
+        }
+    }
+}
+
+
+void Foam::meshTools::constrainToMeshCentre
+(
+    const polyMesh& mesh,
+    pointField& pts
+)
+{
+    const Vector<label>& dirs = mesh.geometricD();
+    const point& min = mesh.bounds().min();
+    const point& max = mesh.bounds().max();
+
+    bool isConstrained = false;
+    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+    {
+        if (dirs[cmpt] == -1)
+        {
+            isConstrained = true;
+            break;
+        }
+    }
+
+    if (isConstrained)
+    {
+        forAll(pts, i)
+        {
+            for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+            {
+                if (dirs[cmpt] == -1)
+                {
+                    pts[i][cmpt] = 0.5*(min[cmpt]+max[cmpt]);
+                }
+            }
+        }
+    }
+}
+
+
+//- Set the constrained components of directions/velocity to zero
+void Foam::meshTools::constrainDirection(const polyMesh& mesh, vector& d)
+{
+    const Vector<label>& dirs = mesh.geometricD();
+
+    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+    {
+        if (dirs[cmpt] == -1)
+        {
+            d[cmpt] = 0.0;
+        }
+    }
+}
+
+
+void Foam::meshTools::constrainDirection(const polyMesh& mesh, vectorField& d)
+{
+    const Vector<label>& dirs = mesh.geometricD();
+
+    bool isConstrained = false;
+    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+    {
+        if (dirs[cmpt] == -1)
+        {
+            isConstrained = true;
+            break;
+        }
+    }
+
+    if (isConstrained)
+    {
+        forAll(d, i)
+        {
+            for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+            {
+                if (dirs[cmpt] == -1)
+                {
+                    d[i][cmpt] = 0.0;
+                }
+            }
+        }
+    }
+}
+
+
 void Foam::meshTools::getParallelEdges
 (
     const primitiveMesh& mesh,
diff --git a/src/meshTools/meshTools/meshTools.H b/src/meshTools/meshTools/meshTools.H
index 304159e045..3710d5af06 100644
--- a/src/meshTools/meshTools/meshTools.H
+++ b/src/meshTools/meshTools/meshTools.H
@@ -50,7 +50,7 @@ namespace Foam
 {
 
 class primitiveMesh;
-
+class polyMesh;
 
 /*---------------------------------------------------------------------------*\
                         Namespace meshTools Declaration
@@ -81,196 +81,212 @@ namespace meshTools
     static const label pXpYpZMask = 1 << pXpYpZ;
 
 
-    //- Check if n is in same direction as normals of all faceLabels
-    bool visNormal
-    (
-        const vector& n,
-        const vectorField& faceNormals,
-        const labelList& faceLabels
-    );
+    // Normal handling
 
-    //- Calculate point normals on a 'box' mesh (all edges aligned with
-    //  coordinate axes)
-    vectorField calcBoxPointNormals(const primitivePatch& pp);
+        //- Check if n is in same direction as normals of all faceLabels
+        bool visNormal
+        (
+            const vector& n,
+            const vectorField& faceNormals,
+            const labelList& faceLabels
+        );
 
-    //- Normalized edge vector
-    vector normEdgeVec(const primitiveMesh&, const label edgeI);
+        //- Calculate point normals on a 'box' mesh (all edges aligned with
+        //  coordinate axes)
+        vectorField calcBoxPointNormals(const primitivePatch& pp);
 
-    //- Write obj representation of point
-    void writeOBJ
-    (
-        Ostream& os,
-        const point& pt
-    );
+        //- Normalized edge vector
+        vector normEdgeVec(const primitiveMesh&, const label edgeI);
 
 
-    //- Write obj representation of faces subset
-    void writeOBJ
-    (
-        Ostream& os,
-        const faceList&,
-        const pointField&,
-        const labelList& faceLabels
-    );
+    // OBJ writing
 
-    //- Write obj representation of faces
-    void writeOBJ
-    (
-        Ostream& os,
-        const faceList&,
-        const pointField&
-    );
+        //- Write obj representation of point
+        void writeOBJ
+        (
+            Ostream& os,
+            const point& pt
+        );
 
-    //- Write obj representation of cell subset
-    void writeOBJ
-    (
-        Ostream& os,
-        const cellList&,
-        const faceList&,
-        const pointField&,
-        const labelList& cellLabels
-    );
+        //- Write obj representation of faces subset
+        void writeOBJ
+        (
+            Ostream& os,
+            const faceList&,
+            const pointField&,
+            const labelList& faceLabels
+        );
 
-    //- Is edge used by cell
-    bool edgeOnCell
-    (
-        const primitiveMesh&,
-        const label cellI,
-        const label edgeI
-    );
+        //- Write obj representation of faces
+        void writeOBJ
+        (
+            Ostream& os,
+            const faceList&,
+            const pointField&
+        );
 
-    //- Is edge used by face
-    bool edgeOnFace
-    (
-        const primitiveMesh&,
-        const label faceI,
-        const label edgeI
-    );
-    
-    //- Is face used by cell
-    bool faceOnCell
-    (
-        const primitiveMesh&,
-        const label cellI,
-        const label faceI
-    );
+        //- Write obj representation of cell subset
+        void writeOBJ
+        (
+            Ostream& os,
+            const cellList&,
+            const faceList&,
+            const pointField&,
+            const labelList& cellLabels
+        );
 
-    //- Return edge among candidates that uses the two vertices.
-    label findEdge
-    (
-        const edgeList& edges,
-        const labelList& candidates,
-        const label v0,
-        const label v1
-    );
 
-    //- Return edge between two vertices. Returns -1 if no edge.
-    label findEdge
-    (
-        const primitiveMesh&,
-        const label v0,
-        const label v1
-    );
+    // Cell/face/edge walking
 
-    //- Return edge shared by two faces. Throws error if no edge found.
-    label getSharedEdge
-    (
-        const primitiveMesh&,
-        const label f0,
-        const label f1
-    );
-    
-    //- Return face shared by two cells. Throws error if none found.
-    label getSharedFace
-    (
-        const primitiveMesh&,
-        const label cell0,
-        const label cell1
-    );
+        //- Is edge used by cell
+        bool edgeOnCell
+        (
+            const primitiveMesh&,
+            const label cellI,
+            const label edgeI
+        );
 
-    //- Get faces on cell using edgeI. Throws error if no two found.
-    void getEdgeFaces
-    (
-        const primitiveMesh&,
-        const label cellI,
-        const label edgeI,
-        label& face0,
-        label& face1
-    );
+        //- Is edge used by face
+        bool edgeOnFace
+        (
+            const primitiveMesh&,
+            const label faceI,
+            const label edgeI
+        );
 
-    //- Return label of other edge (out of candidates edgeLabels)
-    //  connected to vertex but not edgeI. Throws error if none found.
-    label otherEdge
-    (
-        const primitiveMesh&,
-        const labelList& edgeLabels,
-        const label edgeI,
-        const label vertI
-    );
-    
-    //- Return face on cell using edgeI but not faceI. Throws error 
-    //  if none found.
-    label otherFace
-    (
-        const primitiveMesh&,
-        const label cellI,
-        const label faceI,
-        const label edgeI
-    );
+        //- Is face used by cell
+        bool faceOnCell
+        (
+            const primitiveMesh&,
+            const label cellI,
+            const label faceI
+        );
 
-    //- Return cell on other side of face. Throws error 
-    //  if face not internal.
-    label otherCell
-    (
-        const primitiveMesh&,
-        const label cellI,
-        const label faceI
-    );
-    
-    //- Returns label of edge nEdges away from startEdge (in the direction
-    // of startVertI)
-    label walkFace
-    (
-        const primitiveMesh&,
-        const label faceI,
-        const label startEdgeI,
-        const label startVertI,
-        const label nEdges
-    );
+        //- Return edge among candidates that uses the two vertices.
+        label findEdge
+        (
+            const edgeList& edges,
+            const labelList& candidates,
+            const label v0,
+            const label v1
+        );
+
+        //- Return edge between two vertices. Returns -1 if no edge.
+        label findEdge
+        (
+            const primitiveMesh&,
+            const label v0,
+            const label v1
+        );
+
+        //- Return edge shared by two faces. Throws error if no edge found.
+        label getSharedEdge
+        (
+            const primitiveMesh&,
+            const label f0,
+            const label f1
+        );
+
+        //- Return face shared by two cells. Throws error if none found.
+        label getSharedFace
+        (
+            const primitiveMesh&,
+            const label cell0,
+            const label cell1
+        );
+
+        //- Get faces on cell using edgeI. Throws error if no two found.
+        void getEdgeFaces
+        (
+            const primitiveMesh&,
+            const label cellI,
+            const label edgeI,
+            label& face0,
+            label& face1
+        );
+
+        //- Return label of other edge (out of candidates edgeLabels)
+        //  connected to vertex but not edgeI. Throws error if none found.
+        label otherEdge
+        (
+            const primitiveMesh&,
+            const labelList& edgeLabels,
+            const label edgeI,
+            const label vertI
+        );
+
+        //- Return face on cell using edgeI but not faceI. Throws error 
+        //  if none found.
+        label otherFace
+        (
+            const primitiveMesh&,
+            const label cellI,
+            const label faceI,
+            const label edgeI
+        );
+
+        //- Return cell on other side of face. Throws error 
+        //  if face not internal.
+        label otherCell
+        (
+            const primitiveMesh&,
+            const label cellI,
+            const label faceI
+        );
+
+        //- Returns label of edge nEdges away from startEdge (in the direction
+        // of startVertI)
+        label walkFace
+        (
+            const primitiveMesh&,
+            const label faceI,
+            const label startEdgeI,
+            const label startVertI,
+            const label nEdges
+        );
+
+
+    // Constraints on position
+
+        //- Set the constrained components of position to mesh centre
+        void constrainToMeshCentre(const polyMesh&, point&);
+        void constrainToMeshCentre(const polyMesh&, pointField&);
+
+        //- Set the constrained components of directions/velocity to zero
+        void constrainDirection(const polyMesh&, vector&);
+        void constrainDirection(const polyMesh&, vectorField&);
 
     
-    //
     // Hex only functionality.
-    //
 
-    //- Given edge on hex find other 'parallel', non-connected edges.
-    void getParallelEdges
-    (
-        const primitiveMesh&,
-        const label cellI,
-        const label e0,
-        label&,
-        label&,
-        label&
-    );
+        //- Given edge on hex find other 'parallel', non-connected edges.
+        void getParallelEdges
+        (
+            const primitiveMesh&,
+            const label cellI,
+            const label e0,
+            label&,
+            label&,
+            label&
+        );
 
-    //- Given edge on hex find all 'parallel' (i.e. non-connected)
-    //  edges and average direction of them
-    vector edgeToCutDir
-    (
-        const primitiveMesh&,
-        const label cellI,
-        const label edgeI
-    );
+        //- Given edge on hex find all 'parallel' (i.e. non-connected)
+        //  edges and average direction of them
+        vector edgeToCutDir
+        (
+            const primitiveMesh&,
+            const label cellI,
+            const label edgeI
+        );
 
-    //- Reverse of edgeToCutDir: given direction find edge bundle and
-    //  return one of them.
-    label cutDirToEdge
-    (
-        const primitiveMesh&,
-        const label cellI,
-        const vector& cutDir
-    );
+        //- Reverse of edgeToCutDir: given direction find edge bundle and
+        //  return one of them.
+        label cutDirToEdge
+        (
+            const primitiveMesh&,
+            const label cellI,
+            const vector& cutDir
+        );
     
 } // End namespace meshTools
 
diff --git a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
index 0602046703..7f86c0f658 100644
--- a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
+++ b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
@@ -76,10 +76,8 @@ bool Foam::distributedTriSurfaceMesh::read()
     // Distribution type
     distType_ = distributionTypeNames_.read(dict_.lookup("distributionType"));
 
-    if (dict_.found("mergeDistance"))
-    {
-        dict_.lookup("mergeDistance") >> mergeDist_;
-    }
+    // Merge distance
+    mergeDist_ = readScalar(dict_.lookup("mergeDistance"));
 
     return true;
 }
@@ -1341,25 +1339,44 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
     dict_(io, dict)
 {
     read();
+
+    if (debug)
+    {
+        Info<< "Constructed from triSurface:" << endl;
+        writeStats(Info);
+
+        labelList nTris(Pstream::nProcs());
+        nTris[Pstream::myProcNo()] = triSurface::size();
+        Pstream::gatherList(nTris);
+        Pstream::scatterList(nTris);
+
+        Info<< endl<< "\tproc\ttris\tbb" << endl;
+        forAll(nTris, procI)
+        {
+            Info<< '\t' << procI << '\t' << nTris[procI]
+                 << '\t' << procBb_[procI] << endl;
+        }
+        Info<< endl;
+    }
 }
 
 
 Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
 :
-    triSurfaceMesh(io),
-//    triSurfaceMesh
-//    (
-//        IOobject
-//        (
-//            io.name(),
-//            io.db().time().findInstanceDir(io.local()),
-//            io.local(),
-//            io.db(),
-//            io.readOpt(),
-//            io.writeOpt(),
-//            io.registerObject()
-//        )
-//    ),
+    //triSurfaceMesh(io),
+    triSurfaceMesh
+    (
+        IOobject
+        (
+            io.name(),
+            io.time().findInstance(io.local(), word::null),
+            io.local(),
+            io.db(),
+            io.readOpt(),
+            io.writeOpt(),
+            io.registerObject()
+        )
+    ),
     dict_
     (
         IOobject
@@ -1375,6 +1392,26 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
     )
 {
     read();
+
+    if (debug)
+    {
+        Info<< "Read distributedTriSurface from " << io.objectPath()
+            << ':' << endl;
+        writeStats(Info);
+
+        labelList nTris(Pstream::nProcs());
+        nTris[Pstream::myProcNo()] = triSurface::size();
+        Pstream::gatherList(nTris);
+        Pstream::scatterList(nTris);
+
+        Info<< endl<< "\tproc\ttris\tbb" << endl;
+        forAll(nTris, procI)
+        {
+            Info<< '\t' << procI << '\t' << nTris[procI]
+                 << '\t' << procBb_[procI] << endl;
+        }
+        Info<< endl;
+    }
 }
 
 
@@ -1384,21 +1421,21 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
     const dictionary& dict
 )
 :
-    triSurfaceMesh(io, dict),
-//    triSurfaceMesh
-//    (
-//        IOobject
-//        (
-//            io.name(),
-//            io.db().time().findInstanceDir(io.local()),
-//            io.local(),
-//            io.db(),
-//            io.readOpt(),
-//            io.writeOpt(),
-//            io.registerObject()
-//        ),
-//        dict
-//    ),
+    //triSurfaceMesh(io, dict),
+    triSurfaceMesh
+    (
+        IOobject
+        (
+            io.name(),
+            io.time().findInstance(io.local(), word::null),
+            io.local(),
+            io.db(),
+            io.readOpt(),
+            io.writeOpt(),
+            io.registerObject()
+        ),
+        dict
+    ),
     dict_
     (
         IOobject
@@ -1414,6 +1451,26 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
     )
 {
     read();
+
+    if (debug)
+    {
+        Info<< "Read distributedTriSurface from " << io.objectPath()
+            << " and dictionary:" << endl;
+        writeStats(Info);
+
+        labelList nTris(Pstream::nProcs());
+        nTris[Pstream::myProcNo()] = triSurface::size();
+        Pstream::gatherList(nTris);
+        Pstream::scatterList(nTris);
+
+        Info<< endl<< "\tproc\ttris\tbb" << endl;
+        forAll(nTris, procI)
+        {
+            Info<< '\t' << procI << '\t' << nTris[procI]
+                 << '\t' << procBb_[procI] << endl;
+        }
+        Info<< endl;
+    }
 }
 
 
@@ -2416,6 +2473,8 @@ void Foam::distributedTriSurfaceMesh::writeStats(Ostream& os) const
     boundBox bb;
     label nPoints;
     calcBounds(bb, nPoints);
+    reduce(bb.min(), minOp<point>());
+    reduce(bb.max(), maxOp<point>());
 
     os  << "Triangles    : " << returnReduce(triSurface::size(), sumOp<label>())
         << endl
diff --git a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.H b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.H
index 733cd24891..4d5f7dfd17 100644
--- a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.H
+++ b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.H
@@ -91,15 +91,14 @@ public:
 
 private:
 
+    // Private member data
+
         //- Merging distance
         scalar mergeDist_;
 
         //- Decomposition used when independently decomposing surface.
         autoPtr<decompositionMethod> decomposer_;
 
-
-    // Private member data
-
         //- Bounding box settings
         IOdictionary dict_;
 
@@ -317,10 +316,11 @@ public:
             const dictionary& dict
         );
 
-        //- Construct read
+        //- Construct read. Does findInstance to find io.local().
         distributedTriSurfaceMesh(const IOobject& io);
 
-        //- Construct from dictionary (used by searchableSurface)
+        //- Construct from dictionary (used by searchableSurface).
+        //  Does read. Does findInstance to find io.local().
         distributedTriSurfaceMesh
         (
             const IOobject& io,
diff --git a/src/meshTools/searchableSurface/searchableSphere.C b/src/meshTools/searchableSurface/searchableSphere.C
index ba5a1888cf..287ae6743c 100644
--- a/src/meshTools/searchableSurface/searchableSphere.C
+++ b/src/meshTools/searchableSurface/searchableSphere.C
@@ -53,7 +53,14 @@ Foam::pointIndexHit Foam::searchableSphere::findNearest
 
     if (nearestDistSqr > sqr(magN-radius_))
     {
-        info.rawPoint() = centre_ + n/magN*radius_;
+        if (magN < ROOTVSMALL)
+        {
+            info.rawPoint() = centre_ + vector(1,0,0)/magN*radius_;
+        }
+        else
+        {
+            info.rawPoint() = centre_ + n/magN*radius_;
+        }
         info.setHit();
         info.setIndex(0);
     }
diff --git a/src/meshTools/searchableSurface/searchableSurfaceCollection.C b/src/meshTools/searchableSurface/searchableSurfaceCollection.C
new file mode 100644
index 0000000000..ce6ddd6431
--- /dev/null
+++ b/src/meshTools/searchableSurface/searchableSurfaceCollection.C
@@ -0,0 +1,525 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "searchableSurfaceCollection.H"
+#include "addToRunTimeSelectionTable.H"
+#include "SortableList.H"
+#include "Time.H"
+#include "ListOps.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(searchableSurfaceCollection, 0);
+addToRunTimeSelectionTable(searchableSurface, searchableSurfaceCollection, dict);
+
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::searchableSurfaceCollection::findNearest
+(
+    const pointField& samples,
+    scalarField& minDistSqr,
+    List<pointIndexHit>& nearestInfo,
+    labelList& nearestSurf
+) const
+{
+    // Initialise
+    nearestInfo.setSize(samples.size());
+    nearestInfo = pointIndexHit();
+    nearestSurf.setSize(samples.size());
+    nearestSurf = -1;
+
+    List<pointIndexHit> hitInfo(samples.size());
+
+    const scalarField localMinDistSqr(samples.size(), GREAT);
+
+    forAll(subGeom_, surfI)
+    {
+        // Transform then divide
+        tmp<pointField> localSamples = cmptDivide
+        (
+            transform_[surfI].localPosition(samples),
+            scale_[surfI]
+        );
+
+        subGeom_[surfI].findNearest(localSamples, localMinDistSqr, hitInfo);
+
+        forAll(hitInfo, pointI)
+        {
+            if (hitInfo[pointI].hit())
+            {
+                // Rework back into global coordinate sys. Multiply then
+                // transform
+                point globalPt = transform_[surfI].globalPosition
+                (
+                    cmptMultiply
+                    (
+                        hitInfo[pointI].rawPoint(),
+                        scale_[surfI]
+                    )
+                );
+
+                scalar distSqr = magSqr(globalPt - samples[pointI]);
+
+                if (distSqr < minDistSqr[pointI])
+                {
+                    minDistSqr[pointI] = distSqr;
+                    nearestInfo[pointI].setPoint(globalPt);
+                    nearestInfo[pointI].setHit();
+                    nearestInfo[pointI].setIndex(hitInfo[pointI].index());
+                    nearestSurf[pointI] = surfI;
+                }
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceCollection::searchableSurfaceCollection
+(
+    const IOobject& io,
+    const dictionary& dict
+)
+:
+    searchableSurface(io),
+    instance_(dict.size()),
+    scale_(dict.size()),
+    transform_(dict.size()),
+    subGeom_(dict.size())
+{
+    Info<< "SearchableCollection : " << name() << endl;
+
+    label surfI = 0;
+    forAllConstIter(dictionary, dict, iter)
+    {
+        if (dict.isDict(iter().keyword()))
+        {
+            instance_[surfI] = iter().keyword();
+
+            const dictionary& subDict = dict.subDict(instance_[surfI]);
+
+            scale_[surfI] = subDict.lookup("scale");
+            transform_.set
+            (
+                surfI,
+                coordinateSystem::New
+                (
+                    "",
+                    subDict.subDict("transform")
+                )
+            );
+
+            const word subGeomName(subDict.lookup("surface"));
+            //Pout<< "Trying to find " << subGeomName << endl;
+
+            const searchableSurface& s =
+                io.db().lookupObject<searchableSurface>(subGeomName);
+
+            subGeom_.set(surfI, &const_cast<searchableSurface&>(s));
+
+            Info<< "    instance : " << instance_[surfI] << endl;
+            Info<< "    surface  : " << s.name() << endl;
+            Info<< "    scale    : " << scale_[surfI] << endl;
+            Info<< "    coordsys : " << transform_[surfI] << endl;
+
+            surfI++;
+        }
+    }
+    instance_.setSize(surfI);
+    scale_.setSize(surfI);
+    transform_.setSize(surfI);
+    subGeom_.setSize(surfI);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceCollection::~searchableSurfaceCollection()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::wordList& Foam::searchableSurfaceCollection::regions() const
+{
+    if (regions_.size() == 0)
+    {
+        regionOffset_.setSize(subGeom_.size());
+
+        DynamicList<word> allRegions;
+        forAll(subGeom_, surfI)
+        {
+            regionOffset_[surfI] = allRegions.size();
+
+            const wordList& subRegions = subGeom_[surfI].regions();
+
+            forAll(subRegions, i)
+            {
+                //allRegions.append(subRegions[i] + "_" + Foam::name(surfI));
+                allRegions.append(instance_[surfI] + "_" + subRegions[i]);
+            }
+        }
+        regions_.transfer(allRegions.shrink());
+    }
+    return regions_;
+}
+
+
+Foam::label Foam::searchableSurfaceCollection::size() const
+{
+    label n = 0;
+    forAll(subGeom_, surfI)
+    {
+        n += subGeom_[surfI].size();
+    }
+    return n;
+}
+
+
+void Foam::searchableSurfaceCollection::findNearest
+(
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+    List<pointIndexHit>& nearestInfo
+) const
+{
+    // How to scale distance?
+    scalarField minDistSqr(nearestDistSqr);
+
+    labelList nearestSurf;
+    findNearest
+    (
+        samples,
+        minDistSqr,
+        nearestInfo,
+        nearestSurf
+    );
+}
+
+
+void Foam::searchableSurfaceCollection::findLine
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    info.setSize(start.size());
+    info = pointIndexHit();
+
+    // Current nearest (to start) intersection
+    pointField nearest(end);
+
+    List<pointIndexHit> hitInfo(start.size());
+
+    forAll(subGeom_, surfI)
+    {
+        // Starting point
+        tmp<pointField> e0 = cmptDivide
+        (
+            transform_[surfI].localPosition
+            (
+                start
+            ),
+            scale_[surfI]
+        );
+
+        // Current best end point
+        tmp<pointField> e1 = cmptDivide
+        (
+            transform_[surfI].localPosition
+            (
+                nearest
+            ),
+            scale_[surfI]
+        );
+
+        subGeom_[surfI].findLine(e0, e1, hitInfo);
+
+        forAll(hitInfo, pointI)
+        {
+            if (hitInfo[pointI].hit())
+            {
+                // Transform back to global coordinate sys.
+                nearest[pointI] = transform_[surfI].globalPosition
+                (
+                    cmptMultiply
+                    (
+                        hitInfo[pointI].rawPoint(),
+                        scale_[surfI]
+                    )
+                );
+                info[pointI] = hitInfo[pointI];
+                info[pointI].rawPoint() = nearest[pointI];
+            }
+        }
+    }
+
+
+    // Debug check
+    if (false)
+    {
+        forAll(info, pointI)
+        {
+            if (info[pointI].hit())
+            {
+                vector n(end[pointI] - start[pointI]);
+                scalar magN = mag(n);
+
+                if (magN > SMALL)
+                {
+                    n /= mag(n);
+
+                    scalar s = ((info[pointI].rawPoint()-start[pointI])&n);
+
+                    if (s < 0 || s > 1)
+                    {
+                        FatalErrorIn
+                        (
+                            "searchableSurfaceCollection::findLine(..)"
+                        )   << "point:" << info[pointI]
+                            << " s:" << s
+                            << " outside vector "
+                            << " start:" << start[pointI]
+                            << " end:" << end[pointI]
+                            << abort(FatalError);
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::searchableSurfaceCollection::findLineAny
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    // To be done ...
+    findLine(start, end, info);
+}
+
+
+void Foam::searchableSurfaceCollection::findLineAll
+(
+    const pointField& start,
+    const pointField& end,
+    List<List<pointIndexHit> >& info
+) const
+{
+    // To be done. Assume for now only one intersection.
+    List<pointIndexHit> nearestInfo;
+    findLine(start, end, nearestInfo);
+
+    info.setSize(start.size());
+    forAll(info, pointI)
+    {
+        if (nearestInfo[pointI].hit())
+        {
+            info[pointI].setSize(1);
+            info[pointI][0] = nearestInfo[pointI];
+        }
+        else
+        {
+            info[pointI].clear();
+        }
+    }
+}
+
+
+void Foam::searchableSurfaceCollection::getRegion
+(
+    const List<pointIndexHit>& info,
+    labelList& region
+) const
+{
+    if (subGeom_.size() == 0)
+    {}
+    else if (subGeom_.size() == 1)
+    {
+        subGeom_[0].getRegion(info, region);
+    }
+    else
+    {
+        region.setSize(info.size());
+        region = -1;
+
+        // Which region did point come from. Retest for now to see which
+        // surface it originates from - crap solution! Should use global indices
+        // in index inside pointIndexHit to do this better.
+
+        pointField samples(info.size());
+        forAll(info, pointI)
+        {
+            if (info[pointI].hit())
+            {
+                samples[pointI] = info[pointI].hitPoint();
+            }
+            else
+            {
+                samples[pointI] = vector::zero;
+            }
+        }
+        //scalarField minDistSqr(info.size(), SMALL);
+        scalarField minDistSqr(info.size(), GREAT);
+
+        labelList nearestSurf;
+        List<pointIndexHit> nearestInfo;
+        findNearest
+        (
+            samples,
+            minDistSqr,
+            nearestInfo,
+            nearestSurf
+        );
+
+        // Check
+        {
+            forAll(info, pointI)
+            {
+                if (info[pointI].hit() && nearestSurf[pointI] == -1)
+                {
+                    FatalErrorIn
+                    (
+                        "searchableSurfaceCollection::getRegion(..)"
+                    )   << "pointI:" << pointI
+                        << " sample:" << samples[pointI]
+                        << " nearest:" << nearestInfo[pointI]
+                        << " nearestsurf:" << nearestSurf[pointI]
+                        << abort(FatalError);
+                }
+            }
+        }
+
+        forAll(subGeom_, surfI)
+        {
+            // Collect points from my surface
+            labelList indices(findIndices(nearestSurf, surfI));
+
+            labelList surfRegion;
+            subGeom_[surfI].getRegion
+            (
+                IndirectList<pointIndexHit>(info, indices),
+                surfRegion
+            );
+            forAll(indices, i)
+            {
+                region[indices[i]] = regionOffset_[surfI] + surfRegion[i];
+            }
+        }
+    }
+}
+
+
+void Foam::searchableSurfaceCollection::getNormal
+(
+    const List<pointIndexHit>& info,
+    vectorField& normal
+) const
+{
+    if (subGeom_.size() == 0)
+    {}
+    else if (subGeom_.size() == 1)
+    {
+        subGeom_[0].getNormal(info, normal);
+    }
+    else
+    {
+        normal.setSize(info.size());
+
+        // See above - crap retest to find surface point originates from.
+        pointField samples(info.size());
+        forAll(info, pointI)
+        {
+            if (info[pointI].hit())
+            {
+                samples[pointI] = info[pointI].hitPoint();
+            }
+            else
+            {
+                samples[pointI] = vector::zero;
+            }
+        }
+        //scalarField minDistSqr(info.size(), SMALL);
+        scalarField minDistSqr(info.size(), GREAT);
+
+        labelList nearestSurf;
+        List<pointIndexHit> nearestInfo;
+        findNearest
+        (
+            samples,
+            minDistSqr,
+            nearestInfo,
+            nearestSurf
+        );
+
+
+        forAll(subGeom_, surfI)
+        {
+            // Collect points from my surface
+            labelList indices(findIndices(nearestSurf, surfI));
+
+            vectorField surfNormal;
+            subGeom_[surfI].getNormal
+            (
+                IndirectList<pointIndexHit>(info, indices),
+                surfNormal
+            );
+            forAll(indices, i)
+            {
+                normal[indices[i]] = surfNormal[i];
+            }
+        }
+    }
+}
+
+
+void Foam::searchableSurfaceCollection::getVolumeType
+(
+    const pointField& points,
+    List<volumeType>& volType
+) const
+{
+    FatalErrorIn
+    (
+        "searchableSurfaceCollection::getVolumeType(const pointField&"
+        ", List<volumeType>&) const"
+    )   << "Volume type not supported for collection."
+        << exit(FatalError);
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/searchableSurfaceCollection.H b/src/meshTools/searchableSurface/searchableSurfaceCollection.H
new file mode 100644
index 0000000000..c85eaa1026
--- /dev/null
+++ b/src/meshTools/searchableSurface/searchableSurfaceCollection.H
@@ -0,0 +1,212 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::searchableSurfaceCollection
+
+Description
+    Union of transformed searchableSurfaces
+
+SourceFiles
+    searchableSurfaceCollection.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchableSurfaceCollection_H
+#define searchableSurfaceCollection_H
+
+#include "searchableSurface.H"
+#include "treeBoundBox.H"
+#include "coordinateSystem.H"
+#include "UPtrList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+/*---------------------------------------------------------------------------*\
+                           Class searchableSurfaceCollection Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchableSurfaceCollection
+:
+    public searchableSurface
+{
+private:
+
+    // Private Member Data
+
+        // Per instance data
+
+            //- instance name
+            wordList instance_;
+
+            //- scaling vector
+            vectorField scale_;
+
+            //- transformation
+            PtrList<coordinateSystem> transform_;
+
+            UPtrList<searchableSurface> subGeom_;
+
+        //- Region names
+        mutable wordList regions_;
+        //- From individual regions to collection regions
+        mutable labelList regionOffset_;
+
+
+    // Private Member Functions
+
+        //- Find point nearest to sample. Updates minDistSqr. Sets nearestInfo
+        //  and surface index
+        void findNearest
+        (
+            const pointField& samples,
+            scalarField& minDistSqr,
+            List<pointIndexHit>& nearestInfo,
+            labelList& nearestSurf
+        ) const;
+
+
+        //- Disallow default bitwise copy construct
+        searchableSurfaceCollection(const searchableSurfaceCollection&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const searchableSurfaceCollection&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("searchableSurfaceCollection");
+
+
+    // Constructors
+
+        //- Construct from dictionary (used by searchableSurface)
+        searchableSurfaceCollection
+        (
+            const IOobject& io,
+            const dictionary& dict
+        );
+
+    // Destructor
+
+        virtual ~searchableSurfaceCollection();
+
+
+    // Member Functions
+
+        virtual const wordList& regions() const;
+
+        //- Whether supports volume type below
+        virtual bool hasVolumeType() const
+        {
+            return false;
+        }
+
+        //- Range of local indices that can be returned.
+        virtual label size() const;
+
+
+        // Multiple point queries.
+
+            virtual void findNearest
+            (
+                const pointField& sample,
+                const scalarField& nearestDistSqr,
+                List<pointIndexHit>&
+            ) const;
+
+            virtual void findLine
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            virtual void findLineAny
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            //- Get all intersections in order from start to end.
+            virtual void findLineAll
+            (
+                const pointField& start,
+                const pointField& end,
+                List<List<pointIndexHit> >&
+            ) const;
+
+            //- From a set of points and indices get the region
+            virtual void getRegion
+            (
+                const List<pointIndexHit>&,
+                labelList& region
+            ) const;
+
+            //- From a set of points and indices get the normal
+            virtual void getNormal
+            (
+                const List<pointIndexHit>&,
+                vectorField& normal
+            ) const;
+
+            //- Determine type (inside/outside/mixed) for point. unknown if
+            //  cannot be determined (e.g. non-manifold surface)
+            virtual void getVolumeType
+            (
+                const pointField&,
+                List<volumeType>&
+            ) const;
+
+
+        // regIOobject implementation
+
+            bool writeData(Ostream&) const
+            {
+                notImplemented
+                (
+                    "searchableSurfaceCollection::writeData(Ostream&) const"
+                );
+                return false;
+            }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/searchableSurfaces.C b/src/meshTools/searchableSurface/searchableSurfaces.C
index 46c968f153..cec34bd7e7 100644
--- a/src/meshTools/searchableSurface/searchableSurfaces.C
+++ b/src/meshTools/searchableSurface/searchableSurfaces.C
@@ -39,8 +39,6 @@ defineTypeNameAndDebug(searchableSurfaces, 0);
 }
 
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct with length.
@@ -187,14 +185,6 @@ Foam::searchableSurfaces::searchableSurfaces
         // their object name. Maybe have stlTriSurfaceMesh which appends .stl
         // when reading/writing?
         namedIO().rename(key);  // names_[surfI]
-        if (namedIO().local() != word::null)
-        {
-            namedIO().instance() = namedIO().time().findInstance
-            (
-                namedIO().local(),
-                namedIO().name()
-            );
-        }
 
         // Create and hook surface
         set
diff --git a/src/meshTools/searchableSurface/searchableSurfaces.H b/src/meshTools/searchableSurface/searchableSurfaces.H
index ea91fcdee9..af02ae1883 100644
--- a/src/meshTools/searchableSurface/searchableSurfaces.H
+++ b/src/meshTools/searchableSurface/searchableSurfaces.H
@@ -69,6 +69,8 @@ class searchableSurfaces
         labelList allSurfaces_;
 
 
+    // Private Member Functions
+
         //- Disallow default bitwise copy construct
         searchableSurfaces(const searchableSurfaces&);
 
diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.C b/src/meshTools/searchableSurface/triSurfaceMesh.C
index 7b9abc1aec..9958afd1b1 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.C
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.C
@@ -43,6 +43,64 @@ addToRunTimeSelectionTable(searchableSurface, triSurfaceMesh, dict);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+//// Special version of Time::findInstance that does not check headerOk
+//// to search for instances of raw files
+//Foam::word Foam::triSurfaceMesh::findRawInstance
+//(
+//    const Time& runTime,
+//    const fileName& dir,
+//    const word& name
+//)
+//{
+//    // Check current time first
+//    if (isFile(runTime.path()/runTime.timeName()/dir/name))
+//    {
+//        return runTime.timeName();
+//    }
+//    instantList ts = runTime.times();
+//    label instanceI;
+//
+//    for (instanceI = ts.size()-1; instanceI >= 0; --instanceI)
+//    {
+//        if (ts[instanceI].value() <= runTime.timeOutputValue())
+//        {
+//            break;
+//        }
+//    }
+//
+//    // continue searching from here
+//    for (; instanceI >= 0; --instanceI)
+//    {
+//        if (isFile(runTime.path()/ts[instanceI].name()/dir/name))
+//        {
+//            return ts[instanceI].name();
+//        }
+//    }
+//
+//
+//    // not in any of the time directories, try constant
+//
+//    // Note. This needs to be a hard-coded constant, rather than the
+//    // constant function of the time, because the latter points to
+//    // the case constant directory in parallel cases
+//
+//    if (isFile(runTime.path()/runTime.constant()/dir/name))
+//    {
+//        return runTime.constant();
+//    }
+//
+//    FatalErrorIn
+//    (
+//        "searchableSurfaces::findRawInstance"
+//        "(const Time&, const fileName&, const word&)"
+//    )   << "Cannot find file \"" << name << "\" in directory "
+//        << runTime.constant()/dir
+//        << exit(FatalError);
+//
+//    return runTime.constant();
+//}
+
+
 //- Check file existence
 const Foam::fileName& Foam::triSurfaceMesh::checkFile
 (
@@ -149,7 +207,19 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const
 Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
 :
     searchableSurface(io),
-    objectRegistry(io),
+    objectRegistry
+    (
+        IOobject
+        (
+            io.name(),
+            io.instance(),
+            io.local(),
+            io.db(),
+            io.readOpt(),
+            io.writeOpt(),
+            false       // searchableSurface already registered under name
+        )
+    ),
     triSurface(s),
     surfaceClosed_(-1)
 {}
@@ -157,8 +227,34 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
 
 Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
 :
+    // Find instance for triSurfaceMesh
     searchableSurface(io),
-    objectRegistry(io),
+    //(
+    //    IOobject
+    //    (
+    //        io.name(),
+    //        io.time().findInstance(io.local(), word::null),
+    //        io.local(),
+    //        io.db(),
+    //        io.readOpt(),
+    //        io.writeOpt(),
+    //        io.registerObject()
+    //    )
+    //),
+    // Reused found instance in objectRegistry
+    objectRegistry
+    (
+        IOobject
+        (
+            io.name(),
+            static_cast<const searchableSurface&>(*this).instance(),
+            io.local(),
+            io.db(),
+            io.readOpt(),
+            io.writeOpt(),
+            false       // searchableSurface already registered under name
+        )
+    ),
     triSurface
     (
         checkFile
@@ -178,7 +274,32 @@ Foam::triSurfaceMesh::triSurfaceMesh
 )
 :
     searchableSurface(io),
-    objectRegistry(io),
+    //(
+    //    IOobject
+    //    (
+    //        io.name(),
+    //        io.time().findInstance(io.local(), word::null),
+    //        io.local(),
+    //        io.db(),
+    //        io.readOpt(),
+    //        io.writeOpt(),
+    //        io.registerObject()
+    //    )
+    //),
+    // Reused found instance in objectRegistry
+    objectRegistry
+    (
+        IOobject
+        (
+            io.name(),
+            static_cast<const searchableSurface&>(*this).instance(),
+            io.local(),
+            io.db(),
+            io.readOpt(),
+            io.writeOpt(),
+            false       // searchableSurface already registered under name
+        )
+    ),
     triSurface
     (
         checkFile
diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.H b/src/meshTools/searchableSurface/triSurfaceMesh.H
index aa0af5979c..a4231457d5 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.H
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.H
@@ -76,6 +76,14 @@ private:
 
     // Private Member Functions
 
+        ////- Helper: find instance of files without header
+        //static word findRawInstance
+        //(
+        //    const Time&,
+        //    const fileName&,
+        //    const word&
+        //);
+
         //- Check file existence
         static const fileName& checkFile
         (
@@ -105,10 +113,10 @@ public:
         //- Construct from triSurface
         triSurfaceMesh(const IOobject&, const triSurface&);
 
-        //- Construct read
+        //- Construct read.
         triSurfaceMesh(const IOobject& io);
 
-        //- Construct from dictionary (used by searchableSurface)
+        //- Construct from IO and dictionary (used by searchableSurface).
         //  Dictionary may contain a 'scale' entry (eg, 0.001: mm -> m)
         triSurfaceMesh
         (
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C
index c69a51f407..78081050ab 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C
@@ -30,6 +30,7 @@ License
 #include "syncTools.H"
 #include "addToRunTimeSelectionTable.H"
 #include "slicedVolFields.H"
+#include "surfaceFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -84,9 +85,74 @@ bool Foam::isoSurface::isEdgeOfFaceCut
 }
 
 
+// Get neighbour value and position.
+void Foam::isoSurface::getNeighbour
+(
+    const labelList& boundaryRegion,
+    const volScalarField& cVals,
+    const label cellI,
+    const label faceI,
+    scalar& nbrValue,
+    point& nbrPoint
+) const
+{
+    const labelList& own = mesh_.faceOwner();
+    const labelList& nei = mesh_.faceNeighbour();
+    const surfaceScalarField& weights = mesh_.weights();
+
+    if (mesh_.isInternalFace(faceI))
+    {
+        label nbr = (own[faceI] == cellI ? nei[faceI] : own[faceI]);
+        nbrValue = cVals[nbr];
+        nbrPoint = mesh_.cellCentres()[nbr];
+    }
+    else
+    {
+        label bFaceI = faceI-mesh_.nInternalFaces();
+        label patchI = boundaryRegion[bFaceI];
+        const polyPatch& pp = mesh_.boundaryMesh()[patchI];
+        label patchFaceI = faceI-pp.start();
+
+        if (isA<emptyPolyPatch>(pp))
+        {
+            // Assume zero gradient
+            nbrValue = cVals[own[faceI]];
+            nbrPoint = mesh_.faceCentres()[faceI];
+        }
+        else if (pp.coupled())
+        {
+            if (!refCast<const coupledPolyPatch>(pp).separated())
+            {
+                // Behave as internal face:
+                // other side value
+                nbrValue = cVals.boundaryField()[patchI][patchFaceI];
+                // other side cell centre
+                nbrPoint = mesh_.C().boundaryField()[patchI][patchFaceI];
+            }
+            else
+            {
+                // Do some interpolation for now
+                const scalarField& w = weights.boundaryField()[patchI];
+
+                nbrPoint = mesh_.faceCentres()[faceI];
+                nbrValue =
+                    (1-w[patchFaceI])*cVals[own[faceI]]
+                  + w[patchFaceI]*cVals.boundaryField()[patchI][patchFaceI];
+            }
+        }
+        else
+        {
+            nbrValue = cVals.boundaryField()[patchI][patchFaceI];
+            nbrPoint = mesh_.faceCentres()[faceI];
+        }
+    }
+}
+
+
 // Determine for every face/cell whether it (possibly) generates triangles.
 void Foam::isoSurface::calcCutTypes
 (
+    const labelList& boundaryRegion,
     const volScalarField& cVals,
     const scalarField& pVals
 )
@@ -102,7 +168,20 @@ void Foam::isoSurface::calcCutTypes
     {
         // CC edge.
         bool ownLower = (cVals[own[faceI]] < iso_);
-        bool neiLower = (cVals[nei[faceI]] < iso_);
+
+        scalar nbrValue;
+        point nbrPoint;
+        getNeighbour
+        (
+            boundaryRegion,
+            cVals,
+            own[faceI],
+            faceI,
+            nbrValue,
+            nbrPoint
+        );
+
+        bool neiLower = (nbrValue < iso_);
 
         if (ownLower != neiLower)
         {
@@ -120,53 +199,47 @@ void Foam::isoSurface::calcCutTypes
             }
         }
     }
+
     forAll(patches, patchI)
     {
         const polyPatch& pp = patches[patchI];
 
         label faceI = pp.start();
 
-        if (isA<emptyPolyPatch>(pp))
+        forAll(pp, i)
         {
-            // Assume zero gradient so owner and neighbour/boundary value equal
+            bool ownLower = (cVals[own[faceI]] < iso_);
 
-            forAll(pp, i)
+            scalar nbrValue;
+            point nbrPoint;
+            getNeighbour
+            (
+                boundaryRegion,
+                cVals,
+                own[faceI],
+                faceI,
+                nbrValue,
+                nbrPoint
+            );
+
+            bool neiLower = (nbrValue < iso_);
+
+            if (ownLower != neiLower)
             {
-                bool ownLower = (cVals[own[faceI]] < iso_);
-
+                faceCutType_[faceI] = CUT;
+            }
+            else
+            {
+                // Mesh edge.
                 const face f = mesh_.faces()[faceI];
 
-                if (isEdgeOfFaceCut(pVals, f, ownLower, ownLower))
+                if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
                 {
                     faceCutType_[faceI] = CUT;
                 }
-
-                faceI++;
             }
-        }
-        else
-        {
-            forAll(pp, i)
-            {
-                bool ownLower = (cVals[own[faceI]] < iso_);
-                bool neiLower = (cVals.boundaryField()[patchI][i] < iso_);
 
-                if (ownLower != neiLower)
-                {
-                    faceCutType_[faceI] = CUT;
-                }
-                else
-                {
-                    // Mesh edge.
-                    const face f = mesh_.faces()[faceI];
-
-                    if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
-                    {
-                        faceCutType_[faceI] = CUT;
-                    }
-                }
-                faceI++;
-            }
+            faceI++;
         }
     }
 
@@ -330,59 +403,6 @@ Foam::pointIndexHit Foam::isoSurface::collapseSurface
 }
 
 
-// Get neighbour value and position.
-void Foam::isoSurface::getNeighbour
-(
-    const labelList& boundaryRegion,
-    const volScalarField& cVals,
-    const label cellI,
-    const label faceI,
-    scalar& nbrValue,
-    point& nbrPoint
-) const
-{
-    const labelList& own = mesh_.faceOwner();
-    const labelList& nei = mesh_.faceNeighbour();
-
-    if (mesh_.isInternalFace(faceI))
-    {
-        label nbr = (own[faceI] == cellI ? nei[faceI] : own[faceI]);
-        nbrValue = cVals[nbr];
-        nbrPoint = mesh_.cellCentres()[nbr];
-    }
-    else
-    {
-        label bFaceI = faceI-mesh_.nInternalFaces();
-        label patchI = boundaryRegion[bFaceI];
-        const polyPatch& pp = mesh_.boundaryMesh()[patchI];
-        label patchFaceI = faceI-pp.start();
-
-        if (isA<emptyPolyPatch>(pp))
-        {
-            // Assume zero gradient
-            nbrValue = cVals[own[faceI]];
-            nbrPoint = mesh_.faceCentres()[faceI];
-        }
-        else if
-        (
-            pp.coupled()
-        && !refCast<const coupledPolyPatch>(pp).separated()
-        )
-        {
-            // other side value
-            nbrValue = cVals.boundaryField()[patchI][patchFaceI];
-            // other side cell centre
-            nbrPoint = mesh_.C().boundaryField()[patchI][patchFaceI];
-        }
-        else
-        {
-            nbrValue = cVals.boundaryField()[patchI][patchFaceI];
-            nbrPoint = mesh_.faceCentres()[faceI];
-        }
-    }
-}
-
-
 // Determine per cell centre whether all the intersections get collapsed
 // to a single point
 void Foam::isoSurface::calcSnappedCc
@@ -606,13 +626,14 @@ void Foam::isoSurface::calcSnappedPoint
 
         forAll(pFaces, pFaceI)
         {
+            // Create points for all intersections close to point
+            // (i.e. from pyramid edges)
+
             label faceI = pFaces[pFaceI];
             const face& f = mesh_.faces()[faceI];
             label own = mesh_.faceOwner()[faceI];
 
-            // Create points for all intersections close to point
-            // (i.e. from pyramid edges)
-
+            // Get neighbour value
             scalar nbrValue;
             point nbrPoint;
             getNeighbour
@@ -831,6 +852,7 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
     }
 
 
+
     // Determine 'flat hole' situation (see RMT paper).
     // Two unconnected triangles get connected because (some of) the edges
     // separating them get collapsed. Below only checks for duplicate triangles,
@@ -846,22 +868,29 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
         forAll(tris, triI)
         {
             const labelledTri& tri = tris[triI];
-
             const labelList& pFaces = pointFaces[tri[0]];
 
-            // Find the minimum of any duplicates
+            // Find the maximum of any duplicates. Maximum since the tris
+            // below triI
+            // get overwritten so we cannot use them in a comparison.
             label dupTriI = -1;
             forAll(pFaces, i)
             {
-                if (pFaces[i] < triI && tris[pFaces[i]] == tri)
+                label nbrTriI = pFaces[i];
+
+                if (nbrTriI > triI && (tris[nbrTriI] == tri))
                 {
-                    dupTriI = pFaces[i];
+                    //Pout<< "Duplicate : " << triI << " verts:" << tri
+                    //    << " to " << nbrTriI << " verts:" << tris[nbrTriI]
+                    //    << endl;
+                    dupTriI = nbrTriI;
+                    break;
                 }
             }
 
             if (dupTriI == -1)
             {
-                // There is no lower triangle
+                // There is no (higher numbered) duplicate triangle
                 label newTriI = newToOldTri.size();
                 newToOldTri.append(triI);
                 tris[newTriI] = tris[triI];
@@ -876,6 +905,43 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
             Pout<< "isoSurface : merged from " << nTris
                 << " down to " << tris.size() << " unique triangles." << endl;
         }
+
+
+        if (debug)
+        {
+            triSurface surf(tris, geometricSurfacePatchList(0), newPoints);
+
+            forAll(surf, faceI)
+            {
+                const labelledTri& f = surf[faceI];
+                const labelList& fFaces = surf.faceFaces()[faceI];
+
+                forAll(fFaces, i)
+                {
+                    label nbrFaceI = fFaces[i];
+
+                    if (nbrFaceI <= faceI)
+                    {
+                        // lower numbered faces already checked
+                        continue;
+                    }
+
+                    const labelledTri& nbrF = surf[nbrFaceI];
+
+                    if (f == nbrF)
+                    {
+                        FatalErrorIn("validTri(const triSurface&, const label)")
+                            << "Check : "
+                            << " triangle " << faceI << " vertices " << f
+                            << " fc:" << f.centre(surf.points())
+                            << " has the same vertices as triangle " << nbrFaceI
+                            << " vertices " << nbrF
+                            << " fc:" << nbrF.centre(surf.points())
+                            << abort(FatalError);
+                    }
+                }
+            }
+        }
     }
 
     return triSurface(tris, geometricSurfacePatchList(0), newPoints, true);
@@ -1473,45 +1539,23 @@ Foam::isoSurface::isoSurface
 {
     if (debug)
     {
-        Pout<< "isoSurface     :" << nl
-            << "    isoField   : " << cVals.name() << nl
-            << "    isoValue   : " << iso << nl
-            << "    regularise : " << regularise << nl
-            << "    mergeTol   : " << mergeTol << nl
+        Pout<< "isoSurface:" << nl
+            << "    isoField      : " << cVals.name() << nl
+            << "    cell min/max  : "
+            << min(cVals.internalField()) << " / "
+            << max(cVals.internalField()) << nl
+            << "    point min/max : "
+            << min(pVals) << " / "
+            << max(pVals) << nl
+            << "    isoValue      : " << iso << nl
+            << "    regularise    : " << regularise << nl
+            << "    mergeTol      : " << mergeTol << nl
             << endl;
     }
 
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    // Check
-    forAll(patches, patchI)
-    {
-        if (isA<emptyPolyPatch>(patches[patchI]))
-        {
-            FatalErrorIn
-            (
-                "isoSurface::isoSurface\n"
-                "(\n"
-                "    const volScalarField& cVals,\n"
-                "    const scalarField& pVals,\n"
-                "    const scalar iso,\n"
-                "    const bool regularise,\n"
-                "    const scalar mergeTol\n"
-                ")\n"
-            )   << "Iso surfaces not supported on case with empty patches."
-                << exit(FatalError);
-        }
-    }
-
-
-    // Determine if any cut through face/cell
-    calcCutTypes(cVals, pVals);
-
-
-    // Determine if point is on boundary. Points on boundaries are never
-    // snapped. Coupled boundaries are handled explicitly so not marked here.
-    PackedBoolList isBoundaryPoint(mesh_.nPoints());
-
+    // Pre-calculate patch-per-face to avoid whichPatch call.
     labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces());
 
     forAll(patches, patchI)
@@ -1525,33 +1569,13 @@ Foam::isoSurface::isoSurface
             boundaryRegion[faceI-mesh_.nInternalFaces()] = patchI;
             faceI++;
         }
-
-        // Mark all points that are not physically coupled (so anything
-        // but collocated coupled patches)
-        if
-        (
-           !pp.coupled()
-        || refCast<const coupledPolyPatch>(pp).separated()
-        )
-        {
-            label faceI = pp.start();
-
-            forAll(pp, i)
-            {
-                boundaryRegion[faceI-mesh_.nInternalFaces()] = patchI;
-
-                const face& f = mesh_.faces()[faceI];
-
-                forAll(f, fp)
-                {
-                    isBoundaryPoint.set(f[fp], 1);
-                }
-                faceI++;
-            }
-        }
     }
 
 
+    // Determine if any cut through face/cell
+    calcCutTypes(boundaryRegion, cVals, pVals);
+
+
     DynamicList<point> snappedPoints(nCutCells_);
 
     // Per cc -1 or a point inside snappedPoints.
@@ -1584,6 +1608,39 @@ Foam::isoSurface::isoSurface
 
     label nCellSnaps = snappedPoints.size();
 
+
+    // Determine if point is on boundary. Points on boundaries are never
+    // snapped. Coupled boundaries are handled explicitly so not marked here.
+    PackedBoolList isBoundaryPoint(mesh_.nPoints());
+
+    forAll(patches, patchI)
+    {
+        // Mark all boundary points that are not physically coupled (so anything
+        // but collocated coupled patches)
+        const polyPatch& pp = patches[patchI];
+
+        if
+        (
+           !pp.coupled()
+        || refCast<const coupledPolyPatch>(pp).separated()
+        )
+        {
+            label faceI = pp.start();
+
+            forAll(pp, i)
+            {
+                const face& f = mesh_.faces()[faceI];
+
+                forAll(f, fp)
+                {
+                    isBoundaryPoint.set(f[fp], 1);
+                }
+                faceI++;
+            }
+        }
+    }
+
+
     // Per point -1 or a point inside snappedPoints.
     labelList snappedPoint;
     if (regularise)
@@ -1613,7 +1670,7 @@ Foam::isoSurface::isoSurface
 
 
     // Generate field to interpolate. This is identical to the mesh.C()
-    // except on separated coupled patches.
+    // except on separated coupled patches and on empty patches.
     slicedVolVectorField meshC
     (
         IOobject
@@ -1635,10 +1692,12 @@ Foam::isoSurface::isoSurface
         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
         forAll(patches, patchI)
         {
+            const polyPatch& pp = patches[patchI];
+
             if
             (
-                patches[patchI].coupled()
-             && refCast<const coupledPolyPatch>(patches[patchI]).separated()
+                pp.coupled()
+             && refCast<const coupledPolyPatch>(pp).separated()
             )
             {
                 fvPatchVectorField& pfld = const_cast<fvPatchVectorField&>
@@ -1647,9 +1706,32 @@ Foam::isoSurface::isoSurface
                 );
                 pfld.operator==
                 (
-                    patches[patchI].patchSlice(mesh_.faceCentres())
+                    pp.patchSlice(mesh_.faceCentres())
                 );
             }
+            else if (isA<emptyPolyPatch>(pp))
+            {
+                typedef slicedVolVectorField::GeometricBoundaryField bType;
+
+                bType& bfld = const_cast<bType&>(meshC.boundaryField());
+
+                // Clear old value. Cannot resize it since slice.
+                bfld.set(patchI, NULL);
+
+                // Set new value we can change
+                bfld.set
+                (
+                    patchI,
+                    new calculatedFvPatchField<vector>
+                    (
+                        mesh_.boundary()[patchI],
+                        meshC
+                    )
+                );
+
+                // Change to face centres
+                bfld[patchI] = pp.patchSlice(mesh_.faceCentres());
+            }
         }
     }
 
@@ -1676,7 +1758,8 @@ Foam::isoSurface::isoSurface
     if (debug)
     {
         Pout<< "isoSurface : generated " << triMeshCells.size()
-            << " unmerged triangles." << endl;
+            << " unmerged triangles from " << triPoints.size()
+            << " unmerged points." << endl;
     }
 
 
@@ -1773,6 +1856,14 @@ Foam::isoSurface::isoSurface
 
     orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
     //}
+
+
+    if (debug)
+    {
+        fileName stlFile = mesh_.time().path() + ".stl";
+        Pout<< "Dumping surface to " << stlFile << endl;
+        triSurface::write(stlFile);
+    }
 }
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.H b/src/sampling/sampledSurface/isoSurface/isoSurface.H
index 7dc0588a1c..b859524b2c 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.H
@@ -31,9 +31,6 @@ Description
     G.M. Treece, R.W. Prager and A.H. Gee.
 
     Note:
-    - not possible on patches of type 'empty'. There are no values on
-    'empty' patch fields so even the api would have to change
-     (no volScalarField as argument). Too messy.
     - in parallel the regularisation (coarsening) always takes place
     and slightly different surfaces will be created compared to non-parallel.
     The surface will still be continuous though!
@@ -43,6 +40,7 @@ Description
     - uses geometric merge with fraction of bounding box as distance.
     - triangles can be between two cell centres so constant sampling
       does not make sense.
+    - on empty patches behaves like zero gradient.
     - does not do 2D correctly, creates non-flat iso surface.
 
 SourceFiles
@@ -135,9 +133,20 @@ class isoSurface
             const bool neiLower
         ) const;
 
+        void getNeighbour
+        (
+            const labelList& boundaryRegion,
+            const volScalarField& cVals,
+            const label cellI,
+            const label faceI,
+            scalar& nbrValue,
+            point& nbrPoint
+        ) const;
+
         //- Set faceCutType,cellCutType.
         void calcCutTypes
         (
+            const labelList& boundaryRegion,
             const volScalarField& cVals,
             const scalarField& pVals
         );
@@ -156,16 +165,6 @@ class isoSurface
             DynamicList<labelledTri, 64>& localTris
         );
 
-        void getNeighbour
-        (
-            const labelList& boundaryRegion,
-            const volScalarField& cVals,
-            const label cellI,
-            const label faceI,
-            scalar& nbrValue,
-            point& nbrPoint
-        ) const;
-
         //- Determine per cc whether all near cuts can be snapped to single
         //  point.
         void calcSnappedCc
@@ -193,37 +192,39 @@ class isoSurface
         template<class Type>
         Type generatePoint
         (
-            const DynamicList<Type>& snappedPoints,
-
             const scalar s0,
             const Type& p0,
-            const label p0Index,
+            const bool hasSnap0,
+            const Type& snapP0,
 
             const scalar s1,
             const Type& p1,
-            const label p1Index
+            const bool hasSnap1,
+            const Type& snapP1
         ) const;
 
         template<class Type>
         void generateTriPoints
         (
-            const DynamicList<Type>& snapped,
-
             const scalar s0,
             const Type& p0,
-            const label p0Index,
+            const bool hasSnap0,
+            const Type& snapP0,
 
             const scalar s1,
             const Type& p1,
-            const label p1Index,
+            const bool hasSnap1,
+            const Type& snapP1,
 
             const scalar s2,
             const Type& p2,
-            const label p2Index,
+            const bool hasSnap2,
+            const Type& snapP2,
 
             const scalar s3,
             const Type& p3,
-            const label p3Index,
+            const bool hasSnap3,
+            const Type& snapP3,
 
             DynamicList<Type>& points
         ) const;
@@ -244,7 +245,8 @@ class isoSurface
 
             const scalar neiVal,
             const Type& neiPt,
-            const label neiSnap,
+            const bool hasNeiSnap,
+            const Type& neiSnapPt,
 
             DynamicList<Type>& triPoints,
             DynamicList<label>& triMeshCells
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
index e2251aad2a..0760b06a52 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
@@ -33,15 +33,15 @@ License
 template<class Type>
 Type Foam::isoSurface::generatePoint
 (
-    const DynamicList<Type>& snappedPoints,
-
     const scalar s0,
     const Type& p0,
-    const label p0Index,
+    const bool hasSnap0,
+    const Type& snapP0,
 
     const scalar s1,
     const Type& p1,
-    const label p1Index
+    const bool hasSnap1,
+    const Type& snapP1
 ) const
 {
     scalar d = s1-s0;
@@ -50,13 +50,13 @@ Type Foam::isoSurface::generatePoint
     {
         scalar s = (iso_-s0)/d;
 
-        if (s >= 0.5 && s <= 1 && p1Index != -1)
+        if (hasSnap1 && s >= 0.5 && s <= 1)
         {
-            return snappedPoints[p1Index];
+            return snapP1;
         }
-        else if (s >= 0.0 && s <= 0.5 && p0Index != -1)
+        else if (hasSnap0 && s >= 0.0 && s <= 0.5)
         {
-            return snappedPoints[p0Index];
+            return snapP0;
         }
         else
         {
@@ -75,23 +75,25 @@ Type Foam::isoSurface::generatePoint
 template<class Type>
 void Foam::isoSurface::generateTriPoints
 (
-    const DynamicList<Type>& snapped,
-
     const scalar s0,
     const Type& p0,
-    const label p0Index,
+    const bool hasSnap0,
+    const Type& snapP0,
 
     const scalar s1,
     const Type& p1,
-    const label p1Index,
+    const bool hasSnap1,
+    const Type& snapP1,
 
     const scalar s2,
     const Type& p2,
-    const label p2Index,
+    const bool hasSnap2,
+    const Type& snapP2,
 
     const scalar s3,
     const Type& p3,
-    const label p3Index,
+    const bool hasSnap3,
+    const Type& snapP3,
 
     DynamicList<Type>& points
 ) const
@@ -123,29 +125,55 @@ void Foam::isoSurface::generateTriPoints
 
         case 0x0E:
         case 0x01:
-            points.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
-            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            points.append
+            (
+                generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1)
+            );
+            points.append
+            (
+                generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
+            );
+            points.append
+            (
+                generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
+            );
         break;
 
         case 0x0D:
         case 0x02:
-            points.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            points.append
+            (
+                generatePoint(s1,p1,hasSnap1,snapP1,s0,p0,hasSnap0,snapP0)
+            );
+            points.append
+            (
+                generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
+            );
+            points.append
+            (
+                generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
+            );
         break;
 
         case 0x0C:
         case 0x03:
         {
-            Type tp1 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
-            Type tp2 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
+            Type tp1 =
+                generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2);
+            Type tp2 =
+                generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3);
 
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            points.append
+            (
+                generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
+            );
             points.append(tp1);
             points.append(tp2);
             points.append(tp2);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            points.append
+            (
+                generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
+            );
             points.append(tp1);
         }
         break;
@@ -153,23 +181,40 @@ void Foam::isoSurface::generateTriPoints
         case 0x0B:
         case 0x04:
         {
-            points.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
-            points.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
+            points.append
+            (
+                generatePoint(s2,p2,hasSnap2,snapP2,s0,p0,hasSnap0,snapP0)
+            );
+            points.append
+            (
+                generatePoint(s2,p2,hasSnap2,snapP2,s1,p1,hasSnap1,snapP1)
+            );
+            points.append
+            (
+                generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3)
+            );
         }
         break;
 
         case 0x0A:
         case 0x05:
         {
-            Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+            Type tp0 =
+                generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
+            Type tp1 =
+                generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
 
             points.append(tp0);
             points.append(tp1);
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            points.append
+            (
+                generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
+            );
             points.append(tp0);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            points.append
+            (
+                generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
+            );
             points.append(tp1);
         }
         break;
@@ -177,23 +222,40 @@ void Foam::isoSurface::generateTriPoints
         case 0x09:
         case 0x06:
         {
-            Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+            Type tp0 =
+                generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
+            Type tp1 =
+                generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
 
             points.append(tp0);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+            points.append
+            (
+                generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
+            );
             points.append(tp1);
             points.append(tp0);
-            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+            points.append
+            (
+                generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
+            );
             points.append(tp1);
         }
         break;
 
         case 0x07:
         case 0x08:
-            points.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
-            points.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
+            points.append
+            (
+                generatePoint(s3,p3,hasSnap3,snapP3,s0,p0,hasSnap0,snapP0)
+            );
+            points.append
+            (
+                generatePoint(s3,p3,hasSnap3,snapP3,s2,p2,hasSnap2,snapP2)
+            );
+            points.append
+            (
+                generatePoint(s3,p3,hasSnap3,snapP3,s1,p1,hasSnap1,snapP1)
+            );
         break;
     }
 }
@@ -215,7 +277,8 @@ Foam::label Foam::isoSurface::generateFaceTriPoints
 
     const scalar neiVal,
     const Type& neiPt,
-    const label neiSnap,
+    const bool hasNeiSnap,
+    const Type& neiSnapPt,
 
     DynamicList<Type>& triPoints,
     DynamicList<label>& triMeshCells
@@ -234,23 +297,37 @@ Foam::label Foam::isoSurface::generateFaceTriPoints
 
         generateTriPoints
         (
-            snappedPoints,
-
             pVals[pointI],
             pCoords[pointI],
-            snappedPoint[pointI],
+            snappedPoint[pointI] != -1,
+            (
+                snappedPoint[pointI] != -1
+              ? snappedPoints[snappedPoint[pointI]]
+              : pTraits<Type>::zero
+            ),
 
             pVals[nextPointI],
             pCoords[nextPointI],
-            snappedPoint[nextPointI],
+            snappedPoint[nextPointI] != -1,
+            (
+                snappedPoint[nextPointI] != -1
+              ? snappedPoints[snappedPoint[nextPointI]]
+              : pTraits<Type>::zero
+            ),
 
             cVals[own],
             cCoords[own],
-            snappedCc[own],
+            snappedCc[own] != -1,
+            (
+                snappedCc[own] != -1
+              ? snappedPoints[snappedCc[own]]
+              : pTraits<Type>::zero
+            ),
 
             neiVal,
             neiPt,
-            neiSnap,
+            hasNeiSnap,
+            neiSnapPt,
 
             triPoints
         );
@@ -311,25 +388,6 @@ void Foam::isoSurface::generateTriPoints
             << abort(FatalError);
     }
 
-    // Determine neighbouring snap status
-    labelList neiSnappedCc(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (pp.coupled())
-        {
-            label faceI = pp.start();
-            forAll(pp, i)
-            {
-                neiSnappedCc[faceI-mesh_.nInternalFaces()] =
-                    snappedCc[own[faceI]];
-                faceI++;
-            }
-        }
-    }
-    syncTools::swapBoundaryFaceList(mesh_, neiSnappedCc, false);
-
 
 
     // Generate triangle points
@@ -356,7 +414,12 @@ void Foam::isoSurface::generateTriPoints
 
                 cVals[nei[faceI]],
                 cCoords[nei[faceI]],
-                snappedCc[nei[faceI]],
+                snappedCc[nei[faceI]] != -1,
+                (
+                    snappedCc[nei[faceI]] != -1
+                  ? snappedPoints[snappedCc[nei[faceI]]]
+                  : pTraits<Type>::zero
+                ),
 
                 triPoints,
                 triMeshCells
@@ -365,6 +428,34 @@ void Foam::isoSurface::generateTriPoints
     }
 
 
+    // Determine neighbouring snap status
+    boolList neiSnapped(mesh_.nFaces()-mesh_.nInternalFaces(), false);
+    List<Type> neiSnappedPoint(neiSnapped.size(), pTraits<Type>::zero);
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+            forAll(pp, i)
+            {
+                label bFaceI = faceI-mesh_.nInternalFaces();
+                label snappedIndex = snappedCc[own[faceI]];
+
+                if (snappedIndex != -1)
+                {
+                    neiSnapped[bFaceI] = true;
+                    neiSnappedPoint[bFaceI] = snappedPoints[snappedIndex];
+                }
+                faceI++;
+            }
+        }
+    }
+    syncTools::swapBoundaryFaceList(mesh_, neiSnapped, false);
+    syncTools::swapBoundaryFaceList(mesh_, neiSnappedPoint, false);
+
+
     forAll(patches, patchI)
     {
         const polyPatch& pp = patches[patchI];
@@ -372,45 +463,104 @@ void Foam::isoSurface::generateTriPoints
         if
         (
             isA<processorPolyPatch>(pp)
-         && refCast<const processorPolyPatch>(pp).owner()
          && !refCast<const processorPolyPatch>(pp).separated()
         )
         {
-            label faceI = pp.start();
-
-            forAll(pp, i)
+            //if (refCast<const processorPolyPatch>(pp).owner())
             {
-                if (faceCutType_[faceI] != NOTCUT)
+                label faceI = pp.start();
+
+                forAll(pp, i)
                 {
-                    generateFaceTriPoints
-                    (
-                        cVals,
-                        pVals,
+                    if (faceCutType_[faceI] != NOTCUT)
+                    {
+                        label bFaceI = faceI-mesh_.nInternalFaces();
+                        if
+                        (
+                            neiSnapped[bFaceI]
+                         && (neiSnappedPoint[bFaceI]==pTraits<Type>::zero)
+                        )
+                        {
+                            FatalErrorIn("isoSurface::generateTriPoints(..)")
+                                << "problem:" << abort(FatalError);
+                        }
 
-                        cCoords,
-                        pCoords,
 
-                        snappedPoints,
-                        snappedCc,
-                        snappedPoint,
-                        faceI,
+                        generateFaceTriPoints
+                        (
+                            cVals,
+                            pVals,
 
-                        cVals.boundaryField()[patchI][i],
-                        cCoords.boundaryField()[patchI][i],
-                        neiSnappedCc[faceI-mesh_.nInternalFaces()],
+                            cCoords,
+                            pCoords,
 
-                        triPoints,
-                        triMeshCells
-                    );
+                            snappedPoints,
+                            snappedCc,
+                            snappedPoint,
+                            faceI,
+
+                            cVals.boundaryField()[patchI][i],
+                            cCoords.boundaryField()[patchI][i],
+                            neiSnapped[faceI-mesh_.nInternalFaces()],
+                            neiSnappedPoint[faceI-mesh_.nInternalFaces()],
+
+                            triPoints,
+                            triMeshCells
+                        );
+                    }
+                    faceI++;
                 }
-                faceI++;
             }
         }
         else if (isA<emptyPolyPatch>(pp))
         {
-            // Assume zero-gradient. But what about coordinates?
+            // Check if field is there (when generating geometry the
+            // empty patches have been rewritten to be the face centres),
+            // otherwise use zero-gradient.
+
             label faceI = pp.start();
 
+
+            const fvPatchScalarField& fvp = cVals.boundaryField()[patchI];
+
+            // Owner value of cVals
+            scalarField internalVals;
+            if (fvp.size() == 0)
+            {
+                internalVals.setSize(pp.size());
+                forAll(pp, i)
+                {
+                    internalVals[i] = cVals[own[pp.start()+i]];
+                }
+            }
+            const scalarField& bVals =
+            (
+                fvp.size() > 0
+              ? static_cast<const scalarField&>(fvp)
+              : internalVals
+            );
+
+
+            const fvPatchField<Type>& pc = cCoords.boundaryField()[patchI];
+
+            // Owner value of cCoords
+            Field<Type> internalCoords;
+            if (pc.size() == 0)
+            {
+                internalCoords.setSize(pp.size());
+                forAll(pp, i)
+                {
+                    internalCoords[i] = cCoords[own[pp.start()+i]];
+                }
+            }
+            const Field<Type>& bCoords =
+            (
+                pc.size() > 0
+              ? static_cast<const Field<Type>&>(pc)
+              : internalCoords
+            );
+
+
             forAll(pp, i)
             {
                 if (faceCutType_[faceI] != NOTCUT)
@@ -428,9 +578,10 @@ void Foam::isoSurface::generateTriPoints
                         snappedPoint,
                         faceI,
 
-                        cVals[own[faceI]],
-                        cCoords.boundaryField()[patchI][i],
-                        -1, // fc not snapped
+                        bVals[i],
+                        bCoords[i],
+                        false,  // fc not snapped
+                        pTraits<Type>::zero,
 
                         triPoints,
                         triMeshCells
@@ -462,7 +613,8 @@ void Foam::isoSurface::generateTriPoints
 
                         cVals.boundaryField()[patchI][i],
                         cCoords.boundaryField()[patchI][i],
-                        -1, // fc not snapped
+                        false,  // fc not snapped
+                        pTraits<Type>::zero,
 
                         triPoints,
                         triMeshCells
diff --git a/src/turbulenceModels/LES/Allwmake b/src/turbulenceModels/LES/Allwmake
index a80c71ab57..d83439a913 100755
--- a/src/turbulenceModels/LES/Allwmake
+++ b/src/turbulenceModels/LES/Allwmake
@@ -2,7 +2,7 @@
 cd ${0%/*} || exit 1    # run from this directory
 set -x
 
-wmakeLnInclude -f ../incompressible/LES
+wmakeLnInclude -f ../incompressible/LES -sf
 
 wmake libso LESfilters
 wmake libso LESdeltas
diff --git a/src/turbulenceModels/LES/LESdeltas/cubeRootVolDelta/cubeRootVolDelta.C b/src/turbulenceModels/LES/LESdeltas/cubeRootVolDelta/cubeRootVolDelta.C
index a3ce3ab3a9..2002deaa21 100644
--- a/src/turbulenceModels/LES/LESdeltas/cubeRootVolDelta/cubeRootVolDelta.C
+++ b/src/turbulenceModels/LES/LESdeltas/cubeRootVolDelta/cubeRootVolDelta.C
@@ -42,8 +42,7 @@ addToRunTimeSelectionTable(LESdelta, cubeRootVolDelta, dictionary);
 
 void cubeRootVolDelta::calcDelta()
 {
-    const Vector<label>& directions = mesh().directions();
-    label nD = (directions.nComponents + cmptSum(directions))/2;
+    label nD = mesh().nGeometricD();
 
     if (nD == 3)
     {
@@ -55,14 +54,15 @@ void cubeRootVolDelta::calcDelta()
             << "Case is 2D, LES is not strictly applicable\n"
             << endl;
 
+        const Vector<label>& directions = mesh().geometricD();
+
         scalar thickness = 0.0;
         for (direction dir=0; dir<directions.nComponents; dir++)
         {
             if (directions[dir] == -1)
             {
-                boundBox bb(mesh().points(), false);
-
-                thickness = bb.span()[dir];
+                thickness = mesh().bounds().span()[dir];
+                break;
             }
         }
 
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.C b/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.C
index 3f3aa7699b..005ee2eb16 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.C
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.C
@@ -41,8 +41,7 @@ namespace Foam
 
 void Foam::IDDESDelta::calcDelta()
 {
-    const Vector<label>& directions = mesh().directions();
-    label nD = (directions.nComponents + cmptSum(directions))/2;
+    label nD = mesh().nGeometricD();
 
     // initialise hwn as wall distance
     volScalarField hwn = wallDist(mesh()).y();
diff --git a/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C b/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
index 9fc36826fe..b3d7312a06 100644
--- a/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
+++ b/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
@@ -32,8 +32,8 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
-#include "incompressible/RASModel/RASModel.H"
+#include "singlePhaseTransportModel.H"
+#include "RASModel.H"
 #include "MRFZones.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/Make/options b/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/Make/options
index 2a417f9549..1223bdd06f 100644
--- a/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/Make/options
+++ b/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/Make/options
@@ -1,7 +1,9 @@
 EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/RAS \
-    -I$(LIB_SRC)/transportModels
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
     -lincompressibleRASModels \
diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/makeMesh b/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/makeMesh
index f45c6ff0be..b0be10d8e5 100755
--- a/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/makeMesh
+++ b/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/makeMesh
@@ -3,8 +3,9 @@
 m4 < constant/polyMesh/blockMeshDict.m4 > constant/polyMesh/blockMeshDict
 blockMesh
 cellSet
-cp system/faceSetDict_rotorFaces system/faceSetDict
-faceSet
-cp system/faceSetDict_noBoundaryFaces system/faceSetDict
-faceSet
+#- MRF determines its own faceZone if not supplied
+#cp system/faceSetDict_rotorFaces system/faceSetDict
+#faceSet
+#cp system/faceSetDict_noBoundaryFaces system/faceSetDict
+#faceSet
 setsToZones -noFlipMap
diff --git a/tutorials/incompressible/channelFoam/channel395/constant/postChannelDict b/tutorials/incompressible/channelFoam/channel395/constant/postChannelDict
index 89fcd50809..96bbe30d76 100644
--- a/tutorials/incompressible/channelFoam/channel395/constant/postChannelDict
+++ b/tutorials/incompressible/channelFoam/channel395/constant/postChannelDict
@@ -15,10 +15,14 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Seed patches to start layering from
 patches         ( bottomWall );
 
+// Direction in which the layers are
 component       y;
 
+// Is the mesh symmetric? If so average(symmetric fields) or
+// subtract(asymmetric) contributions from both halves
 symmetric       true;
 
 
diff --git a/tutorials/incompressible/channelFoam/channel395/system/controlDict b/tutorials/incompressible/channelFoam/channel395/system/controlDict
index 2ca28c29a1..218a76d861 100644
--- a/tutorials/incompressible/channelFoam/channel395/system/controlDict
+++ b/tutorials/incompressible/channelFoam/channel395/system/controlDict
@@ -43,7 +43,35 @@ timePrecision   6;
 
 runTimeModifiable yes;
 
-functions       ( fieldAverage1 { type fieldAverage ; functionObjectLibs ( "libfieldFunctionObjects.so" ) ; enabled true ; outputControl outputTime ; fields ( U { mean on ; prime2Mean on ; base time ; } p { mean on ; prime2Mean on ; base time ; } ) ; } );
+functions
+(
+    fieldAverage1
+    {
+        type fieldAverage;
+
+        functionObjectLibs ( "libfieldFunctionObjects.so" );
+
+        enabled true;
+
+        outputControl outputTime;
+
+        fields
+        (
+            U
+            {
+                mean on;
+                prime2Mean on;
+                base time;
+            }
+            p
+            {
+                mean on;
+                prime2Mean on;
+                base time;
+            }
+        );
+    }
+);
 
 
 // ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/iglooWithFridges/system/snappyHexMeshDict b/tutorials/mesh/snappyHexMesh/iglooWithFridges/system/snappyHexMeshDict
index 2efe6901aa..19ef2afec1 100644
--- a/tutorials/mesh/snappyHexMesh/iglooWithFridges/system/snappyHexMeshDict
+++ b/tutorials/mesh/snappyHexMesh/iglooWithFridges/system/snappyHexMeshDict
@@ -2,7 +2,7 @@
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
 |  \\    /   O peration     | Version:  1.5                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
@@ -10,138 +10,373 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "system";
-    object      snappyHexMeshDict;
+    object      autoHexMeshDict;
 }
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Which of the steps to run
 castellatedMesh true;
-
 snap            true;
-
 addLayers       true;
 
+
+// Geometry. Definition of all surfaces. All surfaces are of class
+// searchableSurface.
+// Surfaces are used
+// - to specify refinement for any mesh cell intersecting it
+// - to specify refinement for any mesh cell inside/outside/near
+// - to 'snap' the mesh boundary to the surface
 geometry
 {
-    fridgeA
-    {
-        type            searchableBox;
-        min             ( 2 2 0 );
-        max             ( 3 3 2 );
-    }
-
-    fridgeB
-    {
-        type            searchableBox;
-        min             ( 3.5 3 0 );
-        max             ( 4.3 3.8 1.8 );
-    }
-
     igloo
     {
-        type            searchableSphere;
-        centre          ( 3 3 0 );
-        radius          4;
+        type searchableSphere;
+        centre (3 3 0);
+        radius 4;
     }
-}
 
+    box1
+    {
+        type searchableBox;
+        min (0 0 0);
+        max (1 1 1);
+    }
+    fridgeFreezer
+    {
+        type searchableSurfaceCollection;
+    
+        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;
+    
+        seal
+        {
+            surface fridgeFreezer;
+            scale (1.0 1.0 1.0);
+            transform
+            {
+                type    cartesian;
+                origin  (2 2 0);
+                e1      (1 0 0);
+                e3      (0 0 1);
+            }
+        }
+        herring
+        {
+            surface fridgeFreezer;
+            scale (1.0 1.0 1.0);
+            transform
+            {
+                type    cartesian;
+                origin  (3.5 3 0);
+                e1      (1 0 0);
+                e3      (0 0 1);
+            }
+        }
+    }
+};
+
+
+
+// Settings for the castellatedMesh generation.
 castellatedMeshControls
 {
-    maxLocalCells   1000000;
-    maxGlobalCells  2000000;
-    minRefinementCells 0;
+
+    // Refinement parameters
+    // ~~~~~~~~~~~~~~~~~~~~~
+
+    // While refining maximum number of cells per processor. This is basically
+    // the number of cells that fit on a processor. If you choose this too small
+    // it will do just more refinement iterations to obtain a similar mesh.
+    maxLocalCells 1000000;
+
+    // Overall cell limit (approximately). Refinement will stop immediately
+    // upon reaching this number so a refinement level might not complete.
+    // Note that this is the number of cells before removing the part which
+    // is not 'visible' from the keepPoint. The final number of cells might
+    // actually be a lot less.
+    maxGlobalCells 2000000;
+
+    // The surface refinement loop might spend lots of iterations refining just a
+    // few cells. This setting will cause refinement to stop if <= minimumRefine
+    // are selected for refinement. Note: it will at least do one iteration
+    // (unless the number of cells to refine is 0)
+    minRefinementCells 100;
+
+    // Number of buffer layers between different levels.
+    // 1 means normal 2:1 refinement restriction, larger means slower
+    // refinement.
     nCellsBetweenLevels 1;
-    features        ( );
+
+
+
+    // Explicit feature edge refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies a level for any cell intersected by its edges.
+    // This is a featureEdgeMesh, read from constant/triSurface for now.
+    features
+    (
+//        {
+//            file "fridgeA.eMesh";
+//            level 3;
+//        }
+    );
+
+
+
+    // Surface based refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies two levels for every surface. The first is the minimum level,
+    // every cell intersecting a surface gets refined up to the minimum level.
+    // The second level is the maximum level. Cells that 'see' multiple
+    // intersections where the intersections make an
+    // angle > resolveFeatureAngle get refined up to the maximum level.
+
     refinementSurfaces
     {
-        fridgeA
+        twoFridgeFreezers
         {
-            level           ( 2 2 );
+            // Surface-wise min and max refinement level
+            level (2 2);
+
+            regions
+            {
+                // Region-wise override
+                "cook.*"
+                {
+                    level (3 3);
+                }
+            }
         }
 
-        fridgeB
+        "iglo.*"
         {
-            level           ( 2 2 );
-        }
-
-        igloo
-        {
-            level           ( 1 1 );
+            // Surface-wise min and max refinement level
+            level (1 1);
         }
     }
 
+    // Resolve sharp angles on fridges
     resolveFeatureAngle 60;
+
+
+    // Region-wise refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies refinement level for cells in relation to a surface. One of
+    // three modes
+    // - distance. 'levels' specifies per distance to the surface the
+    //   wanted refinement level. The distances need to be specified in
+    //   descending order.
+    // - inside. 'levels' is only one entry and only the level is used. All
+    //   cells inside the surface get refined up to the level. The surface
+    //   needs to be closed for this to be possible.
+    // - outside. Same but cells outside.
+
     refinementRegions
     {
     }
 
-    locationInMesh  ( 3 0.28 0.43 );
+
+    // Mesh selection
+    // ~~~~~~~~~~~~~~
+
+    // After refinement patches get added for all refinementSurfaces and
+    // all cells intersecting the surfaces get put into these patches. The
+    // section reachable from the locationInMesh is kept.
+    // NOTE: This point should never be on a face, always inside a cell, even
+    // after refinement.
+    locationInMesh (3 0.28 0.43);
 }
 
+
+
+// Settings for the snapping.
 snapControls
 {
-    nSmoothPatch    3;
-    tolerance       4;
-    nSolveIter      30;
-    nRelaxIter      5;
+    //- Number of patch smoothing iterations before finding correspondence
+    //  to surface
+    nSmoothPatch 3;
+
+    //- Relative distance for points to be attracted by surface feature point
+    //  or edge. True distance is this factor times local
+    //  maximum edge length.
+    tolerance 4.0;
+
+    //- Number of mesh displacement relaxation iterations.
+    nSolveIter 30;
+
+    //- Maximum number of snapping relaxation iterations. Should stop
+    //  before upon reaching a correct mesh.
+    nRelaxIter 5;
 }
 
+
+
+// Settings for the layer addition.
 addLayersControls
 {
+    // Per final patch (so not geometry!) the layer information
     layers
     {
-        fridgeA_region0
+        "two.*"
         {
-            nSurfaceLayers  1;
+            nSurfaceLayers 3;
         }
-
-        fridgeB_region0
+        "igloo_.*"
         {
-            nSurfaceLayers  1;
-        }
-
-        igloo_region0
-        {
-            nSurfaceLayers  1;
+            nSurfaceLayers 1;
         }
     }
 
-    expansionRatio  1;
+    // Expansion factor for layer mesh
+    expansionRatio 1.0;
+
+    //- Wanted thickness of final added cell layer. If multiple layers
+    //  is the
+    //  thickness of the layer furthest away from the wall.
+    //  Relative to undistorted size of cell outside layer.
     finalLayerRatio 0.5;
-    minThickness    0.25;
-    nGrow           0;
-    featureAngle    60;
-    nRelaxIter      5;
+
+    //- Minimum thickness of cell layer. If for any reason layer
+    //  cannot be above minThickness do not add layer.
+    //  Relative to undistorted size of cell outside layer.
+    minThickness 0.25;
+
+    //- If points get not extruded do nGrow layers of connected faces that are
+    //  also not grown. This helps convergence of the layer addition process
+    //  close to features.
+    nGrow 0;
+
+
+    // Advanced settings
+
+    //- When not to extrude surface. 0 is flat surface, 90 is when two faces
+    //  make straight angle.
+    featureAngle 60;
+
+    //- Maximum number of snapping relaxation iterations. Should stop
+    //  before upon reaching a correct mesh.
+    nRelaxIter 5;
+
+    // Number of smoothing iterations of surface normals
     nSmoothSurfaceNormals 1;
-    nSmoothNormals  3;
+
+    // Number of smoothing iterations of interior mesh movement direction
+    nSmoothNormals 3;
+
+    // Smooth layer thickness over surface patches
     nSmoothThickness 10;
+
+    // Stop layer growth on highly warped cells
     maxFaceThicknessRatio 0.5;
+
+    // Reduce layer growth where ratio thickness to medial
+    // distance is large
     maxThicknessToMedialRatio 0.3;
+
+    // Angle used to pick up medial axis points
     minMedianAxisAngle 130;
+
+    // Create buffer region for new layer terminations
     nBufferCellsNoExtrude 0;
 }
 
+
+
+// Generic mesh quality settings. At any undoable phase these determine
+// where to undo.
 meshQualityControls
 {
-    maxNonOrtho     65;
+    //- Maximum non-orthogonality allowed. Set to 180 to disable.
+    maxNonOrtho 65;
+
+    //- Max skewness allowed. Set to <0 to disable.
     maxBoundarySkewness 20;
     maxInternalSkewness 4;
-    maxConcave      80;
-    minFlatness     0.5;
-    minVol          1e-13;
-    minArea         -1;
-    minTwist        0.05;
-    minDeterminant  0.001;
-    minFaceWeight   0.05;
-    minVolRatio     0.01;
+
+    //- Max concaveness allowed. Is angle (in degrees) below which concavity
+    //  is allowed. 0 is straight face, <0 would be convex face.
+    //  Set to 180 to disable.
+    maxConcave 80;
+
+    //- Minimum projected area v.s. actual area. Set to -1 to disable.
+    minFlatness 0.5;
+
+    //- Minimum pyramid volume. Is absolute volume of cell pyramid.
+    //  Set to very negative number (e.g. -1E30) to disable.
+    minVol 1e-13;
+
+    //- Minimum face area. Set to <0 to disable.
+    minArea -1;
+
+    //- Minimum face twist. Set to <-1 to disable. dot product of face normal
+    //- and face centre triangles normal
+    minTwist 0.05;
+
+    //- minimum normalised cell determinant
+    //- 1 = hex, <= 0 = folded or flattened illegal cell
+    minDeterminant 0.001;
+
+    //- minFaceWeight (0 -> 0.5)
+    minFaceWeight 0.05;
+
+    //- minVolRatio (0 -> 1)
+    minVolRatio 0.01;
+
+    //must be >0 for Fluent compatibility
     minTriangleTwist -1;
-    nSmoothScale    4;
-    errorReduction  0.75;
+
+
+    // Advanced
+
+    //- Number of error distribution iterations
+    nSmoothScale 4;
+    //- amount to scale back displacement at error points
+    errorReduction 0.75;
 }
 
-debug           0;
 
-mergeTolerance  1e-06;
+// Advanced
 
+// Flags for optional output
+// 0 : only write final meshes
+// 1 : write intermediate meshes
+// 2 : write volScalarField with cellLevel for postprocessing
+// 4 : write current intersections as .obj files
+debug 0;
+
+
+// Merge tolerance. Is fraction of overall bounding box of initial mesh.
+// Note: the write tolerance needs to be higher than this.
+mergeTolerance 1E-6;
 
 // ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/motorBike/system/snappyHexMeshDict b/tutorials/mesh/snappyHexMesh/motorBike/system/snappyHexMeshDict
index fa094e1b03..3f0ff0417e 100644
--- a/tutorials/mesh/snappyHexMesh/motorBike/system/snappyHexMeshDict
+++ b/tutorials/mesh/snappyHexMesh/motorBike/system/snappyHexMeshDict
@@ -13,14 +13,21 @@ FoamFile
     location    "system";
     object      snappyHexMeshDict;
 }
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Which of the steps to run
 castellatedMesh true;
-
 snap            true;
-
 addLayers       true;
 
+
+// Geometry. Definition of all surfaces. All surfaces are of class
+// searchableSurface.
+// Surfaces are used
+// - to specify refinement for any mesh cell intersecting it
+// - to specify refinement for any mesh cell inside/outside/near
+// - to 'snap' the mesh boundary to the surface
 geometry
 {
     motorBike.stl
@@ -37,22 +44,89 @@ geometry
     }
 }
 
+
+// Settings for the castellatedMesh generation.
 castellatedMeshControls
 {
+
+    // Refinement parameters
+    // ~~~~~~~~~~~~~~~~~~~~~
+
+    // While refining maximum number of cells per processor. This is basically
+    // the number of cells that fit on a processor. If you choose this too small
+    // it will do just more refinement iterations to obtain a similar mesh.
     maxLocalCells   1000000;
+
+
+    // Overall cell limit (approximately). Refinement will stop immediately
+    // upon reaching this number so a refinement level might not complete.
+    // Note that this is the number of cells before removing the part which
+    // is not 'visible' from the keepPoint. The final number of cells might
+    // actually be a lot less.
     maxGlobalCells  2000000;
+
+    // The surface refinement loop might spend lots of iterations refining just a
+    // few cells. This setting will cause refinement to stop if <= minimumRefine
+    // are selected for refinement. Note: it will at least do one iteration
+    // (unless the number of cells to refine is 0)
     minRefinementCells 10;
+
+    // Number of buffer layers between different levels.
+    // 1 means normal 2:1 refinement restriction, larger means slower
+    // refinement.
     nCellsBetweenLevels 2;
-    features        ( );
+
+
+
+    // Explicit feature edge refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies a level for any cell intersected by its edges.
+    // This is a featureEdgeMesh, read from constant/triSurface for now.
+    features
+    (
+        //{
+        //    file "someLine.eMesh";
+        //    level 2;
+        //}
+    );
+
+
+
+    // Surface based refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies two levels for every surface. The first is the minimum level,
+    // every cell intersecting a surface gets refined up to the minimum level.
+    // The second level is the maximum level. Cells that 'see' multiple
+    // intersections where the intersections make an
+    // angle > resolveFeatureAngle get refined up to the maximum level.
+
     refinementSurfaces
     {
         motorBike
         {
+            // Surface-wise min and max refinement level
             level           ( 5 6 );
         }
     }
 
+    // Resolve sharp angles
     resolveFeatureAngle 30;
+
+
+    // Region-wise refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies refinement level for cells in relation to a surface. One of
+    // three modes
+    // - distance. 'levels' specifies per distance to the surface the
+    //   wanted refinement level. The distances need to be specified in
+    //   descending order.
+    // - inside. 'levels' is only one entry and only the level is used. All
+    //   cells inside the surface get refined up to the level. The surface
+    //   needs to be closed for this to be possible.
+    // - outside. Same but cells outside.
     refinementRegions
     {
         refinementBox
@@ -62,19 +136,45 @@ castellatedMeshControls
         }
     }
 
+
+
+    // Mesh selection
+    // ~~~~~~~~~~~~~~
+
+    // After refinement patches get added for all refinementSurfaces and
+    // all cells intersecting the surfaces get put into these patches. The
+    // section reachable from the locationInMesh is kept.
+    // NOTE: This point should never be on a face, always inside a cell, even
+    // after refinement.
     locationInMesh  ( 3 3 0.43 );
 }
 
+
+// Settings for the snapping.
 snapControls
 {
+    //- Number of patch smoothing iterations before finding correspondence
+    //  to surface
     nSmoothPatch    3;
+
+    //- Relative distance for points to be attracted by surface feature point
+    //  or edge. True distance is this factor times local
+    //  maximum edge length.
     tolerance       4;
+
+    //- Number of mesh displacement relaxation iterations.
     nSolveIter      30;
+
+    //- Maximum number of snapping relaxation iterations. Should stop
+    //  before upon reaching a correct mesh.
     nRelaxIter      5;
 }
 
+
+// Settings for the layer addition.
 addLayersControls
 {
+    // Per final patch (so not geometry!) the layer information
     layers
     {
         minZ
@@ -418,41 +518,126 @@ addLayersControls
         }
     }
 
+    // Expansion factor for layer mesh
     expansionRatio  1;
+
+    //- Wanted thickness of final added cell layer. If multiple layers
+    //  is the
+    //  thickness of the layer furthest away from the wall.
+    //  Relative to undistorted size of cell outside layer.
     finalLayerRatio 0.3;
+
+    //- Minimum thickness of cell layer. If for any reason layer
+    //  cannot be above minThickness do not add layer.
+    //  Relative to undistorted size of cell outside layer.
     minThickness    0.1;
+
+    //- If points get not extruded do nGrow layers of connected faces that are
+    //  also not grown. This helps convergence of the layer addition process
+    //  close to features.
     nGrow           1;
+
+
+    // Advanced settings
+
+    //- When not to extrude surface. 0 is flat surface, 90 is when two faces
+    //  make straight angle.
     featureAngle    30;
+
+    //- Maximum number of snapping relaxation iterations. Should stop
+    //  before upon reaching a correct mesh.
     nRelaxIter      3;
+
+    // Number of smoothing iterations of surface normals
     nSmoothSurfaceNormals 1;
+
+    // Number of smoothing iterations of interior mesh movement direction
     nSmoothNormals  3;
+
+    // Smooth layer thickness over surface patches
     nSmoothThickness 10;
+
+    // Stop layer growth on highly warped cells
     maxFaceThicknessRatio 0.5;
+
+    // Reduce layer growth where ratio thickness to medial
+    // distance is large
     maxThicknessToMedialRatio 0.3;
+
+    // Angle used to pick up medial axis points
     minMedianAxisAngle 130;
+
+    // Create buffer region for new layer terminations
     nBufferCellsNoExtrude 0;
 }
 
+
+// Generic mesh quality settings. At any undoable phase these determine
+// where to undo.
 meshQualityControls
 {
+    //- Maximum non-orthogonality allowed. Set to 180 to disable.
     maxNonOrtho     65;
+
+    //- Max skewness allowed. Set to <0 to disable.
     maxBoundarySkewness 20;
     maxInternalSkewness 4;
+
+    //- Max concaveness allowed. Is angle (in degrees) below which concavity
+    //  is allowed. 0 is straight face, <0 would be convex face.
+    //  Set to 180 to disable.
     maxConcave      80;
+
+    //- Minimum projected area v.s. actual area. Set to -1 to disable.
     minFlatness     0.5;
+
+    //- Minimum pyramid volume. Is absolute volume of cell pyramid.
+    //  Set to very negative number (e.g. -1E30) to disable.
     minVol          1e-13;
+
+    //- Minimum face area. Set to <0 to disable.
     minArea         -1;
+
+    //- Minimum face twist. Set to <-1 to disable. dot product of face normal
+    //- and face centre triangles normal
     minTwist        0.02;
+
+    //- minimum normalised cell determinant
+    //- 1 = hex, <= 0 = folded or flattened illegal cell
     minDeterminant  0.001;
+
+    //- minFaceWeight (0 -> 0.5)
     minFaceWeight   0.02;
+
+    //- minVolRatio (0 -> 1)
     minVolRatio     0.01;
+
+    //must be >0 for Fluent compatibility
     minTriangleTwist -1;
+
+
+    // Advanced
+
+    //- Number of error distribution iterations
     nSmoothScale    4;
+    //- amount to scale back displacement at error points
     errorReduction  0.75;
 }
 
+
+
+// Advanced
+
+// Flags for optional output
+// 0 : only write final meshes
+// 1 : write intermediate meshes
+// 2 : write volScalarField with cellLevel for postprocessing
+// 4 : write current intersections as .obj files
 debug           0;
 
+
+// Merge tolerance. Is fraction of overall bounding box of initial mesh.
+// Note: the write tolerance needs to be higher than this.
 mergeTolerance  1e-06;
 
 
diff --git a/wmake/wmakeLnInclude b/wmake/wmakeLnInclude
index bce2cfbb51..1ed1671e7a 100755
--- a/wmake/wmakeLnInclude
+++ b/wmake/wmakeLnInclude
@@ -145,16 +145,6 @@ find -L . -type l -exec rm \{\} \;
 find .. $findOpt \
     \( -name lnInclude -o -name Make -o -name config -o -name noLink \) -prune \
  -o \( -name '*.[CHh]' -o -name '*.[ch]xx' -o -name '*.[ch]pp' -o -name '*.type' \)  \
- -a ! -name ".#*" \
- -print | \
-    while read src
-    do
-        link=$(readlink ${src##*/})
-        if [ "$link" != "$src" ]
-        then
-            rm $link 2>/dev/null
-            ln $lnOpt $src .
-        fi
-    done
+ -exec ln $lnOpt {} . \;
 
 #------------------------------------------------------------------------------