diff --git a/.gitignore b/.gitignore
index 96c3093b718b9a8ee6d711583018ae836555158a..150288aae2555094adff1f866bb0e4e950f07c9f 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 61091811bf04fd0547f08d2b4cd5c3c9385d4d6d..1223bdd06f48178f0b2eee8032d29c936634f328 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/bin/foamEbrowse b/bin/foamEbrowse
index e303a772643c8b9f74fa5bc60bc565e110a90a20..66fafc1e2b1b46f5384ca1c78f309fb1b785559d 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 0000000000000000000000000000000000000000..1cfd186f52b1b32e1d9d8c77fce0d44e5218d507
--- /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 0000000000000000000000000000000000000000..ce66e0c084a8188bf63b707880344ccc18813143
--- /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 f22e0729dbe8e50694ee43a373205494b0a2a536..c3ce648ba0aa1462b53dff9a1502db02d5a215aa 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/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C
index 666d84a5c1206588bce906d4174fe31a8857cda1..05958edf68ad060e63c51e8f078656fbbb82a5e6 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/finiteVolume/cfdTools/general/MRF/MRFZone.C b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
index 8bc4d3b938109171a57dbb8513f7a4d6ec753537..ab9709e993871b6aeaf16138e3f13408b664a591 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 3ffb4e461a853a7703ee7a5a9038bc99893c5efb..7c489ed25bba06d0e9e44334ff609ced8fbd0f57 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/meshTools/cellDist/wallPoint/wallPoint.H b/src/meshTools/cellDist/wallPoint/wallPoint.H
index 0544b6ce361ea8d9f4143987703f83765b3996cf..82d6b385af96d0efe538574416ec70a3d41f579d 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/searchableSurface/distributedTriSurfaceMesh.C b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
index 870f7e8324bf3665463790e98c2dd72f704fa94d..7f86c0f6581fcd94504b9cb7d39fdc67284a4506 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,12 +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(io),
+    triSurfaceMesh
+    (
+        IOobject
+        (
+            io.name(),
+            io.time().findInstance(io.local(), word::null),
+            io.local(),
+            io.db(),
+            io.readOpt(),
+            io.writeOpt(),
+            io.registerObject()
+        )
+    ),
     dict_
     (
         IOobject
@@ -1362,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;
+    }
 }
 
 
@@ -1371,7 +1421,21 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
     const dictionary& dict
 )
 :
-    triSurfaceMesh(io, 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
@@ -1387,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;
+    }
 }
 
 
@@ -2389,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 733cd24891e21befee9e9318fa301180bc3f8966..4d5f7dfd173edd6ccb3c559221de0564d46edc50 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 ba5a1888cfb343f39605b49b8111083427d7562c..287ae6743c2bb508c0ef64322036fac617123dcd 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
index 80666ef1f96d9541f099a2c6568f3a8f52094502..ce6ddd643122703e9214fc6041a5a76265b098f4 100644
--- a/src/meshTools/searchableSurface/searchableSurfaceCollection.C
+++ b/src/meshTools/searchableSurface/searchableSurfaceCollection.C
@@ -59,43 +59,44 @@ void Foam::searchableSurfaceCollection::findNearest
 
     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
-            ),
+            transform_[surfI].localPosition(samples),
             scale_[surfI]
         );
-        
-        subGeom_[surfI].findNearest(localSamples, minDistSqr, hitInfo);
+
+        subGeom_[surfI].findNearest(localSamples, localMinDistSqr, hitInfo);
 
         forAll(hitInfo, pointI)
         {
             if (hitInfo[pointI].hit())
             {
-                minDistSqr[pointI] = magSqr
-                (   
-                    hitInfo[pointI].hitPoint()
-                  - localSamples()[pointI]
-                );
-
                 // Rework back into global coordinate sys. Multiply then
                 // transform
-                nearestInfo[pointI] = hitInfo[pointI];
-                nearestInfo[pointI].rawPoint() = 
-                    transform_[surfI].globalPosition
+                point globalPt = transform_[surfI].globalPosition
+                (
+                    cmptMultiply
                     (
-                        cmptMultiply
-                        (
-                            nearestInfo[pointI].rawPoint(),
-                            scale_[surfI]
-                        )
-                    );
-                nearestSurf[pointI] = surfI;
+                        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;
+                }
             }
         }
     }
diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.C b/src/meshTools/searchableSurface/triSurfaceMesh.C
index 121762d7c4e219fc05f570535ab903d5eccaa87a..9958afd1b155167a1e08bc9f3926f94cad434532 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.C
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.C
@@ -228,19 +228,19 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
 Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
 :
     // Find instance for triSurfaceMesh
-    searchableSurface
-    (
-        IOobject
-        (
-            io.name(),
-            io.time().findInstance(io.local(), word::null),
-            io.local(),
-            io.db(),
-            io.readOpt(),
-            io.writeOpt(),
-            io.registerObject()
-        )
-    ),
+    searchableSurface(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
     (
@@ -273,19 +273,19 @@ Foam::triSurfaceMesh::triSurfaceMesh
     const dictionary& dict
 )
 :
-    searchableSurface
-    (
-        IOobject
-        (
-            io.name(),
-            io.time().findInstance(io.local(), word::null),
-            io.local(),
-            io.db(),
-            io.readOpt(),
-            io.writeOpt(),
-            io.registerObject()
-        )
-    ),
+    searchableSurface(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
     (
diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.H b/src/meshTools/searchableSurface/triSurfaceMesh.H
index d7aa33cfca43059c1cf5154bd231408b07e37b2c..a4231457d5183097cff2eceb31ca00e944c6a05d 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.H
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.H
@@ -113,11 +113,10 @@ public:
         //- Construct from triSurface
         triSurfaceMesh(const IOobject&, const triSurface&);
 
-        //- Construct read. Does timeInstance search.
+        //- Construct read.
         triSurfaceMesh(const IOobject& io);
 
         //- Construct from IO and dictionary (used by searchableSurface).
-        //  Does timeInstance search.
         //  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 939d77cf1b8e7f536f5e8f20865bb66c160f101c..78081050ab4e715d3ec6f4378e3f91e1e91dc8da 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C
@@ -199,81 +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)
         {
-            // Still needs special treatment?
-
-            forAll(pp, i)
-            {
-                bool ownLower = (cVals[own[faceI]] < iso_);
+            bool ownLower = (cVals[own[faceI]] < iso_);
 
-                scalar nbrValue;
-                point nbrPoint;
-                getNeighbour
-                (
-                    boundaryRegion,
-                    cVals,
-                    own[faceI],
-                    faceI,
-                    nbrValue,
-                    nbrPoint
-                );
+            scalar nbrValue;
+            point nbrPoint;
+            getNeighbour
+            (
+                boundaryRegion,
+                cVals,
+                own[faceI],
+                faceI,
+                nbrValue,
+                nbrPoint
+            );
 
-                bool neiLower = (nbrValue < iso_);
+            bool neiLower = (nbrValue < 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++;
             }
-        }
-        else
-        {
-            forAll(pp, i)
-            {
-                bool ownLower = (cVals[own[faceI]] < iso_);
 
-                scalar nbrValue;
-                point nbrPoint;
-                getNeighbour
-                (
-                    boundaryRegion,
-                    cVals,
-                    own[faceI],
-                    faceI,
-                    nbrValue,
-                    nbrPoint
-                );
-
-                bool neiLower = (nbrValue < 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++;
         }
     }
 
@@ -1589,26 +1555,6 @@ Foam::isoSurface::isoSurface
 
     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);
-        }
-    }
-
     // Pre-calculate patch-per-face to avoid whichPatch call.
     labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces());
 
@@ -1724,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
@@ -1746,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&>
@@ -1758,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());
+            }
         }
     }
 
@@ -1885,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 f6c022c0362f06fdd40c69f0a84d76b9957dcb36..b859524b2cc40f21aad1351b3dfa2d6f4d798a88 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
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
index 7e7206edcc77264aa9109a5f901230e3ec60d63b..0760b06a52561de62d50ceb68b00f79b4b016fe7 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
@@ -514,9 +514,53 @@ void Foam::isoSurface::generateTriPoints
         }
         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)
@@ -534,8 +578,8 @@ void Foam::isoSurface::generateTriPoints
                         snappedPoint,
                         faceI,
 
-                        cVals[own[faceI]],
-                        cCoords.boundaryField()[patchI][i],
+                        bVals[i],
+                        bCoords[i],
                         false,  // fc not snapped
                         pTraits<Type>::zero,
 
diff --git a/src/turbulenceModels/LES/Allwmake b/src/turbulenceModels/LES/Allwmake
index a80c71ab57d28206f7585707735d49f8be229d95..d83439a9130c64ad25c955e5059ee79e1206c5d8 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/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C b/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
index 004a92b45f5c92a81fc537e826508bc7b82bfee0..2029bc066a22712e55bad739612cb3734f2514c5 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 2a417f954936f537595fba5fa461a73a81650574..1223bdd06f48178f0b2eee8032d29c936634f328 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 f45c6ff0be3a4fe281484bbfd4789baf4fe00b32..b0be10d8e5f2b0525cec8f23c93585fde73f5879 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/mesh/snappyHexMesh/iglooWithFridges/system/snappyHexMeshDict b/tutorials/mesh/snappyHexMesh/iglooWithFridges/system/snappyHexMeshDict
index d185a59ce92d82eea0c026aafae488fce4e709e5..19ef2afec1d08e169d49fd7907cc518d0eccd94b 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,9 +10,9 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "system";
-    object      snappyHexMeshDict;
+    object      autoHexMeshDict;
 }
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // Which of the steps to run
@@ -29,27 +29,78 @@ addLayers       true;
 // - to 'snap' the mesh boundary to the surface
 geometry
 {
-    fridgeA
+    igloo
     {
-        type            searchableBox;
-        min             ( 2 2 0 );
-        max             ( 3 3 2 );
+        type searchableSphere;
+        centre (3 3 0);
+        radius 4;
     }
 
-    fridgeB
+    box1
     {
-        type            searchableBox;
-        min             ( 3.5 3 0 );
-        max             ( 4.3 3.8 1.8 );
+        type searchableBox;
+        min (0 0 0);
+        max (1 1 1);
     }
-
-    igloo
+    fridgeFreezer
     {
-        type            searchableSphere;
-        centre          ( 3 3 0 );
-        radius          4;
+        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);
+            }
+        }
+    }
+};
 
 
 
@@ -63,20 +114,20 @@ castellatedMeshControls
     // 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;
+    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;
+    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 0;
+    minRefinementCells 100;
 
     // Number of buffer layers between different levels.
     // 1 means normal 2:1 refinement restriction, larger means slower
@@ -92,10 +143,10 @@ castellatedMeshControls
     // This is a featureEdgeMesh, read from constant/triSurface for now.
     features
     (
-        //{
-        //    file "someLine.eMesh";
-        //    level 2;
-        //}
+//        {
+//            file "fridgeA.eMesh";
+//            level 3;
+//        }
     );
 
 
@@ -108,27 +159,35 @@ castellatedMeshControls
     // 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
         {
             // Surface-wise min and max refinement level
-            level           ( 2 2 );
+            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
     // ~~~~~~~~~~~~~~~~~~~~~~
 
@@ -155,7 +214,7 @@ castellatedMeshControls
     // 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 );
+    locationInMesh (3 0.28 0.43);
 }
 
 
@@ -165,19 +224,19 @@ snapControls
 {
     //- Number of patch smoothing iterations before finding correspondence
     //  to surface
-    nSmoothPatch    3;
+    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;
+    tolerance 4.0;
 
     //- Number of mesh displacement relaxation iterations.
-    nSolveIter      30;
+    nSolveIter 30;
 
     //- Maximum number of snapping relaxation iterations. Should stop
     //  before upon reaching a correct mesh.
-    nRelaxIter      5;
+    nRelaxIter 5;
 }
 
 
@@ -188,25 +247,18 @@ addLayersControls
     // Per final patch (so not geometry!) the layer information
     layers
     {
-        fridgeA_region0
+        "two.*"
         {
-            nSurfaceLayers  1;
+            nSurfaceLayers 3;
         }
-
-        fridgeB_region0
-        {
-            nSurfaceLayers  1;
-        }
-
-        igloo_region0
+        "igloo_.*"
         {
-            nSurfaceLayers  1;
+            nSurfaceLayers 1;
         }
     }
 
-
     // Expansion factor for layer mesh
-    expansionRatio  1;
+    expansionRatio 1.0;
 
     //- Wanted thickness of final added cell layer. If multiple layers
     //  is the
@@ -217,29 +269,29 @@ addLayersControls
     //- 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;
+    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;
+    nGrow 0;
 
 
     // Advanced settings
 
     //- When not to extrude surface. 0 is flat surface, 90 is when two faces
     //  make straight angle.
-    featureAngle    60;
+    featureAngle 60;
 
     //- Maximum number of snapping relaxation iterations. Should stop
     //  before upon reaching a correct mesh.
-    nRelaxIter      5;
+    nRelaxIter 5;
 
     // Number of smoothing iterations of surface normals
     nSmoothSurfaceNormals 1;
 
     // Number of smoothing iterations of interior mesh movement direction
-    nSmoothNormals  3;
+    nSmoothNormals 3;
 
     // Smooth layer thickness over surface patches
     nSmoothThickness 10;
@@ -265,7 +317,7 @@ addLayersControls
 meshQualityControls
 {
     //- Maximum non-orthogonality allowed. Set to 180 to disable.
-    maxNonOrtho     65;
+    maxNonOrtho 65;
 
     //- Max skewness allowed. Set to <0 to disable.
     maxBoundarySkewness 20;
@@ -274,32 +326,31 @@ meshQualityControls
     //- 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;
+    maxConcave 80;
 
     //- Minimum projected area v.s. actual area. Set to -1 to disable.
-    minFlatness     0.5;
+    minFlatness 0.5;
 
     //- Minimum pyramid volume. Is absolute volume of cell pyramid.
-    //  Set to a sensible fraction of the smallest cell volume expected.
     //  Set to very negative number (e.g. -1E30) to disable.
-    minVol          1e-13;
+    minVol 1e-13;
 
     //- Minimum face area. Set to <0 to disable.
-    minArea         -1;
+    minArea -1;
 
     //- Minimum face twist. Set to <-1 to disable. dot product of face normal
     //- and face centre triangles normal
-    minTwist        0.05;
+    minTwist 0.05;
 
     //- minimum normalised cell determinant
     //- 1 = hex, <= 0 = folded or flattened illegal cell
-    minDeterminant  0.001;
+    minDeterminant 0.001;
 
     //- minFaceWeight (0 -> 0.5)
-    minFaceWeight   0.05;
+    minFaceWeight 0.05;
 
     //- minVolRatio (0 -> 1)
-    minVolRatio     0.01;
+    minVolRatio 0.01;
 
     //must be >0 for Fluent compatibility
     minTriangleTwist -1;
@@ -308,13 +359,12 @@ meshQualityControls
     // Advanced
 
     //- Number of error distribution iterations
-    nSmoothScale    4;
+    nSmoothScale 4;
     //- amount to scale back displacement at error points
-    errorReduction  0.75;
+    errorReduction 0.75;
 }
 
 
-
 // Advanced
 
 // Flags for optional output
@@ -322,12 +372,11 @@ meshQualityControls
 // 1 : write intermediate meshes
 // 2 : write volScalarField with cellLevel for postprocessing
 // 4 : write current intersections as .obj files
-debug           0;
+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;
-
+mergeTolerance 1E-6;
 
 // ************************************************************************* //
diff --git a/wmake/wmakeLnInclude b/wmake/wmakeLnInclude
index bce2cfbb51d3cfcb74c52a54432dd386bbdafc49..1ed1671e7a32e3043027fe2a680eb3e3031048ca 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 {} . \;
 
 #------------------------------------------------------------------------------