diff --git a/applications/test/dictionary/calcEntry/calcEntry.H b/applications/test/dictionary/calcEntry/calcEntry.H
index e8a0925f94ebd51f3e4226cf0b6ca40f60b0bc7c..119f15564375fa3bb2c31ab1b1c79014745328a9 100644
--- a/applications/test/dictionary/calcEntry/calcEntry.H
+++ b/applications/test/dictionary/calcEntry/calcEntry.H
@@ -64,7 +64,7 @@ class calcEntry
 public:
 
     //- Runtime type information
-    TypeName("calc");
+    ClassName("calc");
 
 
     // Member Functions
diff --git a/applications/utilities/mesh/conversion/writeMeshObj/writeMeshObj.C b/applications/utilities/mesh/conversion/writeMeshObj/writeMeshObj.C
index e1d571e0816f4f146decba57af55579348a03a2f..d2b23fef161b22f9f746251fc4b5c87a02bb4a29 100644
--- a/applications/utilities/mesh/conversion/writeMeshObj/writeMeshObj.C
+++ b/applications/utilities/mesh/conversion/writeMeshObj/writeMeshObj.C
@@ -32,7 +32,8 @@ Description
 
     patch_YYY_XXX.obj : all face centres of patch YYY
 
-    Optional: patch faces (as polygons) : patchFaces_YYY_XXX.obj
+    Optional: - patch faces (as polygons) : patchFaces_YYY_XXX.obj
+              - non-manifold edges : patchEdges_YYY_XXX.obj
 
 \*---------------------------------------------------------------------------*/
 
@@ -51,7 +52,7 @@ using namespace Foam;
 
 void writeOBJ(const point& pt, Ostream& os)
 {
-    os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
+    os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
 }
 
 // All edges of mesh
@@ -75,8 +76,7 @@ void writePoints(const polyMesh& mesh, const fileName& timeName)
     {
         const edge& e = mesh.edges()[edgeI];
 
-        pointStream << "l " << e.start() + 1 << ' ' << e.end() + 1
-            << endl;
+        pointStream << "l " << e.start() + 1 << ' ' << e.end() + 1 << nl;
     }
 }
 
@@ -277,7 +277,47 @@ void writePatchFaces
             {
                 patchFaceStream << ' ' << f[fp]+1;
             }
-            patchFaceStream << endl;
+            patchFaceStream << nl;
+        }
+    }
+}
+
+
+void writePatchBoundaryEdges
+(
+    const polyMesh& mesh,
+    const fileName& timeName
+)
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        fileName edgeFile
+        (
+            mesh.time().path()
+          / "patchEdges_" + pp.name() + '_' + timeName + ".obj"
+        );
+
+        Info << "Writing patch edges to " << edgeFile << endl;
+
+        OFstream patchEdgeStream(edgeFile);
+
+        forAll(pp.localPoints(), pointI)
+        {
+            writeOBJ(pp.localPoints()[pointI], patchEdgeStream);
+        }
+
+        for (label edgeI = pp.nInternalEdges(); edgeI < pp.nEdges(); edgeI++)
+        {
+            if (pp.edgeFaces()[edgeI].size() == 1)
+            {
+                const edge& e = pp.edges()[edgeI];
+
+                patchEdgeStream<< "l " << e[0]+1 << ' ' << e[1]+1 << nl;
+            }
         }
     }
 }
@@ -339,6 +379,7 @@ int main(int argc, char *argv[])
 {
     timeSelector::addOptions();
     argList::validOptions.insert("patchFaces", "");
+    argList::validOptions.insert("patchEdges", "");
     argList::validOptions.insert("cell", "cellI");
     argList::validOptions.insert("face", "faceI");
     argList::validOptions.insert("point", "pointI");
@@ -351,6 +392,7 @@ int main(int argc, char *argv[])
     runTime.functionObjects().off();
 
     bool patchFaces = args.optionFound("patchFaces");
+    bool patchEdges = args.optionFound("patchEdges");
     bool doCell     = args.optionFound("cell");
     bool doPoint    = args.optionFound("point");
     bool doFace     = args.optionFound("face");
@@ -381,19 +423,23 @@ int main(int argc, char *argv[])
             {
                 writePatchFaces(mesh, runTime.timeName());
             }
-            else if (doCell)
+            if (patchEdges)
+            {
+                writePatchBoundaryEdges(mesh, runTime.timeName());
+            }
+            if (doCell)
             {
                 label cellI = args.optionRead<label>("cell");
 
                 writePoints(mesh, cellI, runTime.timeName());
             }
-            else if (doPoint)
+            if (doPoint)
             {
                 label pointI = args.optionRead<label>("point");
 
                 writePointCells(mesh, pointI, runTime.timeName());
             }
-            else if (doFace)
+            if (doFace)
             {
                 label faceI = args.optionRead<label>("face");
 
@@ -415,7 +461,7 @@ int main(int argc, char *argv[])
 
                 meshTools::writeOBJ(str, faceList(1, f), mesh.points());
             }
-            else if (doCellSet)
+            if (doCellSet)
             {
                 word setName(args.option("cellSet"));
 
@@ -427,7 +473,7 @@ int main(int argc, char *argv[])
                 writePoints(mesh, cells.toc(), runTime.timeName());
 
             }
-            else if (doFaceSet)
+            if (doFaceSet)
             {
                 word setName(args.option("faceSet"));
 
@@ -458,7 +504,16 @@ int main(int argc, char *argv[])
                     faces.toc()
                 );
             }
-            else
+            else if
+            (
+                !patchFaces 
+             && !patchEdges
+             && !doCell
+             && !doPoint
+             && !doFace
+             && !doCellSet
+             && !doFaceSet
+            )
             {
                 // points & edges
                 writePoints(mesh, runTime.timeName());
diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeProperties b/applications/utilities/mesh/generation/extrudeMesh/extrudeProperties
index faf875774d5832f029fd5ffc0652da39302afe00..9a047e7822948eaf9012ccffaf12e0418d577e84 100644
--- a/applications/utilities/mesh/generation/extrudeMesh/extrudeProperties
+++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeProperties
@@ -71,7 +71,8 @@ linearNormalCoeffs
 
 linearDirectionCoeffs
 {
-    direction       (0 0 1);
+    direction       (0 1 0);
+    thickness       0.05;
 }
 
 linearRadialCoeffs
@@ -89,7 +90,7 @@ sigmaRadialCoeffs
 
 // Do front and back need to be merged? Usually only makes sense for 360
 // degree wedges.
-mergeFaces true;
+mergeFaces false;   //true;
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/mesh/manipulation/cellSet/cellSet.C b/applications/utilities/mesh/manipulation/cellSet/cellSet.C
index ca8225d3f2d3d0094db068b45a6a405dd2e42d4e..e8036274261ada966abc16b8210fdf36fcf84623 100644
--- a/applications/utilities/mesh/manipulation/cellSet/cellSet.C
+++ b/applications/utilities/mesh/manipulation/cellSet/cellSet.C
@@ -37,37 +37,6 @@ using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Copy set
-void backup
-(
-    const polyMesh& mesh,
-    const word& fromName,
-    const topoSet& fromSet,
-    const word& toName
-)
-{
-    Info<< "Backing up " << fromName << " into " << toName << endl;
-
-    topoSet backupSet(mesh, toName, fromSet);
-
-    backupSet.write();
-}
-
-
-// Read and copy set
-void backup
-(
-    const polyMesh& mesh,
-    const word& fromName,
-    const word& toName
-)
-{
-    topoSet fromSet(mesh, fromName, IOobject::READ_IF_PRESENT);
-
-    backup(mesh, fromName, fromSet, toName);
-}
-
-
 // Main program:
 
 int main(int argc, char *argv[])
@@ -114,8 +83,6 @@ int main(int argc, char *argv[])
     {
         r = IOobject::NO_READ;
 
-        backup(mesh, setName, setName + "_old");
-
         currentSetPtr.reset
         (
             new cellSet
@@ -151,7 +118,7 @@ int main(int argc, char *argv[])
     if ((r == IOobject::MUST_READ) && (action != topoSetSource::LIST))
     {
         // currentSet has been read so can make copy.
-        backup(mesh, setName, currentSet, setName + "_old");
+        //backup(mesh, setName, currentSet, setName + "_old");
     }
 
     if (action == topoSetSource::CLEAR)
@@ -173,7 +140,16 @@ int main(int argc, char *argv[])
         forAll(topoSetSources, topoSetSourceI)
         {
             // Backup current set.
-            topoSet oldSet(mesh, currentSet.name() + "_old2", currentSet);
+            autoPtr<topoSet> oldSet
+            (
+                topoSet::New
+                (
+                    currentSet.type(),
+                    mesh,
+                    currentSet.name() + "_old2",
+                    currentSet
+                )
+            );
 
             currentSet.clear();
 
diff --git a/applications/utilities/mesh/manipulation/faceSet/faceSet.C b/applications/utilities/mesh/manipulation/faceSet/faceSet.C
index d1700cc2f962259e9bb29e65f2ca81e18a60ade9..8441a1cdbe2ce773f4021ad2e572e7de02959b79 100644
--- a/applications/utilities/mesh/manipulation/faceSet/faceSet.C
+++ b/applications/utilities/mesh/manipulation/faceSet/faceSet.C
@@ -37,37 +37,6 @@ using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Copy set
-void backup
-(
-    const polyMesh& mesh,
-    const word& fromName,
-    const topoSet& fromSet,
-    const word& toName
-)
-{
-    Info<< "Backing up " << fromName << " into " << toName << endl;
-
-    topoSet backupSet(mesh, toName, fromSet);
-
-    backupSet.write();
-}
-
-
-// Read and copy set
-void backup
-(
-    const polyMesh& mesh,
-    const word& fromName,
-    const word& toName
-)
-{
-    topoSet fromSet(mesh, fromName, IOobject::READ_IF_PRESENT);
-
-    backup(mesh, fromName, fromSet, toName);
-}
-
-
 // Main program:
 
 int main(int argc, char *argv[])
@@ -114,8 +83,6 @@ int main(int argc, char *argv[])
     {
         r = IOobject::NO_READ;
 
-        backup(mesh, setName, setName + "_old");
-
         currentSetPtr.reset
         (
             new faceSet
@@ -151,7 +118,7 @@ int main(int argc, char *argv[])
     if ((r == IOobject::MUST_READ) && (action != topoSetSource::LIST))
     {
         // currentSet has been read so can make copy.
-        backup(mesh, setName, currentSet, setName + "_old");
+        //backup(mesh, setName, currentSet, setName + "_old");
     }
 
     if (action == topoSetSource::CLEAR)
@@ -173,7 +140,16 @@ int main(int argc, char *argv[])
         forAll(topoSetSources, topoSetSourceI)
         {
             // Backup current set.
-            topoSet oldSet(mesh, currentSet.name() + "_old2", currentSet);
+            autoPtr<topoSet> oldSet
+            (
+                topoSet::New
+                (
+                    currentSet.type(),
+                    mesh,
+                    currentSet.name() + "_old2",
+                    currentSet
+                )
+            );
 
             currentSet.clear();
 
diff --git a/applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C b/applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C
index 28955c4d6f7831cf7ab0d9b20a82d522c931f36d..b13dfb86ef04bc243b54d7517902d8e7a103ce7d 100644
--- a/applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C
+++ b/applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C
@@ -49,7 +49,7 @@ int main(int argc, char *argv[])
     (
         IOobject
         (
-            mergePolyMesh::defaultRegion,
+            masterRegion,
             runTimeMaster.timeName(),
             runTimeMaster
         )
@@ -63,7 +63,7 @@ int main(int argc, char *argv[])
     (
         IOobject
         (
-            mergePolyMesh::defaultRegion,
+            addRegion,
             runTimeToAdd.timeName(),
             runTimeToAdd
         )
diff --git a/applications/utilities/mesh/manipulation/mergeMeshes/setRoots.H b/applications/utilities/mesh/manipulation/mergeMeshes/setRoots.H
index 7ddd408965db14e24049c8d01167154c32bd6b70..cfbf59cf64af2b013f075e3f60b67c5f0bcbea8a 100644
--- a/applications/utilities/mesh/manipulation/mergeMeshes/setRoots.H
+++ b/applications/utilities/mesh/manipulation/mergeMeshes/setRoots.H
@@ -2,9 +2,11 @@
 
     argList::validArgs.append("master root");
     argList::validArgs.append("master case");
+    argList::validOptions.insert("masterRegion", "name");
 
     argList::validArgs.append("root to add");
     argList::validArgs.append("case to add");
+    argList::validOptions.insert("addRegion", "name");
 
     argList args(argc, argv);
 
@@ -15,9 +17,15 @@
 
     fileName rootDirMaster(args.additionalArgs()[0]);
     fileName caseDirMaster(args.additionalArgs()[1]);
+    word masterRegion = polyMesh::defaultRegion;
+    args.optionReadIfPresent("masterRegion", masterRegion);
 
     fileName rootDirToAdd(args.additionalArgs()[2]);
     fileName caseDirToAdd(args.additionalArgs()[3]);
+    word addRegion = polyMesh::defaultRegion;
+    args.optionReadIfPresent("addRegion", addRegion);
 
-    Info<< "Master:      " << rootDirMaster << " " << caseDirMaster << nl
-        << "mesh to add: " << rootDirToAdd << " " << caseDirToAdd << endl;
+    Info<< "Master:      " << rootDirMaster << " " << caseDirMaster
+        << "  region " << masterRegion << nl
+        << "mesh to add: " << rootDirToAdd << " " << caseDirToAdd
+        << "  region " << addRegion << endl;
diff --git a/applications/utilities/mesh/manipulation/pointSet/pointSet.C b/applications/utilities/mesh/manipulation/pointSet/pointSet.C
index ef57d08d0779ab1da32045e5bdd55e18d52093bf..1aa70191f84cc2e1b16875a2f206ab9fcfeea870 100644
--- a/applications/utilities/mesh/manipulation/pointSet/pointSet.C
+++ b/applications/utilities/mesh/manipulation/pointSet/pointSet.C
@@ -37,37 +37,6 @@ using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Copy set
-void backup
-(
-    const polyMesh& mesh,
-    const word& fromName,
-    const topoSet& fromSet,
-    const word& toName
-)
-{
-    Info<< "Backing up " << fromName << " into " << toName << endl;
-
-    topoSet backupSet(mesh, toName, fromSet);
-
-    backupSet.write();
-}
-
-
-// Read and copy set
-void backup
-(
-    const polyMesh& mesh,
-    const word& fromName,
-    const word& toName
-)
-{
-    topoSet fromSet(mesh, fromName, IOobject::READ_IF_PRESENT);
-
-    backup(mesh, fromName, fromSet, toName);
-}
-
-
 // Main program:
 
 int main(int argc, char *argv[])
@@ -114,8 +83,6 @@ int main(int argc, char *argv[])
     {
         r = IOobject::NO_READ;
 
-        backup(mesh, setName, setName + "_old");
-
         currentSetPtr.reset
         (
             new pointSet
@@ -151,7 +118,7 @@ int main(int argc, char *argv[])
     if ((r == IOobject::MUST_READ) && (action != topoSetSource::LIST))
     {
         // currentSet has been read so can make copy.
-        backup(mesh, setName, currentSet, setName + "_old");
+        //backup(mesh, setName, currentSet, setName + "_old");
     }
 
     if (action == topoSetSource::CLEAR)
@@ -173,7 +140,16 @@ int main(int argc, char *argv[])
         forAll(topoSetSources, topoSetSourceI)
         {
             // Backup current set.
-            topoSet oldSet(mesh, currentSet.name() + "_old2", currentSet);
+            autoPtr<topoSet> oldSet
+            (
+                topoSet::New
+                (
+                    currentSet.type(),
+                    mesh,
+                    currentSet.name() + "_old2",
+                    currentSet
+                )
+            );
 
             currentSet.clear();
 
diff --git a/applications/utilities/mesh/manipulation/setSet/setSet.C b/applications/utilities/mesh/manipulation/setSet/setSet.C
index a1d7a38912a97a41878e6aef075c6c0207d729bd..3732c7e4ec1d9d17f5f55226666ba345b4a0e463 100644
--- a/applications/utilities/mesh/manipulation/setSet/setSet.C
+++ b/applications/utilities/mesh/manipulation/setSet/setSet.C
@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Description
-    Manipulate a cell/face/point set interactively.
+    Manipulate a cell/face/point/ set or zone interactively.
 
 \*---------------------------------------------------------------------------*/
 
@@ -42,6 +42,9 @@ Description
 #include "writePatch.H"
 #include "writePointSet.H"
 #include "IOobjectList.H"
+#include "cellZoneSet.H"
+#include "faceZoneSet.H"
+#include "pointZoneSet.H"
 
 #include <stdio.h>
 
@@ -82,6 +85,7 @@ Istream& selectStream(Istream* is0Ptr, Istream* is1Ptr)
 // Copy set
 void backup
 (
+    const word& setType,
     const polyMesh& mesh,
     const word& fromName,
     const topoSet& fromSet,
@@ -92,9 +96,7 @@ void backup
     {
         Pout<< "    Backing up " << fromName << " into " << toName << endl;
 
-        topoSet backupSet(mesh, toName, fromSet);
-
-        backupSet.write();
+        topoSet::New(setType, mesh, toName, fromSet)().write();
     }
 }
 
@@ -102,14 +104,21 @@ void backup
 // Read and copy set
 void backup
 (
+    const word& setType,
     const polyMesh& mesh,
     const word& fromName,
     const word& toName
 )
 {
-    topoSet fromSet(mesh, fromName, IOobject::READ_IF_PRESENT);
+    autoPtr<topoSet> fromSet = topoSet::New
+    (
+        setType,
+        mesh,
+        fromName,
+        IOobject::READ_IF_PRESENT
+    );
 
-    backup(mesh, fromName, fromSet, toName);
+    backup(setType, mesh, fromName, fromSet(), toName);
 }
 
 
@@ -241,11 +250,13 @@ void printHelp(Ostream& os)
         << "A set command should be of the following form" << endl
         << endl
         << "    cellSet|faceSet|pointSet <setName> <action> <source>"
-        << endl << endl
+        << endl
+        << endl
         << "The <action> is one of" << endl
         << "    list            - prints the contents of the set" << endl
         << "    clear           - clears the set" << endl
         << "    invert          - inverts the set" << endl
+        << "    remove          - remove the set" << endl
         << "    new <source>    - sets to set to the source set" << endl
         << "    add <source>    - adds all elements from the source set" << endl
         << "    delete <source> - deletes      ,," << endl
@@ -270,6 +281,20 @@ void printHelp(Ostream& os)
         << "    cellSet c0 add pointToCell p0 any" << endl
         << "List set:" << endl
         << "    cellSet c0 list" << endl
+        << endl
+        << "Zones can be set using zoneSets from corresponding sets:" << endl
+        << "    cellZoneSet c0Zone new setToZone c0" << endl
+        << "    faceZoneSet f0Zone new setToZone f0" << endl 
+        << endl
+        << "or if orientation is important:" << endl
+        << "    faceZoneSet f0Zone new setsToZone f0 c0" << endl 
+        << endl
+        << "ZoneSets can be manipulated using the general actions:" << endl
+        << "    list            - prints the contents of the set" << endl
+        << "    clear           - clears the set" << endl
+        << "    invert          - inverts the set (undefined orientation)"
+        << endl
+        << "    remove          - remove the set" << endl
         << endl;
 }
 
@@ -312,10 +337,126 @@ void printAllSets(const polyMesh& mesh, Ostream& os)
             os  << '\t' << set.name() << "\tsize:" << set.size() << endl;
         }
     }
+
+    const cellZoneMesh& cellZones = mesh.cellZones();
+    if (cellZones.size())
+    {
+        os  << "cellZones:" << endl;
+        forAll(cellZones, i)
+        {
+            const cellZone& zone = cellZones[i];
+            os  << '\t' << zone.name() << "\tsize:" << zone.size() << endl;
+        }
+    }
+    const faceZoneMesh& faceZones = mesh.faceZones();
+    if (faceZones.size())
+    {
+        os  << "faceZones:" << endl;
+        forAll(faceZones, i)
+        {
+            const faceZone& zone = faceZones[i];
+            os  << '\t' << zone.name() << "\tsize:" << zone.size() << endl;
+        }
+    }
+    const pointZoneMesh& pointZones = mesh.pointZones();
+    if (pointZones.size())
+    {
+        os  << "pointZones:" << endl;
+        forAll(pointZones, i)
+        {
+            const pointZone& zone = pointZones[i];
+            os  << '\t' << zone.name() << "\tsize:" << zone.size() << endl;
+        }
+    }
+
     os  << endl;
 }
 
 
+template<class ZoneType>
+void removeZone
+(
+    ZoneMesh<ZoneType, polyMesh>& zones,
+    const word& setName
+)
+{
+    label zoneID = zones.findZoneID(setName);
+
+    if (zoneID != -1)
+    {
+        Info<< "Removing zone " << setName << " at index " << zoneID << endl;
+        // Shuffle to last position
+        labelList oldToNew(zones.size());
+        label newI = 0;
+        forAll(oldToNew, i)
+        {
+            if (i != zoneID)
+            {
+                oldToNew[i] = newI++;
+            }
+        }
+        oldToNew[zoneID] = newI;
+        zones.reorder(oldToNew);
+        // Remove last element
+        zones.setSize(zones.size()-1);
+        zones.clearAddressing();
+        zones.write();
+    }
+}
+
+
+// Physically remove a set
+void removeSet
+(
+    const polyMesh& mesh,
+    const word& setType,
+    const word& setName
+)
+{
+    // Remove the file
+    IOobjectList objects
+    (
+        mesh,
+        mesh.pointsInstance(),
+        polyMesh::meshSubDir/"sets"
+    );
+
+    if (objects.found(setName))
+    {
+        // Remove file
+        fileName object = objects[setName]->objectPath();
+        Info<< "Removing file " << object << endl;
+        rm(object);
+    }
+
+    // See if zone
+    if (setType == cellZoneSet::typeName)
+    {
+        removeZone
+        (
+            const_cast<cellZoneMesh&>(mesh.cellZones()),
+            setName
+        );
+    }
+    else if (setType == faceZoneSet::typeName)
+    {
+        removeZone
+        (
+            const_cast<faceZoneMesh&>(mesh.faceZones()),
+            setName
+        );
+    }
+    else if (setType == pointZoneSet::typeName)
+    {
+        removeZone
+        (
+            const_cast<pointZoneMesh&>(mesh.pointZones()),
+            setName
+        );
+    }
+}
+
+
 // Read command and execute. Return true if ok, false otherwise.
 bool doCommand
 (
@@ -358,38 +499,29 @@ bool doCommand
 
         IOobject::readOption r;
 
-        if
+        if (action == topoSetSource::REMOVE)
+        {
+            removeSet(mesh, setType, setName);
+        }
+        else if
         (
             (action == topoSetSource::NEW)
          || (action == topoSetSource::CLEAR)
         )
         {
             r = IOobject::NO_READ;
-
-            //backup(mesh, setName, setName + "_old");
-
             currentSetPtr = topoSet::New(setType, mesh, setName, typSize);
         }
         else
         {
             r = IOobject::MUST_READ;
-
             currentSetPtr = topoSet::New(setType, mesh, setName, r);
-
             topoSet& currentSet = currentSetPtr();
-
             // Presize it according to current mesh data.
             currentSet.resize(max(currentSet.size(), typSize));
         }
 
-        if (currentSetPtr.empty())
-        {
-            Pout<< "    Cannot construct/load set "
-                << topoSet::localPath(mesh, setName) << endl;
-
-            ok = false;
-        }
-        else
+        if (currentSetPtr.valid())
         {
             topoSet& currentSet = currentSetPtr();
 
@@ -398,12 +530,6 @@ bool doCommand
                 << "  Action:" << actionName
                 << endl;
 
-            //if ((r == IOobject::MUST_READ) && (action != topoSetSource::LIST))
-            //{
-            //    // currentSet has been read so can make copy.
-            //    backup(mesh, setName, currentSet, setName + "_old");
-            //}
-
             switch (action)
             {
                 case topoSetSource::CLEAR:
@@ -437,15 +563,18 @@ bool doCommand
                         );
 
                         // Backup current set.
-                        topoSet oldSet
+                        autoPtr<topoSet> oldSet
                         (
-                            mesh,
-                            currentSet.name() + "_old2",
-                            currentSet
+                            topoSet::New
+                            (
+                                setType,
+                                mesh,
+                                currentSet.name() + "_old2",
+                                currentSet
+                            )
                         );
 
                         currentSet.clear();
-                        currentSet.resize(oldSet.size());
                         setSource().applyToSet(topoSetSource::NEW, currentSet);
 
                         // Combine new value of currentSet with old one.
@@ -547,7 +676,8 @@ enum commandStatus
 {
     QUIT,           // quit program
     INVALID,        // token is not a valid set manipulation command
-    VALID           // ,,    is a valid     ,,
+    VALIDSETCMD,    // ,,    is a valid     ,,
+    VALIDZONECMD    // ,,    is a valid     zone      ,,
 };
 
 
@@ -654,7 +784,16 @@ commandStatus parseType
      || setType == "pointSet"
     )
     {
-        return VALID;
+        return VALIDSETCMD;
+    }
+    else if
+    (
+        setType == "cellZoneSet"
+     || setType == "faceZoneSet"
+     || setType == "pointZoneSet"
+    )
+    {
+        return VALIDZONECMD;
     }
     else
     {
@@ -664,7 +803,7 @@ commandStatus parseType
             ", IStringStream&)"
         )   << "Illegal command " << setType << endl
             << "Should be one of 'help', 'list', 'time' or a set type :"
-            << " 'cellSet', 'faceSet', 'pointSet'"
+            << " 'cellSet', 'faceSet', 'pointSet', 'faceZoneSet'"
             << endl;
 
         return INVALID;
@@ -682,7 +821,7 @@ commandStatus parseAction(const word& actionName)
         {
             (void)topoSetSource::toAction(actionName);
 
-            stat = VALID;
+            stat = VALIDSETCMD;
         }
         catch (Foam::IOerror& fIOErr)
         {
@@ -724,6 +863,10 @@ int main(int argc, char *argv[])
     // Print some mesh info
     printMesh(runTime, mesh);
 
+    // Print current sets
+    printAllSets(mesh, Pout);
+
+
 
     std::ifstream* fileStreamPtr(NULL);
 
@@ -745,11 +888,10 @@ int main(int argc, char *argv[])
 #if READLINE != 0
     else if (!read_history(historyFile))
     {
-        Info<< "Successfully read history from " << historyFile << endl;
+        Pout<< "Successfully read history from " << historyFile << endl;
     }
 #endif
 
-
     Pout<< "Please type 'help', 'quit' or a set command after prompt." << endl;
 
     bool ok = true;
@@ -816,12 +958,12 @@ int main(int argc, char *argv[])
 
         IStringStream is(rawLine + ' ');
 
-        // Type: cellSet, faceSet, pointSet
+        // Type: cellSet, faceSet, pointSet, faceZoneSet
         is >> setType;
 
         stat = parseType(runTime, mesh, setType, is);
 
-        if (stat == VALID)
+        if (stat == VALIDSETCMD || stat == VALIDZONECMD)
         {
             if (is >> setName)
             {
@@ -831,14 +973,13 @@ int main(int argc, char *argv[])
                 }
             }
         }
-
         ok = true;
 
         if (stat == QUIT)
         {
             break;
         }
-        else if (stat == VALID)
+        else if (stat == VALIDSETCMD || stat == VALIDZONECMD)
         {
             ok = doCommand(mesh, setType, setName, actionName, writeVTK, is);
         }
diff --git a/applications/utilities/mesh/manipulation/splitMesh/splitMesh.C b/applications/utilities/mesh/manipulation/splitMesh/splitMesh.C
index 65e4d33237df440e9b72236179b6d617645ca7db..bc193810bd96c1c43fb4739cf734c5363a54a11f 100644
--- a/applications/utilities/mesh/manipulation/splitMesh/splitMesh.C
+++ b/applications/utilities/mesh/manipulation/splitMesh/splitMesh.C
@@ -270,6 +270,10 @@ int main(int argc, char *argv[])
     {
         mesh.setInstance(oldInstance);
     }
+    else
+    {
+        mesh.setInstance(runTime.timeName());
+    }
 
     Info<< "Writing mesh to " << runTime.timeName() << endl;
     if (!mesh.write())
diff --git a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
index a8d1ff471eba7cd99bc5ed62fe74dfa19883795a..be267f872476d9597ee360dc9ea33013451645ca 100644
--- a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
+++ b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
@@ -43,6 +43,9 @@ Description
     the largest matching region (-sloppyCellZones). This will accept any
     region that covers more than 50% of the zone. It has to be a subset
     so cannot have any cells in any other zone.
+    - useCellZonesOnly does not do a walk and uses the cellZones only. Use
+    this if you don't mind having disconnected domains in a single region.
+    This option requires all cells to be in one (and one only) region.
 
 \*---------------------------------------------------------------------------*/
 
@@ -609,7 +612,7 @@ void getInterfaceSizes
 // Create mesh for region.
 autoPtr<mapPolyMesh> createRegionMesh
 (
-    const regionSplit& cellRegion,
+    const labelList& cellRegion,
     const EdgeMap<label>& interfaceToPatch,
     const fvMesh& mesh,
     const label regionI,
@@ -766,7 +769,7 @@ autoPtr<mapPolyMesh> createRegionMesh
 void createAndWriteRegion
 (
     const fvMesh& mesh,
-    const regionSplit& cellRegion,
+    const labelList& cellRegion,
     const wordList& regionNames,
     const EdgeMap<label>& interfaceToPatch,
     const label regionI,
@@ -1005,7 +1008,8 @@ void createAndWriteRegion
 EdgeMap<label> addRegionPatches
 (
     fvMesh& mesh,
-    const regionSplit& cellRegion,
+    const labelList& cellRegion,
+    const label nCellRegions,
     const edgeList& interfaces,
     const EdgeMap<label>& interfaceSizes,
     const wordList& regionNames
@@ -1016,7 +1020,7 @@ EdgeMap<label> addRegionPatches
 
     Info<< nl << "Adding patches" << nl << endl;
 
-    EdgeMap<label> interfaceToPatch(cellRegion.nRegions());
+    EdgeMap<label> interfaceToPatch(nCellRegions);
 
     forAll(interfaces, interI)
     {
@@ -1108,13 +1112,14 @@ EdgeMap<label> addRegionPatches
 label findCorrespondingRegion
 (
     const labelList& existingZoneID,    // per cell the (unique) zoneID
-    const regionSplit& cellRegion,
+    const labelList& cellRegion,
+    const label nCellRegions,
     const label zoneI,
     const label minOverlapSize
 )
 {
     // Per region the number of cells in zoneI
-    labelList cellsInZone(cellRegion.nRegions(), 0);
+    labelList cellsInZone(nCellRegions, 0);
 
     forAll(cellRegion, cellI)
     {
@@ -1162,7 +1167,8 @@ label findCorrespondingRegion
 //(
 //    const cellZoneMesh& cellZones,
 //    const labelList& existingZoneID,    // per cell the (unique) zoneID
-//    const regionSplit& cellRegion,
+//    const labelList& cellRegion,
+//    const label nCellRegions,
 //    const label zoneI
 //)
 //{
@@ -1227,6 +1233,7 @@ label findCorrespondingRegion
 int main(int argc, char *argv[])
 {
     argList::validOptions.insert("cellZones", "");
+    argList::validOptions.insert("cellZonesOnly", "");
     argList::validOptions.insert("blockedFaces", "faceSet");
     argList::validOptions.insert("makeCellZones", "");
     argList::validOptions.insert("largestOnly", "");
@@ -1249,13 +1256,14 @@ int main(int argc, char *argv[])
             << blockedFacesName << nl << endl;
     }
 
-    bool makeCellZones   = args.optionFound("makeCellZones");
-    bool largestOnly     = args.optionFound("largestOnly");
-    bool insidePoint     = args.optionFound("insidePoint");
-    bool useCellZones    = args.optionFound("cellZones");
-    bool overwrite       = args.optionFound("overwrite");
-    bool detectOnly      = args.optionFound("detectOnly");
-    bool sloppyCellZones = args.optionFound("sloppyCellZones");
+    bool makeCellZones    = args.optionFound("makeCellZones");
+    bool largestOnly      = args.optionFound("largestOnly");
+    bool insidePoint      = args.optionFound("insidePoint");
+    bool useCellZones     = args.optionFound("cellZones");
+    bool useCellZonesOnly = args.optionFound("cellZonesOnly");
+    bool overwrite        = args.optionFound("overwrite");
+    bool detectOnly       = args.optionFound("detectOnly");
+    bool sloppyCellZones  = args.optionFound("sloppyCellZones");
 
     if (insidePoint && largestOnly)
     {
@@ -1370,13 +1378,37 @@ int main(int argc, char *argv[])
     }
 
 
-    // Do the topological walk to determine regions
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Determine per cell the region it belongs to
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // regionSplit is the labelList with the region per cell.
-    regionSplit cellRegion(mesh, blockedFace);
+    // cellRegion is the labelList with the region per cell.
+    labelList cellRegion;
+    label nCellRegions = 0;
+    if (useCellZonesOnly)
+    {
+        label unzonedCellI = findIndex(zoneID, -1);
+        if (unzonedCellI != -1)
+        {
+            FatalErrorIn(args.executable())
+                << "For the cellZonesOnly option all cells "
+                << "have to be in a cellZone." << endl
+                << "Cell " << unzonedCellI
+                << " at" << mesh.cellCentres()[unzonedCellI]
+                << " is not in a cellZone. There might be more unzoned cells."
+                << exit(FatalError);
+        }
+        cellRegion = zoneID;
+        nCellRegions = gMax(cellRegion)+1;
+    }
+    else
+    {
+        // Do a topological walk to determine regions
+        regionSplit regions(mesh, blockedFace);
+        nCellRegions = regions.nRegions();
+        cellRegion.transfer(regions);
+    }
 
-    Info<< endl << "Number of regions:" << cellRegion.nRegions() << nl << endl;
+    Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
 
 
     // Write to manual decomposition option
@@ -1429,7 +1461,7 @@ int main(int argc, char *argv[])
     // Sizes per region
     // ~~~~~~~~~~~~~~~~
 
-    labelList regionSizes(cellRegion.nRegions(), 0);
+    labelList regionSizes(nCellRegions, 0);
 
     forAll(cellRegion, cellI)
     {
@@ -1489,9 +1521,9 @@ int main(int argc, char *argv[])
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Region per zone
-    labelList regionToZone(cellRegion.nRegions(), -1);
+    labelList regionToZone(nCellRegions, -1);
     // Name of region
-    wordList regionNames(cellRegion.nRegions());
+    wordList regionNames(nCellRegions);
     // Zone to region
     labelList zoneToRegion(cellZones.size(), -1);
 
@@ -1506,6 +1538,7 @@ int main(int argc, char *argv[])
             (
                 zoneID,
                 cellRegion,
+                nCellRegions,
                 zoneI,
                 label(0.5*zoneSizes[zoneI]) // minimum overlap
             );
@@ -1532,6 +1565,7 @@ int main(int argc, char *argv[])
             (
                 zoneID,
                 cellRegion,
+                nCellRegions,
                 zoneI,
                 1               // minimum overlap
             );
@@ -1660,7 +1694,7 @@ int main(int argc, char *argv[])
     mesh.clearOut();
 
 
-    if (cellRegion.nRegions() == 1)
+    if (nCellRegions == 1)
     {
         Info<< "Only one region. Doing nothing." << endl;
     }
@@ -1671,7 +1705,7 @@ int main(int argc, char *argv[])
 
         // Check if region overlaps with existing zone. If so keep.
 
-        for (label regionI = 0; regionI < cellRegion.nRegions(); regionI++)
+        for (label regionI = 0; regionI < nCellRegions; regionI++)
         {
             label zoneI = regionToZone[regionI];
 
@@ -1759,6 +1793,7 @@ int main(int argc, char *argv[])
             (
                 mesh,
                 cellRegion,
+                nCellRegions,
                 interfaces,
                 interfaceSizes,
                 regionNames
@@ -1837,7 +1872,7 @@ int main(int argc, char *argv[])
         else
         {
             // Split all
-            for (label regionI = 0; regionI < cellRegion.nRegions(); regionI++)
+            for (label regionI = 0; regionI < nCellRegions; regionI++)
             {
                 Info<< nl
                     << "Region " << regionI << nl
diff --git a/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C b/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C
index 338c8dbcb1c6ffd1da88073613122d3d85b1a39d..f00cc8716fe0a1db3043dd93cc672e3b287a1515 100644
--- a/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C
+++ b/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C
@@ -334,7 +334,7 @@ int main(int argc, char *argv[])
     }
     else
     {
-        mesh.setInstance(oldInstance);
+        subsetter.subMesh().setInstance(oldInstance);
     }
 
     Info<< "Writing subsetted mesh and fields to time " << runTime.timeName()
diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
index e310cc758f0809cc6ae1e49fe95add7c5614bae3..d3b9bd049249ba3211935864d8438eb32a204517 100644
--- a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
+++ b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
@@ -34,6 +34,7 @@ License
 #include "Map.H"
 #include "globalMeshData.H"
 #include "DynamicList.H"
+#include "fvFieldDecomposer.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -306,12 +307,29 @@ bool domainDecomposition::writeDecomposition()
 
         forAll (curPatchSizes, patchi)
         {
+            // Get the face labels consistent with the field mapping
+            // (reuse the patch field mappers)
+            const polyPatch& meshPatch =
+                meshPatches[curBoundaryAddressing[patchi]];
+
+            fvFieldDecomposer::patchFieldDecomposer patchMapper
+            (
+                SubList<label>
+                (
+                    curFaceLabels,
+                    curPatchSizes[patchi],
+                    curPatchStarts[patchi]
+                ),
+                meshPatch.start()
+            );
+
+            // Map existing patches
             procPatches[nPatches] =
                 meshPatches[curBoundaryAddressing[patchi]].clone
                 (
                     procMesh.boundaryMesh(),
                     nPatches,
-                    curPatchSizes[patchi],
+                    patchMapper.directAddressing(),
                     curPatchStarts[patchi]
                 ).ptr();
 
diff --git a/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C b/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C
index 054c09cac5ba27beee744540cea6708bf1356dd9..8d189dfe28b9df13ea30516a8bfbf390585268c0 100644
--- a/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C
+++ b/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C
@@ -38,13 +38,14 @@ Description
 
 int main(int argc, char *argv[])
 {
+#   include "addRegionOption.H"
     timeSelector::addOptions();
     argList::validArgs.append("fieldName");
     argList::validArgs.append("patchName");
 #   include "setRootCase.H"
 #   include "createTime.H"
     instantList timeDirs = timeSelector::select0(runTime, args);
-#   include "createMesh.H"
+#   include "createNamedMesh.H"
 
     word fieldName(args.additionalArgs()[0]);
     word patchName(args.additionalArgs()[1]);
diff --git a/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryIO.C b/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryIO.C
index 896dee8d1450fe511512ee2ffebdc168d2c54089..64810f54f476d99ce3b83f8868ebca4a60b4ffab 100644
--- a/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryIO.C
+++ b/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryIO.C
@@ -81,7 +81,12 @@ bool Foam::primitiveEntry::expandVariable
     word varName = w(1, w.size()-1);
 
     // lookup the variable name in the given dictionary....
-    const entry* ePtr = dict.lookupEntryPtr(varName, true, true);
+    // Note: allow wildcards to match? For now disabled since following
+    // would expand internalField to wildcard match and not expected
+    // internalField:
+    //      internalField XXX;
+    //      boundaryField { ".*" {YYY;} movingWall {value $internalField;}
+    const entry* ePtr = dict.lookupEntryPtr(varName, true, false);
 
     // ...if defined insert its tokens into this
     if (ePtr != NULL)
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
index f471fe1b6e83f0989aba4f211f6e09e53e4cd207..ee2d966836c27c04e1e2edf97d0e6d13c329076d 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
@@ -445,6 +445,19 @@ Foam::coupledPolyPatch::coupledPolyPatch
 {}
 
 
+Foam::coupledPolyPatch::coupledPolyPatch
+(
+    const coupledPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    polyPatch(pp, bm, index, mapAddressing, newStart)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::coupledPolyPatch::~coupledPolyPatch()
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H
index de1d9cc9f98c4fe60502606c3842c1057e7c8c9c..f823148d3a48b020a9d98a57e22fa3f5a4d5c853 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H
@@ -221,6 +221,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        coupledPolyPatch
+        (
+            const coupledPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
 
     // Destructor
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/generic/genericPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/generic/genericPolyPatch.C
index 47730b06deb4c0e6880bfd26f80d114f4d581721..fc15b4ebed8a2f61c4042da0b27bf68937d96092 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/generic/genericPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/generic/genericPolyPatch.C
@@ -94,6 +94,21 @@ Foam::genericPolyPatch::genericPolyPatch
 {}
 
 
+Foam::genericPolyPatch::genericPolyPatch
+(
+    const genericPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    polyPatch(pp, bm, index, mapAddressing, newStart),
+    actualTypeName_(pp.actualTypeName_),
+    dict_(pp.dict_)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::genericPolyPatch::~genericPolyPatch()
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/generic/genericPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/generic/genericPolyPatch.H
index 53bde568734f2f4b877b89c32f431b7d9f02e9ca..ae8f7cbedb08f8f0eba9c1f09480ed1e5981b552 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/generic/genericPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/generic/genericPolyPatch.H
@@ -106,6 +106,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        genericPolyPatch
+        (
+            const genericPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -128,6 +138,22 @@ public:
             );
         }
 
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new genericPolyPatch(*this, bm, index, mapAddressing, newStart)
+            );
+        }
+
 
     // Destructor
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
index 265b61cc4bdf461bd691281c3c62bc12deb08dd6..8c6b20207913c5eac723b45455257560e4c62b92 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
@@ -826,6 +826,26 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
 {}
 
 
+Foam::cyclicPolyPatch::cyclicPolyPatch
+(
+    const cyclicPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    coupledPolyPatch(pp, bm, index, mapAddressing, newStart),
+    coupledPointsPtr_(NULL),
+    coupledEdgesPtr_(NULL),
+    featureCos_(pp.featureCos_),
+    transform_(pp.transform_),
+    rotationAxis_(pp.rotationAxis_),
+    rotationCentre_(pp.rotationCentre_),
+    separationVector_(pp.separationVector_)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::cyclicPolyPatch::~cyclicPolyPatch()
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
index 8bd2bff42f27ae16de9ff025a1cc1342468f4122..ea67348891994a77bc356d3621e4d28a4216efeb 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
@@ -234,6 +234,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        cyclicPolyPatch
+        (
+            const cyclicPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -256,6 +266,22 @@ public:
             );
         }
 
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new cyclicPolyPatch(*this, bm, index, mapAddressing, newStart)
+            );
+        }
+
 
     // Destructor
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/empty/emptyPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/empty/emptyPolyPatch.C
index 01288857583dd0c1a07574b121e1d24403a8411b..a4cf9749e4554886d53856c23d52ce1143751a55 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/empty/emptyPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/empty/emptyPolyPatch.C
@@ -87,4 +87,17 @@ Foam::emptyPolyPatch::emptyPolyPatch
 {}
 
 
+Foam::emptyPolyPatch::emptyPolyPatch
+(
+    const emptyPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    polyPatch(pp, bm, index, mapAddressing, newStart)
+{}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/empty/emptyPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/empty/emptyPolyPatch.H
index 0525e14f89a49e288e31fcc84759d615e1c5f211..286e9eeb492b623cedfe5f45144f69522cffedf8 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/empty/emptyPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/empty/emptyPolyPatch.H
@@ -92,6 +92,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        emptyPolyPatch
+        (
+            const emptyPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -113,6 +123,22 @@ public:
                 new emptyPolyPatch(*this, bm, index, newSize, newStart)
             );
         }
+
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new emptyPolyPatch(*this, bm, index, mapAddressing, newStart)
+            );
+        }
 };
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C
index 70eafc70a88d8c4ed96982789ae03e8f916fc193..6e4df9f4e5cb5b6f63c780aa83d209e2986e7ad1 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C
@@ -124,6 +124,26 @@ Foam::processorPolyPatch::processorPolyPatch
 {}
 
 
+Foam::processorPolyPatch::processorPolyPatch
+(
+    const processorPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+     const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    coupledPolyPatch(pp, bm, index, mapAddressing, newStart),
+    myProcNo_(pp.myProcNo_),
+    neighbProcNo_(pp.neighbProcNo_),
+    neighbFaceCentres_(),
+    neighbFaceAreas_(),
+    neighbFaceCellCentres_(),
+    neighbPointsPtr_(NULL),
+    neighbEdgesPtr_(NULL)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::processorPolyPatch::~processorPolyPatch()
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H
index 246f0e466ab4763a900d0ead37fce14da1d84dd7..a8948fee4376804ce4322973d00499160b94c6d1 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H
@@ -158,6 +158,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        processorPolyPatch
+        (
+            const processorPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -178,7 +188,7 @@ public:
             (
                 new processorPolyPatch
                 (
-                    refCast<const processorPolyPatch>(*this),
+                    *this,
                     bm,
                     index,
                     newSize,
@@ -187,6 +197,29 @@ public:
             );
         }
 
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new processorPolyPatch
+                (
+                    *this,
+                    bm,
+                    index,
+                    mapAddressing,
+                    newStart
+                )
+            );
+        }
+
 
     // Destructor
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/symmetry/symmetryPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/symmetry/symmetryPolyPatch.C
index 99468ea477672d9fbac6170de8381e8871cfb606..ebcd04f113afd4d6ba7860396c902133c5389269 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/symmetry/symmetryPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/symmetry/symmetryPolyPatch.C
@@ -87,4 +87,17 @@ Foam::symmetryPolyPatch::symmetryPolyPatch
 {}
 
 
+Foam::symmetryPolyPatch::symmetryPolyPatch
+(
+    const symmetryPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    polyPatch(pp, bm, index, mapAddressing, newStart)
+{}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/symmetry/symmetryPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/symmetry/symmetryPolyPatch.H
index b078076bb979fa30a625fb60e77ac64fb4c226a9..bb826d8647642f6ed3bbc46095deb7a37b0aca84 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/symmetry/symmetryPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/symmetry/symmetryPolyPatch.H
@@ -92,6 +92,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        symmetryPolyPatch
+        (
+            const symmetryPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -113,6 +123,22 @@ public:
                 new symmetryPolyPatch(*this, bm, index, newSize, newStart)
             );
         }
+
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new symmetryPolyPatch(*this, bm, index, mapAddressing, newStart)
+            );
+        }
 };
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C
index 9795b3d33289889ebb6899270b94d9347b711bf0..78a784906700817f4f0fa47a7f7461c12138d330 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C
@@ -157,4 +157,19 @@ Foam::wedgePolyPatch::wedgePolyPatch
 }
 
 
+Foam::wedgePolyPatch::wedgePolyPatch
+(
+    const wedgePolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    polyPatch(pp, bm, index, mapAddressing, newStart)
+{
+    initTransforms();
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.H
index f513037ec7ec52e025f91b6a3beff0795032b158..c4ca6add0deec239e53bf67f6e61dd799a570c99 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.H
@@ -113,6 +113,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        wedgePolyPatch
+        (
+            const wedgePolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -135,6 +145,22 @@ public:
             );
         }
 
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new wedgePolyPatch(*this, bm, index, mapAddressing, newStart)
+            );
+        }
+
 
     // Member functions
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/derived/wall/wallPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/derived/wall/wallPolyPatch.C
index 13197c86dde79d18893f7a39dde1fa6ff4131371..930759df7e64d6d3a8c5fcc21e5bad70fe29e091 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/derived/wall/wallPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/derived/wall/wallPolyPatch.C
@@ -87,4 +87,17 @@ Foam::wallPolyPatch::wallPolyPatch
 {}
 
 
+Foam::wallPolyPatch::wallPolyPatch
+(
+    const wallPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    polyPatch(pp, bm, index, mapAddressing, newStart)
+{}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/derived/wall/wallPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/derived/wall/wallPolyPatch.H
index 7950eee150d202c39169b4682dcc446230abee51..d5af8724418bf38f140aa3670513b9a4a0aa1e66 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/derived/wall/wallPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/derived/wall/wallPolyPatch.H
@@ -92,6 +92,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        wallPolyPatch
+        (
+            const wallPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -113,6 +123,22 @@ public:
                 new wallPolyPatch(*this, bm, index, newSize, newStart)
             );
         }
+
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new wallPolyPatch(*this, bm, index, mapAddressing, newStart)
+            );
+        }
 };
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
index 7c3d206db46746dcd55c6e7eec2ef2da74e98807..be0aa378a67e00fffba1b354cd2e65b885d072a1 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
@@ -167,6 +167,33 @@ Foam::polyPatch::polyPatch
 {}
 
 
+Foam::polyPatch::polyPatch
+(
+    const polyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    patchIdentifier(pp, index),
+    primitivePatch
+    (
+        faceSubList
+        (
+            bm.mesh().faces(),
+            mapAddressing.size(),
+            newStart
+        ),
+        bm.mesh().points()
+    ),
+    start_(newStart),
+    boundaryMesh_(bm),
+    faceCellsPtr_(NULL),
+    mePtr_(NULL)
+{}
+
+
 Foam::polyPatch::polyPatch(const polyPatch& p)
 :
     patchIdentifier(p),
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
index b43e4db9ff5f355672432d03e04a391aee18ac8c..f8c8fa72e72f511f4cdd29d0ead265c3ba66905b 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
@@ -199,6 +199,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        polyPatch
+        (
+            const polyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct as copy
         polyPatch(const polyPatch&);
 
@@ -224,6 +234,22 @@ public:
             );
         }
 
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new polyPatch(*this, bm, index, mapAddressing, newStart)
+            );
+        }
+
 
     // Selectors
 
diff --git a/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C b/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C
index e56730bd9eec8a55d1a3bc796e849fed0d461fd9..796e1cc9101d1703c9598dde880273340ac7b14b 100644
--- a/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C
+++ b/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C
@@ -134,16 +134,24 @@ void Foam::faceZone::calcCellLayers() const
 
         forAll (mf, faceI)
         {
+            label ownCellI = own[mf[faceI]];
+            label neiCellI =
+            (
+                zoneMesh().mesh().isInternalFace(mf[faceI])
+              ? nei[mf[faceI]]
+              : -1
+            );
+
             if (!faceFlip[faceI])
             {
                 // Face is oriented correctly, no flip needed
-                mc[faceI] = nei[mf[faceI]];
-                sc[faceI] = own[mf[faceI]];
+                mc[faceI] = neiCellI;
+                sc[faceI] = ownCellI;
             }
             else
             {
-                mc[faceI] = own[mf[faceI]];
-                sc[faceI] = nei[mf[faceI]];
+                mc[faceI] = ownCellI;
+                sc[faceI] = neiCellI;
             }
         }
         //Info << "masterCells: " << mc << endl;
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index 3d0274bc3300aafa112098a5e530b77fb817c266..c9bc6c71697613f50e3e3690eed4e19b8b17d4ef 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -2509,138 +2509,156 @@ void Foam::autoLayerDriver::addLayers
 (
     const layerParameters& layerParams,
     const dictionary& motionDict,
+    const labelList& patchIDs,
     const label nAllowableErrors,
-    motionSmoother& meshMover,
     decompositionMethod& decomposer,
     fvMeshDistribute& distributor
 )
 {
     fvMesh& mesh = meshRefiner_.mesh();
-    const indirectPrimitivePatch& pp = meshMover.patch();
-    const labelList& meshPoints = pp.meshPoints();
 
-    // Precalculate mesh edge labels for patch edges
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    labelList meshEdges(pp.nEdges());
-    forAll(pp.edges(), edgeI)
-    {
-        const edge& ppEdge = pp.edges()[edgeI];
-        label v0 = meshPoints[ppEdge[0]];
-        label v1 = meshPoints[ppEdge[1]];
-        meshEdges[edgeI] = meshTools::findEdge
+    autoPtr<indirectPrimitivePatch> pp
+    (
+        meshRefinement::makePatch
         (
-            mesh.edges(),
-            mesh.pointEdges()[v0],
-            v0,
-            v1
-        );
-    }
+            mesh,
+            patchIDs
+        )
+    );
 
-    // Displacement for all pp.localPoints.
-    vectorField patchDisp(pp.nPoints(), vector::one);
+    // Construct iterative mesh mover.
+    Info<< "Constructing mesh displacer ..." << endl;
 
-    // Number of layers for all pp.localPoints. Note: only valid if
-    // extrudeStatus = EXTRUDE.
-    labelList patchNLayers(pp.nPoints(), 0);
+    autoPtr<motionSmoother> meshMover
+    (
+        new motionSmoother
+        (
+            mesh,
+            pp(),
+            patchIDs,
+            meshRefinement::makeDisplacementField
+            (
+                pointMesh::New(mesh),
+                patchIDs
+            ),
+            motionDict
+        )
+    );
 
-    // Whether to add edge for all pp.localPoints.
-    List<extrudeMode> extrudeStatus(pp.nPoints(), EXTRUDE);
 
+    // Undistorted edge length
+    const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
+    const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
 
-    // Get number of layer per point from number of layers per patch
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    setNumLayers
-    (
-        layerParams.numLayers(),    // per patch the num layers
-        meshMover.adaptPatchIDs(),  // patches that are being moved
-        pp,                         // indirectpatch for all faces moving
+    // Point-wise extrusion data
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        patchDisp,
-        patchNLayers,
-        extrudeStatus
-    );
+    // Displacement for all pp.localPoints.
+    vectorField patchDisp(pp().nPoints(), vector::one);
 
-    // Disable extrusion on non-manifold points
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Number of layers for all pp.localPoints. Note: only valid if
+    // extrudeStatus = EXTRUDE.
+    labelList patchNLayers(pp().nPoints(), 0);
 
-    handleNonManifolds
-    (
-        pp,
-        meshEdges,
+    // Whether to add edge for all pp.localPoints.
+    List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
 
-        patchDisp,
-        patchNLayers,
-        extrudeStatus
-    );
+    {
+        // Get number of layer per point from number of layers per patch
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Disable extrusion on feature angles
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        setNumLayers
+        (
+            layerParams.numLayers(),    // per patch the num layers
+            meshMover().adaptPatchIDs(),// patches that are being moved
+            pp,                         // indirectpatch for all faces moving
 
-    handleFeatureAngle
-    (
-        pp,
-        meshEdges,
-        layerParams.featureAngle()*constant::math::pi/180.0,
+            patchDisp,
+            patchNLayers,
+            extrudeStatus
+        );
 
-        patchDisp,
-        patchNLayers,
-        extrudeStatus
-    );
+        // Precalculate mesh edge labels for patch edges
+        labelList meshEdges(pp().meshEdges(mesh.edges(), mesh.pointEdges()));
 
-    // Disable extrusion on warped faces
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // Disable extrusion on non-manifold points
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Undistorted edge length
-    const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
-    const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
+        handleNonManifolds
+        (
+            pp,
+            meshEdges,
 
-    handleWarpedFaces
-    (
-        pp,
-        layerParams.maxFaceThicknessRatio(),
-        edge0Len,
-        cellLevel,
+            patchDisp,
+            patchNLayers,
+            extrudeStatus
+        );
 
-        patchDisp,
-        patchNLayers,
-        extrudeStatus
-    );
+        // Disable extrusion on feature angles
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    //// Disable extrusion on cells with multiple patch faces
-    //// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    //
-    //handleMultiplePatchFaces
-    //(
-    //    pp,
-    //
-    //    patchDisp,
-    //    patchNLayers,
-    //    extrudeStatus
-    //);
+        handleFeatureAngle
+        (
+            pp,
+            meshEdges,
+            layerParams.featureAngle()*constant::math::pi/180.0,
+
+            patchDisp,
+            patchNLayers,
+            extrudeStatus
+        );
 
+        // Disable extrusion on warped faces
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Grow out region of non-extrusion
-    for (label i = 0; i < layerParams.nGrow(); i++)
-    {
-        growNoExtrusion
+        handleWarpedFaces
         (
             pp,
+            layerParams.maxFaceThicknessRatio(),
+            edge0Len,
+            cellLevel,
+
             patchDisp,
             patchNLayers,
             extrudeStatus
         );
+
+        //// Disable extrusion on cells with multiple patch faces
+        //// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        //
+        //handleMultiplePatchFaces
+        //(
+        //    pp,
+        //
+        //    patchDisp,
+        //    patchNLayers,
+        //    extrudeStatus
+        //);
+
+
+        // Grow out region of non-extrusion
+        for (label i = 0; i < layerParams.nGrow(); i++)
+        {
+            growNoExtrusion
+            (
+                pp,
+                patchDisp,
+                patchNLayers,
+                extrudeStatus
+            );
+        }
     }
 
+
     // Determine (wanted) point-wise layer thickness and expansion ratio
-    scalarField thickness(pp.nPoints());
-    scalarField minThickness(pp.nPoints());
-    scalarField expansionRatio(pp.nPoints());
+    scalarField thickness(pp().nPoints());
+    scalarField minThickness(pp().nPoints());
+    scalarField expansionRatio(pp().nPoints());
     calculateLayerThickness
     (
         pp,
-        meshMover.adaptPatchIDs(),
+        meshMover().adaptPatchIDs(),
         layerParams.expansionRatio(),
 
         layerParams.relativeSizes(),        // thickness relative to cellsize?
@@ -2663,11 +2681,11 @@ void Foam::autoLayerDriver::addLayers
 
         // Find maximum length of a patch name, for a nicer output
         label maxPatchNameLen = 0;
-        forAll(meshMover.adaptPatchIDs(), i)
+        forAll(meshMover().adaptPatchIDs(), i)
         {
-            label patchI = meshMover.adaptPatchIDs()[i];
+            label patchI = meshMover().adaptPatchIDs()[i];
             word patchName = patches[patchI].name();
-            maxPatchNameLen = max(maxPatchNameLen,label(patchName.size()));
+            maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
         }
 
         Info<< nl 
@@ -2678,9 +2696,9 @@ void Foam::autoLayerDriver::addLayers
             << setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
             << setw(0) << " -----    ------ --------- -------" << endl;
 
-        forAll(meshMover.adaptPatchIDs(), i)
+        forAll(meshMover().adaptPatchIDs(), i)
         {
-            label patchI = meshMover.adaptPatchIDs()[i];
+            label patchI = meshMover().adaptPatchIDs()[i];
 
             const labelList& meshPoints = patches[patchI].meshPoints();
 
@@ -2691,7 +2709,7 @@ void Foam::autoLayerDriver::addLayers
 
             forAll(meshPoints, patchPointI)
             {
-                label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
+                label ppPointI = pp().meshPointMap()[meshPoints[patchPointI]];
 
                 //maxThickness = max(maxThickness, thickness[ppPointI]);
                 //minThickness = min(minThickness, thickness[ppPointI]);
@@ -2761,7 +2779,7 @@ void Foam::autoLayerDriver::addLayers
             IOobject::NO_WRITE,
             false
         ),
-        meshMover.pMesh(),
+        meshMover().pMesh(),
         dimensionedScalar("pointMedialDist", dimless, 0.0)
     );
 
@@ -2776,7 +2794,7 @@ void Foam::autoLayerDriver::addLayers
             IOobject::NO_WRITE,
             false
         ),
-        meshMover.pMesh(),
+        meshMover().pMesh(),
         dimensionedVector("dispVec", dimless, vector::zero)
     );
 
@@ -2791,7 +2809,7 @@ void Foam::autoLayerDriver::addLayers
             IOobject::NO_WRITE,
             false
         ),
-        meshMover.pMesh(),
+        meshMover().pMesh(),
         dimensionedScalar("medialRatio", dimless, 0.0)
     );
 
@@ -2871,7 +2889,7 @@ void Foam::autoLayerDriver::addLayers
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
         {
-            pointField oldPatchPos(pp.localPoints());
+            pointField oldPatchPos(pp().localPoints());
 
             //// Laplacian displacement shrinking.
             //shrinkMeshDistance
@@ -2886,7 +2904,7 @@ void Foam::autoLayerDriver::addLayers
             // Medial axis based shrinking
             shrinkMeshMedialDistance
             (
-                meshMover,
+                meshMover(),
                 meshQualityDict,
 
                 layerParams.nSmoothThickness(),
@@ -2908,7 +2926,7 @@ void Foam::autoLayerDriver::addLayers
             );
 
             // Update patchDisp (since not all might have been honoured)
-            patchDisp = oldPatchPos - pp.localPoints();
+            patchDisp = oldPatchPos - pp().localPoints();
         }
 
         // Truncate displacements that are too small (this will do internal
@@ -2916,7 +2934,7 @@ void Foam::autoLayerDriver::addLayers
         faceSet dummySet(mesh, "wrongPatchFaces", 0);
         truncateDisplacement
         (
-            meshMover,
+            meshMover(),
             minThickness,
             dummySet,
             patchDisp,
@@ -2931,7 +2949,7 @@ void Foam::autoLayerDriver::addLayers
             dumpDisplacement
             (
                 mesh.time().path()/"layer",
-                pp,
+                pp(),
                 patchDisp,
                 extrudeStatus
             );
@@ -2955,11 +2973,11 @@ void Foam::autoLayerDriver::addLayers
         // Determine per point/per face number of layers to extrude. Also
         // handles the slow termination of layers when going switching layers
 
-        labelList nPatchPointLayers(pp.nPoints(),-1);
-        labelList nPatchFaceLayers(pp.localFaces().size(),-1);
+        labelList nPatchPointLayers(pp().nPoints(),-1);
+        labelList nPatchFaceLayers(pp().localFaces().size(),-1);
         setupLayerInfoTruncation
         (
-            meshMover,
+            meshMover(),
             patchNLayers,
             extrudeStatus,
             layerParams.nBufferCellsNoExtrude(),
@@ -2998,7 +3016,7 @@ void Foam::autoLayerDriver::addLayers
         addLayer.setRefinement
         (
             invExpansionRatio,
-            pp,
+            pp(),
             labelList(0),       // exposed patchIDs, not used for adding layers
             nPatchFaceLayers,   // layers per face
             nPatchPointLayers,  // layers per point
@@ -3050,8 +3068,8 @@ void Foam::autoLayerDriver::addLayers
         addLayer.updateMesh
         (
             map,
-            identity(pp.size()),
-            identity(pp.nPoints())
+            identity(pp().size()),
+            identity(pp().nPoints())
         );
 
         // Collect layer faces and cells for outside loop.
@@ -3098,7 +3116,7 @@ void Foam::autoLayerDriver::addLayers
         (
             addLayer,
             meshQualityDict,
-            pp,
+            pp(),
             newMesh,
 
             patchDisp,
@@ -3107,7 +3125,7 @@ void Foam::autoLayerDriver::addLayers
         );
 
         Info<< "Extruding " << countExtrusion(pp, extrudeStatus)
-            << " out of " << returnReduce(pp.size(), sumOp<label>())
+            << " out of " << returnReduce(pp().size(), sumOp<label>())
             << " faces. Removed extrusion at " << nTotChanged << " faces."
             << endl;
 
@@ -3117,8 +3135,8 @@ void Foam::autoLayerDriver::addLayers
         }
 
         // Reset mesh points and start again
-        meshMover.movePoints(oldPoints);
-        meshMover.correct();
+        meshMover().movePoints(oldPoints);
+        meshMover().correct();
 
         Info<< endl;
     }
@@ -3173,6 +3191,7 @@ void Foam::autoLayerDriver::addLayers
         (
             false,
             false,
+            scalarField(mesh.nCells(), 1.0),
             decomposer,
             distributor
         );
@@ -3211,7 +3230,7 @@ void Foam::autoLayerDriver::doLayers
     fvMeshDistribute& distributor
 )
 {
-    fvMesh& mesh = meshRefiner_.mesh();
+    const fvMesh& mesh = meshRefiner_.mesh();
 
     Info<< nl
         << "Shrinking and layer addition phase" << nl
@@ -3245,59 +3264,80 @@ void Foam::autoLayerDriver::doLayers
     }
     else
     {
-        autoPtr<indirectPrimitivePatch> ppPtr
+        // Check that outside of mesh is not multiply connected.
+        checkMeshManifold();
+
+        // Check initial mesh
+        Info<< "Checking initial mesh ..." << endl;
+        labelHashSet wrongFaces(mesh.nFaces()/100);
+        motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
+        const label nInitErrors = returnReduce
         (
-            meshRefinement::makePatch
-            (
-                mesh,
-                patchIDs
-            )
+            wrongFaces.size(),
+            sumOp<label>()
         );
-        indirectPrimitivePatch& pp = ppPtr();
 
-        // Construct iterative mesh mover.
-        Info<< "Constructing mesh displacer ..." << endl;
-
-        {
-            const pointMesh& pMesh = pointMesh::New(mesh);
-
-            motionSmoother meshMover
-            (
-                mesh,
-                pp,
-                patchIDs,
-                meshRefinement::makeDisplacementField(pMesh, patchIDs),
-                motionDict
-            );
-
-            // Check that outside of mesh is not multiply connected.
-            checkMeshManifold();
+        Info<< "Detected " << nInitErrors << " illegal faces"
+            << " (concave, zero area or negative cell pyramid volume)"
+            << endl;
 
-            // Check initial mesh
-            Info<< "Checking initial mesh ..." << endl;
-            labelHashSet wrongFaces(mesh.nFaces()/100);
-            motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
-            const label nInitErrors = returnReduce
-            (
-                wrongFaces.size(),
-                sumOp<label>()
-            );
 
-            Info<< "Detected " << nInitErrors << " illegal faces"
-                << " (concave, zero area or negative cell pyramid volume)"
+        // Balance
+        if (Pstream::parRun())
+        {
+            Info<< nl
+                << "Doing initial balancing" << nl
+                << "-----------------------" << nl
                 << endl;
-
-            // Do all topo changes
-            addLayers
+        
+            scalarField cellWeights(mesh.nCells(), 1);
+            forAll(numLayers, patchI)
+            {
+                if (numLayers[patchI] > 0)
+                {
+                    const polyPatch& pp = mesh.boundaryMesh()[patchI];
+                    forAll(pp.faceCells(), i)
+                    {
+                        cellWeights[pp.faceCells()[i]] += numLayers[patchI];
+                    }
+                }
+            }
+        
+            // Balance mesh (and meshRefinement). No restriction on face zones
+            // and baffles.
+            autoPtr<mapDistributePolyMesh> map = meshRefiner_.balance
             (
-                layerParams,
-                motionDict,
-                nInitErrors,
-                meshMover,
+                false,
+                false,
+                cellWeights,
                 decomposer,
                 distributor
             );
+        
+            //{
+            //    globalIndex globalCells(mesh.nCells());
+            //
+            //    Info<< "** Distribution after balancing:" << endl;
+            //    for (label procI = 0; procI < Pstream::nProcs(); procI++)
+            //    {
+            //        Info<< "    " << procI << '\t'
+            //            << globalCells.localSize(procI) << endl;
+            //    }
+            //    Info<< endl;
+            //}
         }
+
+
+        // Do all topo changes
+        addLayers
+        (
+            layerParams,
+            motionDict,
+            patchIDs,
+            nInitErrors,
+            decomposer,
+            distributor
+        );
     }
 }
 
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
index 12569570069144b603074002ea3ad7af8c3bb91f..512766c53937a9739df92997c99dfcc5cce41677 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
@@ -523,8 +523,8 @@ public:
             (
                 const layerParameters& layerParams,
                 const dictionary& motionDict,
+                const labelList& patchIDs,
                 const label nAllowableErrors,
-                motionSmoother& meshMover,
                 decompositionMethod& decomposer,
                 fvMeshDistribute& distributor
             );
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index 3ea198e6c4ebc2b216c6e2471495f27d504d1049..41fe3853142d559b3a1245c3c9cf21d59e1739e5 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -35,7 +35,6 @@ License
 #include "refinementSurfaces.H"
 #include "shellSurfaces.H"
 #include "mapDistributePolyMesh.H"
-#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -210,7 +209,8 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
                 "feature refinement iteration " + name(iter),
                 decomposer_,
                 distributor_,
-                cellsToRefine
+                cellsToRefine,
+                refineParams.maxLoadUnbalance()
             );
         }
     }
@@ -310,7 +310,8 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
             "surface refinement iteration " + name(iter),
             decomposer_,
             distributor_,
-            cellsToRefine
+            cellsToRefine,
+            refineParams.maxLoadUnbalance()
         );
     }
     return iter;
@@ -492,7 +493,8 @@ Foam::label Foam::autoRefineDriver::shellRefine
             "shell refinement iteration " + name(iter),
             decomposer_,
             distributor_,
-            cellsToRefine
+            cellsToRefine,
+            refineParams.maxLoadUnbalance()
         );
     }
     meshRefiner_.userFaceData().clear();
@@ -776,11 +778,14 @@ void Foam::autoRefineDriver::doRefine
             const_cast<Time&>(mesh.time())++;
         }
 
-        // Do final balancing. Keep zoned faces on one processor.
+        // Do final balancing. Keep zoned faces on one processor since the
+        // snap phase will convert them to baffles and this only works for
+        // internal faces.
         meshRefiner_.balance
         (
             true,
             false,
+            scalarField(mesh.nCells(), 1), // dummy weights
             decomposer_,
             distributor_
         );
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
index ce3b82bfc6c91e80b5ac269a0194330f31c6d72e..8c79ab9c6c4b30ab37e6e8ff134bf67f3cae2f46 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
@@ -43,7 +43,8 @@ Foam::refinementParameters::refinementParameters
     minRefineCells_(readLabel(dict.lookup("minimumRefine"))),
     curvature_(readScalar(dict.lookup("curvature"))),
     nBufferLayers_(readLabel(dict.lookup("nBufferLayers"))),
-    keepPoints_(dict.lookup("keepPoints"))
+    keepPoints_(dict.lookup("keepPoints")),
+    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0))
 {}
 
 
@@ -53,7 +54,8 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
     maxLocalCells_(readLabel(dict.lookup("maxLocalCells"))),
     minRefineCells_(readLabel(dict.lookup("minRefinementCells"))),
     nBufferLayers_(readLabel(dict.lookup("nCellsBetweenLevels"))),
-    keepPoints_(pointField(1, dict.lookup("locationInMesh")))
+    keepPoints_(pointField(1, dict.lookup("locationInMesh"))),
+    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0))
 {
     scalar featAngle(readScalar(dict.lookup("resolveFeatureAngle")));
 
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H b/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H
index 1e12dd30a943c74f5087617a83c3a90b4c1621c1..ae09e3fc0e7e9473ff86805a4bc525d3d0f3c43b 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H
@@ -73,6 +73,9 @@ class refinementParameters
         //- Areas to keep
         const pointField keepPoints_;
 
+        //- Allowed load unbalance
+        scalar maxLoadUnbalance_;
+
 
     // Private Member Functions
 
@@ -134,6 +137,12 @@ public:
                 return keepPoints_;
             }
 
+            //- Allowed load unbalance
+            scalar maxLoadUnbalance() const
+            {
+                return maxLoadUnbalance_;
+            }
+
 
         // Other
 
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index 7e90beae6f6699b8231773e40e3c76f2f9a06f5d..1064aa3682347095f87f4b2fa455cf3ac60fcb51 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -552,13 +552,16 @@ void Foam::meshRefinement::calcLocalRegions
     const globalIndex& globalCells,
     const labelList& globalRegion,
     const Map<label>& coupledRegionToMaster,
+    const scalarField& cellWeights,
 
     Map<label>& globalToLocalRegion,
-    pointField& localPoints
+    pointField& localPoints,
+    scalarField& localWeights
 ) const
 {
     globalToLocalRegion.resize(globalRegion.size());
     DynamicList<point> localCc(globalRegion.size()/2);
+    DynamicList<scalar> localWts(globalRegion.size()/2);
 
     forAll(globalRegion, cellI)
     {
@@ -573,6 +576,7 @@ void Foam::meshRefinement::calcLocalRegions
                 // I am master. Allocate region for me.
                 globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
                 localCc.append(mesh_.cellCentres()[cellI]);
+                localWts.append(cellWeights[cellI]);
             }
         }
         else
@@ -581,11 +585,13 @@ void Foam::meshRefinement::calcLocalRegions
             if (globalToLocalRegion.insert(globalRegion[cellI], localCc.size()))
             {
                 localCc.append(mesh_.cellCentres()[cellI]);
+                localWts.append(cellWeights[cellI]);
             }
         }
     }
 
     localPoints.transfer(localCc);
+    localWeights.transfer(localWts);
 
     if (localPoints.size() != globalToLocalRegion.size())
     {
@@ -924,6 +930,7 @@ Foam::label Foam::meshRefinement::countHits() const
 // Determine distribution to move connected regions onto one processor.
 Foam::labelList Foam::meshRefinement::decomposeCombineRegions
 (
+    const scalarField& cellWeights,
     const boolList& blockedFace,
     const List<labelPair>& explicitConnections,
     decompositionMethod& decomposer
@@ -965,14 +972,17 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
 
     Map<label> globalToLocalRegion;
     pointField localPoints;
+    scalarField localWeights;
     calcLocalRegions
     (
         globalCells,
         globalRegion,
         coupledRegionToMaster,
+        cellWeights,
 
         globalToLocalRegion,
-        localPoints
+        localPoints,
+        localWeights
     );
 
 
@@ -984,7 +994,7 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
 
     if (isA<geomDecomp>(decomposer))
     {
-        regionDistribution = decomposer.decompose(localPoints);
+        regionDistribution = decomposer.decompose(localPoints, localWeights);
     }
     else
     {
@@ -998,7 +1008,12 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
             regionRegions
         );
 
-        regionDistribution = decomposer.decompose(regionRegions, localPoints);
+        regionDistribution = decomposer.decompose
+        (
+            regionRegions,
+            localPoints,
+            localWeights
+        );
     }
 
 
@@ -1058,6 +1073,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
 (
     const bool keepZoneFaces,
     const bool keepBaffles,
+    const scalarField& cellWeights,
     decompositionMethod& decomposer,
     fvMeshDistribute& distributor
 )
@@ -1151,6 +1167,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
 
                 distribution = decomposeCombineRegions
                 (
+                    cellWeights,
                     blockedFace,
                     couples,
                     decomposer
@@ -1170,13 +1187,21 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
             else
             {
                 // Normal decomposition
-                distribution = decomposer.decompose(mesh_.cellCentres());
+                distribution = decomposer.decompose
+                (
+                    mesh_.cellCentres(),
+                    cellWeights
+                );
             }
         }
         else
         {
             // Normal decomposition
-            distribution = decomposer.decompose(mesh_.cellCentres());
+            distribution = decomposer.decompose
+            (
+                mesh_.cellCentres(),
+                cellWeights
+            );
         }
 
         if (debug)
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 93ce3ceeb8fb3043d0b8435ef22209d6ee073711..5b603ddbac9660d1e92a3eda8968686f7700247a 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -214,9 +214,11 @@ private:
                 const globalIndex& globalCells,
                 const labelList& globalRegion,
                 const Map<label>& coupledRegionToMaster,
+                const scalarField& cellWeights,
 
                 Map<label>& globalToLocalRegion,
-                pointField& localPoints
+                pointField& localPoints,
+                scalarField& localWeights
             ) const;
 
             //- Convert region into global index.
@@ -574,6 +576,7 @@ public:
             //  (e.g. to keep baffles together)
             labelList decomposeCombineRegions
             (
+                const scalarField& cellWeights,
                 const boolList& blockedFace,
                 const List<labelPair>& explicitConnections,
                 decompositionMethod&
@@ -587,6 +590,7 @@ public:
             (
                 const bool keepZoneFaces,
                 const bool keepBaffles,
+                const scalarField& cellWeights,
                 decompositionMethod& decomposer,
                 fvMeshDistribute& distributor
             );
@@ -648,7 +652,8 @@ public:
                 const string& msg,
                 decompositionMethod& decomposer,
                 fvMeshDistribute& distributor,
-                const labelList& cellsToRefine
+                const labelList& cellsToRefine,
+                const scalar maxLoadUnbalance
             );
 
 
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index e35607a5980172c0d90377eae1f38efeeca6614f..c77a634d301c06ab803298c0a4a06ee4d18ebbaa 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -37,6 +37,7 @@ License
 #include "mapDistributePolyMesh.H"
 #include "featureEdgeMesh.H"
 #include "Cloud.H"
+//#include "globalIndex.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -1244,40 +1245,122 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::refine
 }
 
 
-// Do refinement of consistent set of cells followed by truncation and
-// load balancing.
+//// Do refinement of consistent set of cells followed by truncation and
+//// load balancing.
+//Foam::autoPtr<Foam::mapDistributePolyMesh>
+//Foam::meshRefinement::refineAndBalance
+//(
+//    const string& msg,
+//    decompositionMethod& decomposer,
+//    fvMeshDistribute& distributor,
+//    const labelList& cellsToRefine
+//)
+//{
+//    // Do all refinement
+//    refine(cellsToRefine);
+//
+//    if (debug)
+//    {
+//        Pout<< "Writing refined but unbalanced " << msg
+//            << " mesh to time " << timeName() << endl;
+//        write
+//        (
+//            debug,
+//            mesh_.time().path()
+//           /timeName()
+//        );
+//        Pout<< "Dumped debug data in = "
+//            << mesh_.time().cpuTimeIncrement() << " s" << endl;
+//
+//        // test all is still synced across proc patches
+//        checkData();
+//    }
+//
+//    Info<< "Refined mesh in = "
+//        << mesh_.time().cpuTimeIncrement() << " s" << endl;
+//    printMeshInfo(debug, "After refinement " + msg);
+//
+//
+//    // Load balancing
+//    // ~~~~~~~~~~~~~~
+//
+//    autoPtr<mapDistributePolyMesh> distMap;
+//
+//    if (Pstream::nProcs() > 1)
+//    {
+//        scalarField cellWeights(mesh_.nCells(), 1);
+//
+//        distMap = balance
+//        (
+//            false,  //keepZoneFaces
+//            false,  //keepBaffles
+//            cellWeights,
+//            decomposer,
+//            distributor
+//        );
+//
+//        Info<< "Balanced mesh in = "
+//            << mesh_.time().cpuTimeIncrement() << " s" << endl;
+//
+//        printMeshInfo(debug, "After balancing " + msg);
+//
+//
+//        if (debug)
+//        {
+//            Pout<< "Writing balanced " << msg
+//                << " mesh to time " << timeName() << endl;
+//            write
+//            (
+//                debug,
+//                mesh_.time().path()/timeName()
+//            );
+//            Pout<< "Dumped debug data in = "
+//                << mesh_.time().cpuTimeIncrement() << " s" << endl;
+//
+//            // test all is still synced across proc patches
+//            checkData();
+//        }
+//    }
+//
+//    return distMap;
+//}
+
+
+// Do load balancing followed by refinement of consistent set of cells.
 Foam::autoPtr<Foam::mapDistributePolyMesh>
 Foam::meshRefinement::refineAndBalance
 (
     const string& msg,
     decompositionMethod& decomposer,
     fvMeshDistribute& distributor,
-    const labelList& cellsToRefine
+    const labelList& initCellsToRefine,
+    const scalar maxLoadUnbalance
 )
 {
-    // Do all refinement
-    refine(cellsToRefine);
-
-    if (debug)
-    {
-        Pout<< "Writing refined but unbalanced " << msg
-            << " mesh to time " << timeName() << endl;
-        write
-        (
-            debug,
-            mesh_.time().path()
-           /timeName()
-        );
-        Pout<< "Dumped debug data in = "
-            << mesh_.time().cpuTimeIncrement() << " s" << endl;
-
-        // test all is still synced across proc patches
-        checkData();
-    }
+    labelList cellsToRefine(initCellsToRefine);
 
-    Info<< "Refined mesh in = "
-        << mesh_.time().cpuTimeIncrement() << " s" << endl;
-    printMeshInfo(debug, "After refinement " + msg);
+    //{
+    //    globalIndex globalCells(mesh_.nCells());
+    //
+    //    Info<< "** Distribution before balancing/refining:" << endl;
+    //    for (label procI = 0; procI < Pstream::nProcs(); procI++)
+    //    {
+    //        Info<< "    " << procI << '\t'
+    //            << globalCells.localSize(procI) << endl;
+    //    }
+    //    Info<< endl;
+    //}
+    //{
+    //    globalIndex globalCells(cellsToRefine.size());
+    //
+    //    Info<< "** Cells to be refined:" << endl;
+    //    for (label procI = 0; procI < Pstream::nProcs(); procI++)
+    //    {
+    //        Info<< "    " << procI << '\t'
+    //            << globalCells.localSize(procI) << endl;
+    //    }
+    //    Info<< endl;
+    //}
 
 
     // Load balancing
@@ -1287,19 +1370,61 @@ Foam::meshRefinement::refineAndBalance
 
     if (Pstream::nProcs() > 1)
     {
-        distMap = balance
+        // First check if we need to balance at all. Precalculate number of
+        // cells after refinement and see what maximum difference is.
+        scalar nNewCells = scalar(mesh_.nCells() + 7*cellsToRefine.size());
+        scalar nIdealNewCells =
+            returnReduce(nNewCells, sumOp<scalar>())
+          / Pstream::nProcs();
+        scalar unbalance = returnReduce
         (
-            false,  //keepZoneFaces
-            false,  //keepBaffles
-            decomposer,
-            distributor
+            mag(1.0-nNewCells/nIdealNewCells),
+            maxOp<scalar>()
         );
 
-        Info<< "Balanced mesh in = "
-            << mesh_.time().cpuTimeIncrement() << " s" << endl;
+        if (unbalance <= maxLoadUnbalance)
+        {
+            Info<< "Skipping balancing since max unbalance " << unbalance
+                << " in = "
+                << mesh_.time().cpuTimeIncrement() << " s" << endl;
+        }
+        else
+        {
+            scalarField cellWeights(mesh_.nCells(), 1);
+            forAll(cellsToRefine, i)
+            {
+                cellWeights[cellsToRefine[i]] += 7;
+            }
 
-        printMeshInfo(debug, "After balancing " + msg);
+            distMap = balance
+            (
+                false,  //keepZoneFaces
+                false,  //keepBaffles
+                cellWeights,
+                decomposer,
+                distributor
+            );
 
+            // Update cells to refine
+            distMap().distributeCellIndices(cellsToRefine);
+
+            Info<< "Balanced mesh in = "
+                << mesh_.time().cpuTimeIncrement() << " s" << endl;
+        }
+
+        //{
+        //    globalIndex globalCells(mesh_.nCells());
+        //
+        //    Info<< "** Distribution after balancing:" << endl;
+        //    for (label procI = 0; procI < Pstream::nProcs(); procI++)
+        //    {
+        //        Info<< "    " << procI << '\t'
+        //            << globalCells.localSize(procI) << endl;
+        //    }
+        //    Info<< endl;
+        //}
+
+        printMeshInfo(debug, "After balancing " + msg);
 
         if (debug)
         {
@@ -1318,6 +1443,46 @@ Foam::meshRefinement::refineAndBalance
         }
     }
 
+
+    // Refinement
+    // ~~~~~~~~~~
+
+    refine(cellsToRefine);
+
+    if (debug)
+    {
+        Pout<< "Writing refined " << msg
+            << " mesh to time " << timeName() << endl;
+        write
+        (
+            debug,
+            mesh_.time().path()
+           /timeName()
+        );
+        Pout<< "Dumped debug data in = "
+            << mesh_.time().cpuTimeIncrement() << " s" << endl;
+
+        // test all is still synced across proc patches
+        checkData();
+    }
+
+    Info<< "Refined mesh in = "
+        << mesh_.time().cpuTimeIncrement() << " s" << endl;
+
+    //{
+    //    globalIndex globalCells(mesh_.nCells());
+    //
+    //    Info<< "** After refinement distribution:" << endl;
+    //    for (label procI = 0; procI < Pstream::nProcs(); procI++)
+    //    {
+    //        Info<< "    " << procI << '\t'
+    //            << globalCells.localSize(procI) << endl;
+    //    }
+    //    Info<< endl;
+    //}
+
+    printMeshInfo(debug, "After refinement " + msg);
+
     return distMap;
 }
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/totalTemperature/totalTemperatureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/totalTemperature/totalTemperatureFvPatchScalarField.C
index ca77c9aed47381c493bd2ee8b8a3e84717a1b502..8c29e5a0e591bbd7bfb562623c836b2ccd9eaeed 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/totalTemperature/totalTemperatureFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/totalTemperature/totalTemperatureFvPatchScalarField.C
@@ -188,7 +188,7 @@ void Foam::totalTemperatureFvPatchScalarField::write(Ostream& os) const
     {
         os.writeKeyword("phi") << phiName_ << token::END_STATEMENT << nl;
     }
-    if (phiName_ != "psi")
+    if (psiName_ != "psi")
     {
         os.writeKeyword("psi") << psiName_ << token::END_STATEMENT << nl;
     }
diff --git a/src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C b/src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C
index bbafb189faa97eba4faea0aa5e69ad1d163e70fd..9c0cf511dce6faa06860e2d7fd85501439cf4c04 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C
+++ b/src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C
@@ -329,7 +329,7 @@ surfaceDisplacementPointPatchVectorField
     surfacesDict_(dict.subDict("geometry")),
     projectMode_(projectModeNames_.read(dict.lookup("projectMode"))),
     projectDir_(dict.lookup("projectDirection")),
-    wedgePlane_(readLabel(dict.lookup("wedgePlane"))),
+    wedgePlane_(dict.lookupOrDefault(dict.lookup("wedgePlane"), -1)),
     frozenPointsZone_(dict.lookupOrDefault("frozenPointsZone", word::null))
 {
     if (velocity_.x() < 0 || velocity_.y() < 0 || velocity_.z() < 0)
diff --git a/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C b/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C
index 8f37cd91abadbe0bc26a3f7a49985a961902721c..4944c82ff1c5cf8f94513bb1010780b7a27182bb 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C
+++ b/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C
@@ -326,7 +326,7 @@ surfaceSlipDisplacementPointPatchVectorField
     surfacesDict_(dict.subDict("geometry")),
     projectMode_(projectModeNames_.read(dict.lookup("projectMode"))),
     projectDir_(dict.lookup("projectDirection")),
-    wedgePlane_(readLabel(dict.lookup("wedgePlane"))),
+    wedgePlane_(dict.lookupOrDefault("wedgePlane", -1)),
     frozenPointsZone_(dict.lookupOrDefault("frozenPointsZone", word::null))
 {}
 
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index 3f8eda4f307215936632f3a4eb1d709e3fa76e71..46d47a3e56b6aee9439063c5980b804a8d7c408c 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -72,6 +72,9 @@ $(topoSets)/cellSet.C
 $(topoSets)/topoSet.C
 $(topoSets)/faceSet.C
 $(topoSets)/pointSet.C
+$(topoSets)/cellZoneSet.C
+$(topoSets)/faceZoneSet.C
+$(topoSets)/pointZoneSet.C
 
 sets/topoSetSource/topoSetSource.C
 
@@ -114,6 +117,18 @@ $(pointSources)/surfaceToPoint/surfaceToPoint.C
 $(pointSources)/zoneToPoint/zoneToPoint.C
 $(pointSources)/nearestToPoint/nearestToPoint.C
 
+faceZoneSources = sets/faceZoneSources
+$(faceZoneSources)/faceZoneToFaceZone/faceZoneToFaceZone.C
+$(faceZoneSources)/setsToFaceZone/setsToFaceZone.C
+$(faceZoneSources)/setToFaceZone/setToFaceZone.C
+
+cellZoneSources = sets/cellZoneSources
+$(cellZoneSources)/setToCellZone/setToCellZone.C
+
+pointZoneSources = sets/pointZoneSources
+$(pointZoneSources)/setToPointZone/setToPointZone.C
+
+
 surfaceSets/surfaceSets.C
 
 triSurface/orientedSurface/orientedSurface.C
diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C
index 20dc01bc60aa4394bc07a7e3bd507f6a4d01576f..899203b54055f35dfe3bd9305f44df82cd0ae453 100644
--- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C
+++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C
@@ -50,6 +50,34 @@ namespace Foam
 
     const NamedEnum<directMappedPatchBase::sampleMode, 3>
         directMappedPatchBase::sampleModeNames_;
+
+
+    //- Private class for finding nearest
+    //  - point+local index
+    //  - sqr(distance)
+    //  - processor
+    typedef Tuple2<pointIndexHit, Tuple2<scalar, label> > nearInfo;
+
+    class nearestEqOp
+    {
+
+    public:
+
+        void operator()(nearInfo& x, const nearInfo& y) const
+        {
+            if (y.first().hit())
+            {
+                if (!x.first().hit())
+                {
+                    x = y;
+                }
+                else if (y.second().first() < x.second().first())
+                {
+                    x = y;
+                }
+            }
+        }
+    };
 }
 
 
@@ -70,7 +98,7 @@ void Foam::directMappedPatchBase::collectSamples
     labelListList globalFaces(Pstream::nProcs());
 
     globalFc[Pstream::myProcNo()] = patch_.faceCentres();
-    globalSamples[Pstream::myProcNo()] = globalFc[Pstream::myProcNo()]+offset_;
+    globalSamples[Pstream::myProcNo()] = globalFc[Pstream::myProcNo()]+offsets_;
     globalFaces[Pstream::myProcNo()] = identity(patch_.size());
 
     // Distribute to all processors
@@ -365,17 +393,17 @@ void Foam::directMappedPatchBase::calcMapping() const
 
     if
     (
-        offset_ == vector::zero
+        gAverage(mag(offsets_)) <= ROOTVSMALL
      && mode_ == NEARESTPATCHFACE
      && sampleRegion_ == patch_.boundaryMesh().mesh().name()
      && samplePatch_ == patch_.name()
     )
     {
         WarningIn("directMappedPatchBase::calcMapping() const")
-            << "Invalid offset " << offset_ << endl
+            << "Invalid offset " << offsets_ << endl
             << "Offset is the vector added to the patch face centres to"
             << " find the patch face supplying the data." << endl
-            << "Setting it to " << offset_
+            << "Setting it to " << offsets_
             << " on the same patch, on the same region"
             << " will find the faces themselves which does not make sense"
             << " for anything but testing." << endl
@@ -383,7 +411,7 @@ void Foam::directMappedPatchBase::calcMapping() const
             << "sampleRegion_:" << sampleRegion_ << endl
             << "mode_:" << sampleModeNames_[mode_] << endl
             << "samplePatch_:" << samplePatch_ << endl
-            << "offset_:" << offset_ << endl;
+            << "offsets_:" << offsets_ << endl;
     }
 
 
@@ -447,34 +475,38 @@ void Foam::directMappedPatchBase::calcMapping() const
     }
 
 
-    // Check that actual offset vector (sampleLocations - patchFc) is more or
-    // less constant.
-    if (Pstream::master())
-    {
-        const scalarField magOffset(mag(sampleLocations - patchFc));
-        const scalar avgOffset(average(magOffset));
-
-        forAll(magOffset, sampleI)
-        {
-            if (mag(magOffset[sampleI]-avgOffset) > max(SMALL, 0.001*avgOffset))
-            {
-                WarningIn("directMappedPatchBase::calcMapping() const")
-                    << "The actual cell/face centres picked up using offset "
-                    << offset_ << " are not" << endl
-                    << "    on a single plane."
-                    << " This might give numerical problems." << endl
-                    << "    At patchface " << patchFc[sampleI]
-                    << " the sampled cell/face " << sampleLocations[sampleI]
-                    << endl
-                    << "    is not on a plane " << avgOffset
-                    << " offset from the patch." << endl
-                    << "    You might want to shift your plane offset."
-                    << " Set the debug flag to get a dump of sampled cells."
-                    << endl;
-                break;
-            }
-        }
-    }
+    //// Check that actual offset vector (sampleLocations - patchFc) is more or
+    //// less constant.
+    //if (Pstream::master())
+    //{
+    //    const scalarField magOffset(mag(sampleLocations - patchFc));
+    //    const scalar avgOffset(average(magOffset));
+    //
+    //    forAll(magOffset, sampleI)
+    //    {
+    //        if
+    //        (
+    //            mag(magOffset[sampleI]-avgOffset)
+    //          > max(SMALL, 0.001*avgOffset)
+    //        )
+    //        {
+    //            WarningIn("directMappedPatchBase::calcMapping() const")
+    //                << "The actual cell/face centres picked up using offset "
+    //                << offsets_ << " are not" << endl
+    //                << "    on a single plane."
+    //                << " This might give numerical problems." << endl
+    //                << "    At patchface " << patchFc[sampleI]
+    //                << " the sampled cell/face " << sampleLocations[sampleI]
+    //                << endl
+    //                << "    is not on a plane " << avgOffset
+    //                << " offset from the patch." << endl
+    //                << "    You might want to shift your plane offset."
+    //                << " Set the debug flag to get a dump of sampled cells."
+    //                << endl;
+    //            break;
+    //        }
+    //    }
+    //}
 
 
     // Determine schedule.
@@ -556,14 +588,36 @@ void Foam::directMappedPatchBase::calcMapping() const
 
 Foam::directMappedPatchBase::directMappedPatchBase
 (
-     const polyPatch& pp
+    const polyPatch& pp
 )
 :
     patch_(pp),
     sampleRegion_(patch_.boundaryMesh().mesh().name()),
     mode_(NEARESTPATCHFACE),
     samplePatch_("none"),
+    uniformOffset_(true),
     offset_(vector::zero),
+    offsets_(pp.size(), offset_),
+    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
+    mapPtr_(NULL)
+{}
+
+
+Foam::directMappedPatchBase::directMappedPatchBase
+(
+    const polyPatch& pp,
+    const word& sampleRegion,
+    const sampleMode mode,
+    const word& samplePatch,
+    const vectorField& offsets
+)
+:
+    patch_(pp),
+    sampleRegion_(sampleRegion),
+    mode_(mode),
+    samplePatch_(samplePatch),
+    uniformOffset_(false),
+    offsets_(offsets),
     sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
     mapPtr_(NULL)
 {}
@@ -586,7 +640,19 @@ Foam::directMappedPatchBase::directMappedPatchBase
     ),
     mode_(sampleModeNames_.read(dict.lookup("sampleMode"))),
     samplePatch_(dict.lookup("samplePatch")),
-    offset_(dict.lookup("offset")),
+    uniformOffset_(dict.found("offset")),
+    offset_
+    (
+        uniformOffset_
+      ? point(dict.lookup("offset"))
+      : vector::zero
+    ),
+    offsets_
+    (
+        uniformOffset_
+      ? pointField(patch_.size(), offset_)
+      : dict.lookup("offsets")
+    ),
     sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
     mapPtr_(NULL)
 {}
@@ -602,7 +668,28 @@ Foam::directMappedPatchBase::directMappedPatchBase
     sampleRegion_(dmp.sampleRegion_),
     mode_(dmp.mode_),
     samplePatch_(dmp.samplePatch_),
+    uniformOffset_(dmp.uniformOffset_),
     offset_(dmp.offset_),
+    offsets_(dmp.offsets_),
+    sameRegion_(dmp.sameRegion_),
+    mapPtr_(NULL)
+{}
+
+
+Foam::directMappedPatchBase::directMappedPatchBase
+(
+    const polyPatch& pp,
+    const directMappedPatchBase& dmp,
+    const unallocLabelList& mapAddressing
+)
+:
+    patch_(pp),
+    sampleRegion_(dmp.sampleRegion_),
+    mode_(dmp.mode_),
+    samplePatch_(dmp.samplePatch_),
+    uniformOffset_(dmp.uniformOffset_),
+    offset_(dmp.offset_),
+    offsets_(dmp.offsets_, mapAddressing),
     sameRegion_(dmp.sameRegion_),
     mapPtr_(NULL)
 {}
@@ -660,7 +747,14 @@ void Foam::directMappedPatchBase::write(Ostream& os) const
         << token::END_STATEMENT << nl;
     os.writeKeyword("samplePatch") << samplePatch_
         << token::END_STATEMENT << nl;
-    os.writeKeyword("offset") << offset_ << token::END_STATEMENT << nl;
+    if (uniformOffset_)
+    {
+        os.writeKeyword("offset") << offset_ << token::END_STATEMENT << nl;
+    }
+    else
+    {
+        os.writeKeyword("offsets") << offsets_ << token::END_STATEMENT << nl;
+    }
 }
 
 
diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.H b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.H
index a3d936757eaa5fc190dfb7890ff2f0194406f9a8..3d6c0724290180d58ffca64d5fc8552ac2b6e727 100644
--- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.H
+++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.H
@@ -87,12 +87,18 @@ private:
         //- What to sample
         const sampleMode mode_;
 
-        //- Patch (only if NEARESTBOUNDARY)
+        //- Patch (only if NEARESTPATCHFACE)
         const word samplePatch_;
 
-        //- Offset vector
+        //- For backwards compatibility : reading/writing of uniform offset.
+        const bool uniformOffset_;
+
+        //- Offset vector (uniform)
         const vector offset_;
 
+        //- Offset vector
+        const vectorField offsets_;
+
         //- Same region
         const bool sameRegion_;
 
@@ -130,36 +136,6 @@ private:
         void calcMapping() const;
 
 
-    // Private class
-
-        //- Private class for finding nearest
-        //  - point+local index
-        //  - sqr(distance)
-        //  - processor
-        typedef Tuple2<pointIndexHit, Tuple2<scalar, label> > nearInfo;
-
-        class nearestEqOp
-        {
-
-        public:
-
-            void operator()(nearInfo& x, const nearInfo& y) const
-            {
-                if (y.first().hit())
-                {
-                    if (!x.first().hit())
-                    {
-                        x = y;
-                    }
-                    else if (y.second().first() < x.second().first())
-                    {
-                        x = y;
-                    }
-                }
-            }
-        };
-
-
 public:
 
     //- Runtime type information
@@ -168,15 +144,33 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from patch
         directMappedPatchBase(const polyPatch&);
 
+        //- Construct from components
+        directMappedPatchBase
+        (
+            const polyPatch& pp,
+            const word& sampleRegion,
+            const sampleMode sampleMode,
+            const word& samplePatch,
+            const vectorField& offset
+        );
+
         //- Construct from dictionary
         directMappedPatchBase(const polyPatch&, const dictionary&);
 
         //- Construct as copy, resetting patch
         directMappedPatchBase(const polyPatch&, const directMappedPatchBase&);
 
+        //- Construct as copy, resetting patch, map original data
+        directMappedPatchBase
+        (
+            const polyPatch&,
+            const directMappedPatchBase&,
+            const unallocLabelList& mapAddressing
+        );
+
 
     //- Destructor
     virtual ~directMappedPatchBase();
@@ -210,6 +204,12 @@ public:
             return offset_;
         }
 
+        //- Offset vector (from patch faces to destination mesh objects)
+        const vectorField& offsets() const
+        {
+            return offsets_;
+        }
+
         //- Return reference to the parallel distribution map
         const mapDistribute& map() const
         {
diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C
index 0058c84d0df42bd5448ff9e3ba15fd4f1674c95e..e16b8204f7c4ac72211ed3783173bb4c9312e385 100644
--- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C
+++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C
@@ -57,6 +57,31 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 {}
 
 
+Foam::directMappedPolyPatch::directMappedPolyPatch
+(
+    const word& name,
+    const label size,
+    const label start,
+    const label index,
+    const word& sampleRegion,
+    const directMappedPatchBase::sampleMode mode,
+    const word& samplePatch,
+    const vectorField& offset,
+    const polyBoundaryMesh& bm
+)
+:
+    polyPatch(name, size, start, index, bm),
+    directMappedPatchBase
+    (
+        static_cast<const polyPatch&>(*this),
+        sampleRegion,
+        mode,
+        samplePatch,
+        offset
+    )
+{}
+
+
 Foam::directMappedPolyPatch::directMappedPolyPatch
 (
     const word& name,
@@ -95,6 +120,20 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 {}
 
 
+Foam::directMappedPolyPatch::directMappedPolyPatch
+(
+    const directMappedPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    polyPatch(pp, bm, index, mapAddressing, newStart),
+    directMappedPatchBase(*this, pp, mapAddressing)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::directMappedPolyPatch::~directMappedPolyPatch()
diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H
index fdd4baec85c9547ded1b7911c765d696c0d10946..8cb5907c9514da934a6e05e6409bf31e9fece070 100644
--- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H
+++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H
@@ -101,6 +101,20 @@ public:
             const polyBoundaryMesh& bm
         );
 
+        //- Construct from components
+        directMappedPolyPatch
+        (
+            const word& name,
+            const label size,
+            const label start,
+            const label index,
+            const word& sampleRegion,
+            const directMappedPatchBase::sampleMode mode,
+            const word& samplePatch,
+            const vectorField& offset,
+            const polyBoundaryMesh& bm
+        );
+
         //- Construct from dictionary
         directMappedPolyPatch
         (
@@ -128,6 +142,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        directMappedPolyPatch
+        (
+            const directMappedPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -150,6 +174,29 @@ public:
             );
         }
 
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new directMappedPolyPatch
+                (
+                    *this,
+                    bm,
+                    index,
+                    mapAddressing,
+                    newStart
+                )
+            );
+        }
+
 
     //- Destructor
     virtual ~directMappedPolyPatch();
diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.C b/src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.C
index daaee3bcb4adf2b275627ede0a7489de6ca1115a..6d6038430d21ad7189b8b901e1dd4309ef37f735 100644
--- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.C
+++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.C
@@ -62,6 +62,31 @@ Foam::directMappedWallPolyPatch::directMappedWallPolyPatch
 {}
 
 
+Foam::directMappedWallPolyPatch::directMappedWallPolyPatch
+(
+    const word& name,
+    const label size,
+    const label start,
+    const label index,
+    const word& sampleRegion,
+    const directMappedPatchBase::sampleMode mode,
+    const word& samplePatch,
+    const vectorField& offset,
+    const polyBoundaryMesh& bm
+)
+:
+    wallPolyPatch(name, size, start, index, bm),
+    directMappedPatchBase
+    (
+        static_cast<const polyPatch&>(*this),
+        sampleRegion,
+        mode,
+        samplePatch,
+        offset
+    )
+{}
+
+
 Foam::directMappedWallPolyPatch::directMappedWallPolyPatch
 (
     const word& name,
@@ -100,6 +125,20 @@ Foam::directMappedWallPolyPatch::directMappedWallPolyPatch
 {}
 
 
+Foam::directMappedWallPolyPatch::directMappedWallPolyPatch
+(
+    const directMappedWallPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    wallPolyPatch(pp, bm, index, mapAddressing, newStart),
+    directMappedPatchBase(*this, pp, mapAddressing)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::directMappedWallPolyPatch::~directMappedWallPolyPatch()
diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.H b/src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.H
index 457552e42140f9a92a49b1ca3fd40ba99d1df1c9..0a2f5bdfd37f2536ed6feae5161df29d8d755614 100644
--- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.H
+++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.H
@@ -101,6 +101,20 @@ public:
             const polyBoundaryMesh& bm
         );
 
+        //- Construct from components
+        directMappedWallPolyPatch
+        (
+            const word& name,
+            const label size,
+            const label start,
+            const label index,
+            const word& sampleRegion,
+            const directMappedPatchBase::sampleMode mode,
+            const word& samplePatch,
+            const vectorField& offset,
+            const polyBoundaryMesh& bm
+        );
+
         //- Construct from dictionary
         directMappedWallPolyPatch
         (
@@ -128,6 +142,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        directMappedWallPolyPatch
+        (
+            const directMappedWallPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -157,6 +181,29 @@ public:
             );
         }
 
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new directMappedWallPolyPatch
+                (
+                    *this,
+                    bm,
+                    index,
+                    mapAddressing,
+                    newStart
+                )
+            );
+        }
+
 
     //- Destructor
     virtual ~directMappedWallPolyPatch();
diff --git a/src/meshTools/searchableSurface/searchableSurfaceCollection.H b/src/meshTools/searchableSurface/searchableSurfaceCollection.H
index dd4c5fdce215c5ccf3743a97541f9075114c668e..3f7ee84baa4012d2482007fa75434a8d1dbc4ca2 100644
--- a/src/meshTools/searchableSurface/searchableSurfaceCollection.H
+++ b/src/meshTools/searchableSurface/searchableSurfaceCollection.H
@@ -125,6 +125,31 @@ public:
 
     // Member Functions
 
+        //- scaling vector per subsurface
+        const vectorField& scale() const
+        {
+            return scale_;
+        }
+
+        //- scaling vector per subsurface
+        vectorField& scale()
+        {
+            return scale_;
+        }
+
+        //- coordinate system per subsurface
+        const PtrList<coordinateSystem>& transform() const
+        {
+            return transform_;
+        }
+
+        //- coordinate system per subsurface
+        PtrList<coordinateSystem>& transform()
+        {
+            return transform_;
+        }
+
+
         virtual const wordList& regions() const;
 
         //- Whether supports volume type below
diff --git a/src/meshTools/sets/cellSources/cellToCell/cellToCell.C b/src/meshTools/sets/cellSources/cellToCell/cellToCell.C
index 222d14965f636b8b345401c394139113ce2b354d..bbbcf7785ed85fa91b42ace67e743383db2c8607 100644
--- a/src/meshTools/sets/cellSources/cellToCell/cellToCell.C
+++ b/src/meshTools/sets/cellSources/cellToCell/cellToCell.C
@@ -106,7 +106,7 @@ void Foam::cellToCell::applyToSet
 {
     if ((action == topoSetSource::ADD) || (action == topoSetSource::NEW))
     {
-        Pout<< "    Adding all elements of cellSet " << setName_ << " ..."
+        Info<< "    Adding all elements of cellSet " << setName_ << " ..."
             << endl;
 
         // Load the set
@@ -116,7 +116,7 @@ void Foam::cellToCell::applyToSet
     }
     else if (action == topoSetSource::DELETE)
     {
-        Pout<< "    Removing all elements of cellSet " << setName_ << " ..."
+        Info<< "    Removing all elements of cellSet " << setName_ << " ..."
             << endl;
 
         // Load the set
diff --git a/src/meshTools/sets/cellZoneSources/setToCellZone/setToCellZone.C b/src/meshTools/sets/cellZoneSources/setToCellZone/setToCellZone.C
new file mode 100644
index 0000000000000000000000000000000000000000..5ff8597a863cf196a2746b1dd44173fdd4e0ae1e
--- /dev/null
+++ b/src/meshTools/sets/cellZoneSources/setToCellZone/setToCellZone.C
@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "setToCellZone.H"
+#include "polyMesh.H"
+#include "cellZoneSet.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(setToCellZone, 0);
+
+addToRunTimeSelectionTable(topoSetSource, setToCellZone, word);
+
+addToRunTimeSelectionTable(topoSetSource, setToCellZone, istream);
+
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::setToCellZone::usage_
+(
+    setToCellZone::typeName,
+    "\n    Usage: setToCellZone <cellSet>\n\n"
+    "    Select all cells in the cellSet.\n\n"
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::setToCellZone::setToCellZone
+(
+    const polyMesh& mesh,
+    const word& setName
+)
+:
+    topoSetSource(mesh),
+    setName_(setName)
+{}
+
+
+// Construct from dictionary
+Foam::setToCellZone::setToCellZone
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    topoSetSource(mesh),
+    setName_(dict.lookup("set"))
+{}
+
+
+// Construct from Istream
+Foam::setToCellZone::setToCellZone
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetSource(mesh),
+    setName_(checkIs(is))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::setToCellZone::~setToCellZone()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::setToCellZone::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if (!isA<cellZoneSet>(set))
+    {
+        WarningIn
+        (
+            "setToCellZone::applyToSet(const topoSetSource::setAction"
+            ", topoSet"
+        )   << "Operation only allowed on a cellZoneSet." << endl;
+    }
+    else
+    {
+        cellZoneSet& fzSet = refCast<cellZoneSet>(set);
+
+        if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+        {
+            Info<< "    Adding all cells from cellSet " << setName_
+                << " ..." << endl;
+
+            // Load the sets
+            cellSet fSet(mesh_, setName_);
+
+            // Start off from copy
+            DynamicList<label> newAddressing(fzSet.addressing());
+
+            forAllConstIter(cellSet, fSet, iter)
+            {
+                label cellI = iter.key();
+
+                if (!fzSet.found(cellI))
+                {
+                    newAddressing.append(cellI);
+                }
+            }
+
+            fzSet.addressing().transfer(newAddressing);
+            fzSet.updateSet();
+        }
+        else if (action == topoSetSource::DELETE)
+        {
+            Info<< "    Removing all cells from cellSet " << setName_
+                << " ..." << endl;
+
+            // Load the set
+            cellSet loadedSet(mesh_, setName_);
+
+            // Start off empty
+            DynamicList<label> newAddressing(fzSet.addressing().size());
+
+            forAll(fzSet.addressing(), i)
+            {
+                if (!loadedSet.found(fzSet.addressing()[i]))
+                {
+                    newAddressing.append(fzSet.addressing()[i]);
+                }
+            }
+            fzSet.addressing().transfer(newAddressing);
+            fzSet.updateSet();
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/cellZoneSources/setToCellZone/setToCellZone.H b/src/meshTools/sets/cellZoneSources/setToCellZone/setToCellZone.H
new file mode 100644
index 0000000000000000000000000000000000000000..cd60837d1ab4fda56009f0f16e4ecbc8e98eceb4
--- /dev/null
+++ b/src/meshTools/sets/cellZoneSources/setToCellZone/setToCellZone.H
@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::setToCellZone
+
+Description
+    A topoSetSource to select cells based on usage in a cellSet.
+
+SourceFiles
+    setToCellZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef setToCellZone_H
+#define setToCellZone_H
+
+#include "topoSetSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class setToCellZone Declaration
+\*---------------------------------------------------------------------------*/
+
+class setToCellZone
+:
+    public topoSetSource
+{
+    // Private data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- Name of set to use
+        word setName_;
+
+public:
+
+    //- Runtime type information
+    TypeName("setToCellZone");
+
+    // Constructors
+
+        //- Construct from components
+        setToCellZone
+        (
+            const polyMesh& mesh,
+            const word& setName
+        );
+
+        //- Construct from dictionary
+        setToCellZone
+        (
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct from Istream
+        setToCellZone
+        (
+            const polyMesh& mesh,
+            Istream&
+        );
+
+
+    // Destructor
+
+        virtual ~setToCellZone();
+
+
+    // Member Functions
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet&
+        ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/faceSources/cellToFace/cellToFace.C b/src/meshTools/sets/faceSources/cellToFace/cellToFace.C
index 59d3e2e9ebb95fdf5c7e8fc8c7c5d653fe1d8ad3..5b868c4e8c3b007990786af9668c1e5fac2825e3 100644
--- a/src/meshTools/sets/faceSources/cellToFace/cellToFace.C
+++ b/src/meshTools/sets/faceSources/cellToFace/cellToFace.C
@@ -217,14 +217,14 @@ void Foam::cellToFace::applyToSet
 {
     if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
     {
-        Pout<< "    Adding faces according to cellSet " << setName_
+        Info<< "    Adding faces according to cellSet " << setName_
             << " ..." << endl;
 
         combine(set, true);
     }
     else if (action == topoSetSource::DELETE)
     {
-        Pout<< "    Removing faces according to cellSet " << setName_
+        Info<< "    Removing faces according to cellSet " << setName_
             << " ..." << endl;
 
         combine(set, false);
diff --git a/src/meshTools/sets/faceZoneSources/faceZoneToFaceZone/faceZoneToFaceZone.C b/src/meshTools/sets/faceZoneSources/faceZoneToFaceZone/faceZoneToFaceZone.C
new file mode 100644
index 0000000000000000000000000000000000000000..49fa39eb09376aefd22d67b1a9f9880249216968
--- /dev/null
+++ b/src/meshTools/sets/faceZoneSources/faceZoneToFaceZone/faceZoneToFaceZone.C
@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "faceZoneToFaceZone.H"
+#include "polyMesh.H"
+#include "faceZoneSet.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(faceZoneToFaceZone, 0);
+
+addToRunTimeSelectionTable(topoSetSource, faceZoneToFaceZone, word);
+
+addToRunTimeSelectionTable(topoSetSource, faceZoneToFaceZone, istream);
+
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::faceZoneToFaceZone::usage_
+(
+    faceZoneToFaceZone::typeName,
+    "\n    Usage: faceZoneToFaceZone <faceZone>\n\n"
+    "    Select all faces in the faceZone\n\n"
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::faceZoneToFaceZone::faceZoneToFaceZone
+(
+    const polyMesh& mesh,
+    const word& setName
+)
+:
+    topoSetSource(mesh),
+    setName_(setName)
+{}
+
+
+// Construct from dictionary
+Foam::faceZoneToFaceZone::faceZoneToFaceZone
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    topoSetSource(mesh),
+    setName_(dict.lookup("zone"))
+{}
+
+
+// Construct from Istream
+Foam::faceZoneToFaceZone::faceZoneToFaceZone
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetSource(mesh),
+    setName_(checkIs(is))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::faceZoneToFaceZone::~faceZoneToFaceZone()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::faceZoneToFaceZone::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if (!isA<faceZoneSet>(set))
+    {
+        WarningIn
+        (
+            "faceZoneToFaceZone::applyToSet(const topoSetSource::setAction"
+            ", topoSet"
+        )   << "Operation only allowed on a faceZoneSet." << endl;
+    }
+    else
+    {
+        faceZoneSet& fSet = refCast<faceZoneSet>(set);
+
+        if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+        {
+            Info<< "    Adding all faces from faceZone " << setName_ << " ..."
+                << endl;
+
+            // Load the set
+            faceZoneSet loadedSet(mesh_, setName_);
+
+            DynamicList<label> newAddressing(fSet.addressing());
+            DynamicList<bool> newFlipMap(fSet.flipMap());
+
+            forAll(loadedSet.addressing(), i)
+            {
+                if (!fSet.found(loadedSet.addressing()[i]))
+                {
+                    newAddressing.append(loadedSet.addressing()[i]);
+                    newFlipMap.append(loadedSet.flipMap()[i]);
+                }
+            }
+            fSet.addressing().transfer(newAddressing);
+            fSet.flipMap().transfer(newFlipMap);
+            fSet.updateSet();
+        }
+        else if (action == topoSetSource::DELETE)
+        {
+            Info<< "    Removing all faces from faceZone " << setName_ << " ..."
+                << endl;
+
+            // Load the set
+            faceZoneSet loadedSet(mesh_, setName_);
+
+            DynamicList<label> newAddressing(fSet.addressing().size());
+            DynamicList<bool> newFlipMap(fSet.flipMap().size());
+
+            forAll(fSet.addressing(), i)
+            {
+                if (!loadedSet.found(fSet.addressing()[i]))
+                {
+                    newAddressing.append(fSet.addressing()[i]);
+                    newFlipMap.append(fSet.flipMap()[i]);
+                }
+            }
+            fSet.addressing().transfer(newAddressing);
+            fSet.flipMap().transfer(newFlipMap);
+            fSet.updateSet();
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/faceZoneSources/faceZoneToFaceZone/faceZoneToFaceZone.H b/src/meshTools/sets/faceZoneSources/faceZoneToFaceZone/faceZoneToFaceZone.H
new file mode 100644
index 0000000000000000000000000000000000000000..d8a0145d4d8b5a2da2d6bf50790b5ba8f5ca49fb
--- /dev/null
+++ b/src/meshTools/sets/faceZoneSources/faceZoneToFaceZone/faceZoneToFaceZone.H
@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::faceZoneToFaceZone
+
+Description
+    A topoSetSource to select faces based on usage in another faceSet.
+
+SourceFiles
+    faceZoneToFaceZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faceZoneToFaceZone_H
+#define faceZoneToFaceZone_H
+
+#include "topoSetSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class faceZoneToFaceZone Declaration
+\*---------------------------------------------------------------------------*/
+
+class faceZoneToFaceZone
+:
+    public topoSetSource
+{
+    // Private data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- Name of set to use
+        word setName_;
+
+public:
+
+    //- Runtime type information
+    TypeName("faceZoneToFaceZone");
+
+    // Constructors
+
+        //- Construct from components
+        faceZoneToFaceZone
+        (
+            const polyMesh& mesh,
+            const word& setName
+        );
+
+        //- Construct from dictionary
+        faceZoneToFaceZone
+        (
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct from Istream
+        faceZoneToFaceZone
+        (
+            const polyMesh& mesh,
+            Istream&
+        );
+
+
+    // Destructor
+
+        virtual ~faceZoneToFaceZone();
+
+
+    // Member Functions
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet&
+        ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/faceZoneSources/setToFaceZone/setToFaceZone.C b/src/meshTools/sets/faceZoneSources/setToFaceZone/setToFaceZone.C
new file mode 100644
index 0000000000000000000000000000000000000000..03382923165209827b23bd6e38fea13812f2becf
--- /dev/null
+++ b/src/meshTools/sets/faceZoneSources/setToFaceZone/setToFaceZone.C
@@ -0,0 +1,175 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "setToFaceZone.H"
+#include "polyMesh.H"
+#include "faceZoneSet.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(setToFaceZone, 0);
+
+addToRunTimeSelectionTable(topoSetSource, setToFaceZone, word);
+
+addToRunTimeSelectionTable(topoSetSource, setToFaceZone, istream);
+
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::setToFaceZone::usage_
+(
+    setToFaceZone::typeName,
+    "\n    Usage: setToFaceZone <faceSet>\n\n"
+    "    Select all faces in the faceSet."
+    " Sets flipMap.\n\n"
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::setToFaceZone::setToFaceZone
+(
+    const polyMesh& mesh,
+    const word& setName
+)
+:
+    topoSetSource(mesh),
+    setName_(setName)
+{}
+
+
+// Construct from dictionary
+Foam::setToFaceZone::setToFaceZone
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    topoSetSource(mesh),
+    setName_(dict.lookup("faceSet"))
+{}
+
+
+// Construct from Istream
+Foam::setToFaceZone::setToFaceZone
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetSource(mesh),
+    setName_(checkIs(is))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::setToFaceZone::~setToFaceZone()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::setToFaceZone::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if (!isA<faceZoneSet>(set))
+    {
+        WarningIn
+        (
+            "setToFaceZone::applyToSet(const topoSetSource::setAction"
+            ", topoSet"
+        )   << "Operation only allowed on a faceZoneSet." << endl;
+    }
+    else
+    {
+        faceZoneSet& fzSet = refCast<faceZoneSet>(set);
+
+        if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+        {
+            Info<< "    Adding all faces from faceSet " << setName_
+                << " ..." << endl;
+
+            // Load the sets
+            faceSet fSet(mesh_, setName_);
+
+            // Start off from copy
+            DynamicList<label> newAddressing(fzSet.addressing());
+            DynamicList<bool> newFlipMap(fzSet.flipMap());
+
+            forAllConstIter(faceSet, fSet, iter)
+            {
+                label faceI = iter.key();
+
+                if (!fzSet.found(faceI))
+                {
+                    newAddressing.append(faceI);
+                    newFlipMap.append(false);
+                }
+            }
+
+            fzSet.addressing().transfer(newAddressing);
+            fzSet.flipMap().transfer(newFlipMap);
+            fzSet.updateSet();
+        }
+        else if (action == topoSetSource::DELETE)
+        {
+            Info<< "    Removing all faces from faceSet " << setName_
+                << " ..." << endl;
+
+            // Load the set
+            faceSet loadedSet(mesh_, setName_);
+
+            // Start off empty
+            DynamicList<label> newAddressing(fzSet.addressing().size());
+            DynamicList<bool> newFlipMap(fzSet.flipMap().size());
+
+            forAll(fzSet.addressing(), i)
+            {
+                if (!loadedSet.found(fzSet.addressing()[i]))
+                {
+                    newAddressing.append(fzSet.addressing()[i]);
+                    newFlipMap.append(fzSet.flipMap()[i]);
+                }
+            }
+            fzSet.addressing().transfer(newAddressing);
+            fzSet.flipMap().transfer(newFlipMap);
+            fzSet.updateSet();
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/faceZoneSources/setToFaceZone/setToFaceZone.H b/src/meshTools/sets/faceZoneSources/setToFaceZone/setToFaceZone.H
new file mode 100644
index 0000000000000000000000000000000000000000..63c6e04efaf823dc6474052f79464d079d586b11
--- /dev/null
+++ b/src/meshTools/sets/faceZoneSources/setToFaceZone/setToFaceZone.H
@@ -0,0 +1,116 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::setToFaceZone
+
+Description
+    A topoSetSource to select faces based on usage in a faceSet. Sets flipMap
+    to true.
+
+SourceFiles
+    setToFaceZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef setToFaceZone_H
+#define setToFaceZone_H
+
+#include "topoSetSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class setToFaceZone Declaration
+\*---------------------------------------------------------------------------*/
+
+class setToFaceZone
+:
+    public topoSetSource
+{
+    // Private data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- Name of set to use
+        word setName_;
+
+public:
+
+    //- Runtime type information
+    TypeName("setToFaceZone");
+
+    // Constructors
+
+        //- Construct from components
+        setToFaceZone
+        (
+            const polyMesh& mesh,
+            const word& setName
+        );
+
+        //- Construct from dictionary
+        setToFaceZone
+        (
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct from Istream
+        setToFaceZone
+        (
+            const polyMesh& mesh,
+            Istream&
+        );
+
+
+    // Destructor
+
+        virtual ~setToFaceZone();
+
+
+    // Member Functions
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet&
+        ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.C b/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.C
new file mode 100644
index 0000000000000000000000000000000000000000..11e2eadd483d9db9bc7c9a6bcc38fd7c32e5c2f9
--- /dev/null
+++ b/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.C
@@ -0,0 +1,222 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "setsToFaceZone.H"
+#include "polyMesh.H"
+#include "faceZoneSet.H"
+#include "cellSet.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(setsToFaceZone, 0);
+
+addToRunTimeSelectionTable(topoSetSource, setsToFaceZone, word);
+
+addToRunTimeSelectionTable(topoSetSource, setsToFaceZone, istream);
+
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::setsToFaceZone::usage_
+(
+    setsToFaceZone::typeName,
+    "\n    Usage: setsToFaceZone <faceSet> <slaveCellSet>\n\n"
+    "    Select all faces in the faceSet."
+    " Orientated so slave side is in cellSet.\n\n"
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::setsToFaceZone::setsToFaceZone
+(
+    const polyMesh& mesh,
+    const word& faceSetName,
+    const word& cellSetName
+)
+:
+    topoSetSource(mesh),
+    faceSetName_(faceSetName),
+    cellSetName_(cellSetName)
+{}
+
+
+// Construct from dictionary
+Foam::setsToFaceZone::setsToFaceZone
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    topoSetSource(mesh),
+    faceSetName_(dict.lookup("faceSet")),
+    cellSetName_(dict.lookup("cellSet"))
+{}
+
+
+// Construct from Istream
+Foam::setsToFaceZone::setsToFaceZone
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetSource(mesh),
+    faceSetName_(checkIs(is)),
+    cellSetName_(checkIs(is))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::setsToFaceZone::~setsToFaceZone()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::setsToFaceZone::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if (!isA<faceZoneSet>(set))
+    {
+        WarningIn
+        (
+            "setsToFaceZone::applyToSet(const topoSetSource::setAction"
+            ", topoSet"
+        )   << "Operation only allowed on a faceZoneSet." << endl;
+    }
+    else
+    {
+        faceZoneSet& fzSet = refCast<faceZoneSet>(set);
+
+        if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+        {
+            Info<< "    Adding all faces from faceSet " << faceSetName_
+                << " ..." << endl;
+
+            // Load the sets
+            faceSet fSet(mesh_, faceSetName_);
+            cellSet cSet(mesh_, cellSetName_);
+
+            // Start off from copy
+            DynamicList<label> newAddressing(fzSet.addressing());
+            DynamicList<bool> newFlipMap(fzSet.flipMap());
+
+            forAllConstIter(faceSet, fSet, iter)
+            {
+                label faceI = iter.key();
+
+                if (!fzSet.found(faceI))
+                {
+                    bool flip = false;
+
+                    label own = mesh_.faceOwner()[faceI];
+                    bool ownFound = cSet.found(own);
+
+                    if (mesh_.isInternalFace(faceI))
+                    {
+                        label nei = mesh_.faceNeighbour()[faceI];
+                        bool neiFound = cSet.found(nei);
+
+                        if (ownFound && !neiFound)
+                        {
+                            flip = false;
+                        }
+                        else if(!ownFound && neiFound)
+                        {
+                            flip = true;
+                        }
+                        else
+                        {
+                            WarningIn
+                            (
+                                "setsToFaceZone::applyToSet"
+                                "(const topoSetSource::setAction, topoSet)"
+                            )   << "One of owner or neighbour of internal face "
+                                << faceI << " should be in cellSet "
+                                << cSet.name()
+                                << " to be able to determine orientation."
+                                << endl
+                                << "Face:" << faceI << " own:" << own
+                                << " OwnInCellSet:" << ownFound
+                                << " nei:" << nei
+                                << " NeiInCellSet:" << neiFound
+                                << endl;
+                        }
+                    }
+                    else
+                    {
+                        flip = !ownFound;
+                    }
+
+                    newAddressing.append(faceI);
+                    newFlipMap.append(flip);
+                }
+            }
+
+            fzSet.addressing().transfer(newAddressing);
+            fzSet.flipMap().transfer(newFlipMap);
+            fzSet.updateSet();
+        }
+        else if (action == topoSetSource::DELETE)
+        {
+            Info<< "    Removing all faces from faceSet " << faceSetName_
+                << " ..." << endl;
+
+            // Load the set
+            faceZoneSet loadedSet(mesh_, faceSetName_);
+
+            // Start off empty
+            DynamicList<label> newAddressing(fzSet.addressing().size());
+            DynamicList<bool> newFlipMap(fzSet.flipMap().size());
+
+            forAll(fzSet.addressing(), i)
+            {
+                if (!loadedSet.found(fzSet.addressing()[i]))
+                {
+                    newAddressing.append(fzSet.addressing()[i]);
+                    newFlipMap.append(fzSet.flipMap()[i]);
+                }
+            }
+            fzSet.addressing().transfer(newAddressing);
+            fzSet.flipMap().transfer(newFlipMap);
+            fzSet.updateSet();
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.H b/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.H
new file mode 100644
index 0000000000000000000000000000000000000000..512c8243d8c210bbcc83cc59a66d30f8dfacf746
--- /dev/null
+++ b/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.H
@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::setsToFaceZone
+
+Description
+    A topoSetSource to select faces based on usage in a faceSet and cellSet
+
+SourceFiles
+    setsToFaceZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef setsToFaceZone_H
+#define setsToFaceZone_H
+
+#include "topoSetSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class setsToFaceZone Declaration
+\*---------------------------------------------------------------------------*/
+
+class setsToFaceZone
+:
+    public topoSetSource
+{
+    // Private data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- Name of set to use
+        word faceSetName_;
+
+        //- Name of set to use
+        word cellSetName_;
+
+public:
+
+    //- Runtime type information
+    TypeName("setsToFaceZone");
+
+    // Constructors
+
+        //- Construct from components
+        setsToFaceZone
+        (
+            const polyMesh& mesh,
+            const word& faceSetName,
+            const word& cellSetName
+        );
+
+        //- Construct from dictionary
+        setsToFaceZone
+        (
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct from Istream
+        setsToFaceZone
+        (
+            const polyMesh& mesh,
+            Istream&
+        );
+
+
+    // Destructor
+
+        virtual ~setsToFaceZone();
+
+
+    // Member Functions
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet&
+        ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/pointSources/cellToPoint/cellToPoint.C b/src/meshTools/sets/pointSources/cellToPoint/cellToPoint.C
index 1753f84848f8e90f007a9e2a462fb2837eaa4be1..e96f10945642a1af0f2418bda3ae55350099ec08 100644
--- a/src/meshTools/sets/pointSources/cellToPoint/cellToPoint.C
+++ b/src/meshTools/sets/pointSources/cellToPoint/cellToPoint.C
@@ -151,13 +151,13 @@ void Foam::cellToPoint::applyToSet
 {
     if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
     {
-        Pout<< "    Adding from " << setName_ << " ..." << endl;
+        Info<< "    Adding from " << setName_ << " ..." << endl;
 
         combine(set, true);
     }
     else if (action == topoSetSource::DELETE)
     {
-        Pout<< "    Removing from " << setName_ << " ..." << endl;
+        Info<< "    Removing from " << setName_ << " ..." << endl;
 
         combine(set, false);
     }
diff --git a/src/meshTools/sets/pointZoneSources/setToPointZone/setToPointZone.C b/src/meshTools/sets/pointZoneSources/setToPointZone/setToPointZone.C
new file mode 100644
index 0000000000000000000000000000000000000000..80b76f5eed4fd55ca9cc114f17879c5980f89d58
--- /dev/null
+++ b/src/meshTools/sets/pointZoneSources/setToPointZone/setToPointZone.C
@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "setToPointZone.H"
+#include "polyMesh.H"
+#include "pointZoneSet.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(setToPointZone, 0);
+
+addToRunTimeSelectionTable(topoSetSource, setToPointZone, word);
+
+addToRunTimeSelectionTable(topoSetSource, setToPointZone, istream);
+
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::setToPointZone::usage_
+(
+    setToPointZone::typeName,
+    "\n    Usage: setToPointZone <pointSet>\n\n"
+    "    Select all points in the pointSet.\n\n"
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::setToPointZone::setToPointZone
+(
+    const polyMesh& mesh,
+    const word& setName
+)
+:
+    topoSetSource(mesh),
+    setName_(setName)
+{}
+
+
+// Construct from dictionary
+Foam::setToPointZone::setToPointZone
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    topoSetSource(mesh),
+    setName_(dict.lookup("set"))
+{}
+
+
+// Construct from Istream
+Foam::setToPointZone::setToPointZone
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetSource(mesh),
+    setName_(checkIs(is))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::setToPointZone::~setToPointZone()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::setToPointZone::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if (!isA<pointZoneSet>(set))
+    {
+        WarningIn
+        (
+            "setToPointZone::applyToSet(const topoSetSource::setAction"
+            ", topoSet"
+        )   << "Operation only allowed on a pointZoneSet." << endl;
+    }
+    else
+    {
+        pointZoneSet& fzSet = refCast<pointZoneSet>(set);
+
+        if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+        {
+            Info<< "    Adding all points from pointSet " << setName_
+                << " ..." << endl;
+
+            // Load the sets
+            pointSet fSet(mesh_, setName_);
+
+            // Start off from copy
+            DynamicList<label> newAddressing(fzSet.addressing());
+
+            forAllConstIter(pointSet, fSet, iter)
+            {
+                label pointI = iter.key();
+
+                if (!fzSet.found(pointI))
+                {
+                    newAddressing.append(pointI);
+                }
+            }
+
+            fzSet.addressing().transfer(newAddressing);
+            fzSet.updateSet();
+        }
+        else if (action == topoSetSource::DELETE)
+        {
+            Info<< "    Removing all points from pointSet " << setName_
+                << " ..." << endl;
+
+            // Load the set
+            pointSet loadedSet(mesh_, setName_);
+
+            // Start off empty
+            DynamicList<label> newAddressing(fzSet.addressing().size());
+
+            forAll(fzSet.addressing(), i)
+            {
+                if (!loadedSet.found(fzSet.addressing()[i]))
+                {
+                    newAddressing.append(fzSet.addressing()[i]);
+                }
+            }
+            fzSet.addressing().transfer(newAddressing);
+            fzSet.updateSet();
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/pointZoneSources/setToPointZone/setToPointZone.H b/src/meshTools/sets/pointZoneSources/setToPointZone/setToPointZone.H
new file mode 100644
index 0000000000000000000000000000000000000000..d13a3fe68f58d1f68e645ef0b890c63d0af26253
--- /dev/null
+++ b/src/meshTools/sets/pointZoneSources/setToPointZone/setToPointZone.H
@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::setToPointZone
+
+Description
+    A topoSetSource to select points based on usage in a pointSet.
+
+SourceFiles
+    setToPointZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef setToPointZone_H
+#define setToPointZone_H
+
+#include "topoSetSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class setToPointZone Declaration
+\*---------------------------------------------------------------------------*/
+
+class setToPointZone
+:
+    public topoSetSource
+{
+    // Private data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- Name of set to use
+        word setName_;
+
+public:
+
+    //- Runtime type information
+    TypeName("setToPointZone");
+
+    // Constructors
+
+        //- Construct from components
+        setToPointZone
+        (
+            const polyMesh& mesh,
+            const word& setName
+        );
+
+        //- Construct from dictionary
+        setToPointZone
+        (
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct from Istream
+        setToPointZone
+        (
+            const polyMesh& mesh,
+            Istream&
+        );
+
+
+    // Destructor
+
+        virtual ~setToPointZone();
+
+
+    // Member Functions
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet&
+        ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/topoSetSource/topoSetSource.C b/src/meshTools/sets/topoSetSource/topoSetSource.C
index 7579ffffad5c33ba25c2efd399cb528606e96b8c..8135859129cd9ace2980ee02d1e6e65e11c96ef6 100644
--- a/src/meshTools/sets/topoSetSource/topoSetSource.C
+++ b/src/meshTools/sets/topoSetSource/topoSetSource.C
@@ -101,7 +101,7 @@ autoPtr<topoSetSource> topoSetSource::New
 Foam::HashTable<Foam::string>* Foam::topoSetSource::usageTablePtr_ = NULL;
 
 template<>
-const char* Foam::NamedEnum<Foam::topoSetSource::setAction, 7>::names[] =
+const char* Foam::NamedEnum<Foam::topoSetSource::setAction, 8>::names[] =
 {
     "clear",
     "new",
@@ -109,11 +109,12 @@ const char* Foam::NamedEnum<Foam::topoSetSource::setAction, 7>::names[] =
     "add",
     "delete",
     "subset",
-    "list"
+    "list",
+    "remove"
 };
 
 
-const Foam::NamedEnum<Foam::topoSetSource::setAction, 7>
+const Foam::NamedEnum<Foam::topoSetSource::setAction, 8>
     Foam::topoSetSource::actionNames_;
 
 
diff --git a/src/meshTools/sets/topoSetSource/topoSetSource.H b/src/meshTools/sets/topoSetSource/topoSetSource.H
index 26fa39af9eb0041ab9c20e88eeed55a3704f8a5c..b1c68a14c14e058a6ceb5b8354ef29f2ba9cc6f4 100644
--- a/src/meshTools/sets/topoSetSource/topoSetSource.H
+++ b/src/meshTools/sets/topoSetSource/topoSetSource.H
@@ -77,7 +77,8 @@ public:
             ADD,
             DELETE,
             SUBSET,
-            LIST
+            LIST,
+            REMOVE
         };
 
 protected:
@@ -120,7 +121,7 @@ protected:
 
 private:
 
-        static const NamedEnum<setAction, 7> actionNames_;
+        static const NamedEnum<setAction, 8> actionNames_;
 
         static const string illegalSource_;
 
diff --git a/src/meshTools/sets/topoSets/cellSet.C b/src/meshTools/sets/topoSets/cellSet.C
index 4cff9f2908164ccbf8da66ccb31c4013ab22fc66..dabf987a12057f1e5f6cbe11bb62c66d55b97995 100644
--- a/src/meshTools/sets/topoSets/cellSet.C
+++ b/src/meshTools/sets/topoSets/cellSet.C
@@ -41,6 +41,7 @@ defineTypeNameAndDebug(cellSet, 0);
 
 addToRunTimeSelectionTable(topoSet, cellSet, word);
 addToRunTimeSelectionTable(topoSet, cellSet, size);
+addToRunTimeSelectionTable(topoSet, cellSet, set);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -78,6 +79,18 @@ cellSet::cellSet
 {}
 
 
+cellSet::cellSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const topoSet& set,
+    writeOption w
+)
+:
+    topoSet(mesh, name, set, w)
+{}
+
+
 cellSet::cellSet
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/topoSets/cellSet.H b/src/meshTools/sets/topoSets/cellSet.H
index 627e80e1420bbdcb7d262be4b69a2472a69b8b60..b931ac22f142fa030a4bab64d9b2cc808733d581 100644
--- a/src/meshTools/sets/topoSets/cellSet.H
+++ b/src/meshTools/sets/topoSets/cellSet.H
@@ -88,6 +88,15 @@ public:
             writeOption w=NO_WRITE
         );
 
+        //- Construct from existing set
+        cellSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const topoSet&,
+            writeOption w=NO_WRITE
+        );
+
         //- Construct from labelHashSet
         cellSet
         (
diff --git a/src/meshTools/sets/topoSets/cellZoneSet.C b/src/meshTools/sets/topoSets/cellZoneSet.C
new file mode 100644
index 0000000000000000000000000000000000000000..613a23bc8c7d59eed456b26e6db1544181c54261
--- /dev/null
+++ b/src/meshTools/sets/topoSets/cellZoneSet.C
@@ -0,0 +1,314 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "cellZoneSet.H"
+#include "mapPolyMesh.H"
+#include "polyMesh.H"
+#include "processorPolyPatch.H"
+#include "cyclicPolyPatch.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(cellZoneSet, 0);
+
+addToRunTimeSelectionTable(topoSet, cellZoneSet, word);
+addToRunTimeSelectionTable(topoSet, cellZoneSet, size);
+addToRunTimeSelectionTable(topoSet, cellZoneSet, set);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void cellZoneSet::updateSet()
+{
+    labelList order;
+    sortedOrder(addressing_, order);
+    inplaceReorder(order, addressing_);
+
+    cellSet::clearStorage();
+    cellSet::resize(2*addressing_.size());
+    forAll(addressing_, i)
+    {
+        cellSet::insert(addressing_[i]);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+cellZoneSet::cellZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    readOption r,
+    writeOption w
+)
+:
+    cellSet(mesh, name, 1000),  // do not read cellSet
+    mesh_(mesh),
+    addressing_(0)
+{
+    const cellZoneMesh& cellZones = mesh.cellZones();
+    label zoneID = cellZones.findZoneID(name);
+
+    if
+    (
+        (r == IOobject::MUST_READ)
+     || (r == IOobject::READ_IF_PRESENT && zoneID != -1)
+    )
+    {
+        const cellZone& fz = cellZones[zoneID];
+        addressing_ = fz;
+    }
+
+    updateSet();
+
+    check(mesh.nCells());
+}
+
+
+cellZoneSet::cellZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const label size,
+    writeOption w
+)
+:
+    cellSet(mesh, name, size, w),
+    mesh_(mesh),
+    addressing_(0)
+{
+    updateSet();
+}
+
+
+cellZoneSet::cellZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const topoSet& set,
+    writeOption w
+)
+:
+    cellSet(mesh, name, set.size(), w),
+    mesh_(mesh),
+    addressing_(refCast<const cellZoneSet>(set).addressing())
+{
+    updateSet();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+cellZoneSet::~cellZoneSet()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void cellZoneSet::invert(const label maxLen)
+{
+    label n = 0;
+
+    for (label cellI = 0; cellI < maxLen; cellI++)
+    {
+        if (!found(cellI))
+        {
+            addressing_[n] = cellI;
+            n++;
+        }
+    }
+    addressing_.setSize(n);
+    updateSet();
+}
+
+
+void cellZoneSet::subset(const topoSet& set)
+{
+    DynamicList<label> newAddressing(addressing_.size());
+
+    const cellZoneSet& fSet = refCast<const cellZoneSet>(set);
+
+    forAll(fSet.addressing(), i)
+    {
+        label cellI = fSet.addressing()[i];
+
+        if (found(cellI))
+        {
+            newAddressing.append(cellI);
+        }
+    }
+
+    addressing_.transfer(newAddressing);
+    updateSet();
+}
+
+
+void cellZoneSet::addSet(const topoSet& set)
+{
+    DynamicList<label> newAddressing(addressing_);
+
+    const cellZoneSet& fSet = refCast<const cellZoneSet>(set);
+
+    forAll(fSet.addressing(), i)
+    {
+        label cellI = fSet.addressing()[i];
+
+        if (!found(cellI))
+        {
+            newAddressing.append(cellI);
+        }
+    }
+
+    addressing_.transfer(newAddressing);
+    updateSet();
+}
+
+
+void cellZoneSet::deleteSet(const topoSet& set)
+{
+    DynamicList<label> newAddressing(addressing_.size());
+
+    const cellZoneSet& fSet = refCast<const cellZoneSet>(set);
+
+    forAll(addressing_, i)
+    {
+        label cellI = addressing_[i];
+
+        if (!fSet.found(cellI))
+        {
+            // Not found in fSet so add
+            newAddressing.append(cellI);
+        }
+    }
+
+    addressing_.transfer(newAddressing);
+    updateSet();
+}
+
+
+void cellZoneSet::sync(const polyMesh& mesh)
+{}
+
+
+label cellZoneSet::maxSize(const polyMesh& mesh) const
+{
+    return mesh.nCells();
+}
+
+
+//- Write using given format, version and compression
+bool cellZoneSet::writeObject
+(
+    IOstream::streamFormat s,
+    IOstream::versionNumber v,
+    IOstream::compressionType c
+) const
+{
+    // Write shadow cellSet
+    word oldTypeName = typeName;
+    const_cast<word&>(type()) = cellSet::typeName;
+    bool ok = cellSet::writeObject(s, v, c);
+    const_cast<word&>(type()) = oldTypeName;
+
+    // Modify cellZone
+    cellZoneMesh& cellZones = const_cast<polyMesh&>(mesh_).cellZones();
+    label zoneID = cellZones.findZoneID(name());
+
+    if (zoneID == -1)
+    {
+        zoneID = cellZones.size();
+
+        cellZones.setSize(zoneID+1);
+        cellZones.set
+        (
+            zoneID,
+            new cellZone
+            (
+                name(),
+                addressing_,
+                zoneID,
+                cellZones
+            )
+        );
+    }
+    else
+    {
+        cellZones[zoneID] = addressing_;
+    }
+    cellZones.clearAddressing();
+
+    return ok && cellZones.write();
+}
+
+
+void cellZoneSet::updateMesh(const mapPolyMesh& morphMap)
+{
+    // cellZone
+    labelList newAddressing(addressing_.size());
+
+    label n = 0;
+    forAll(addressing_, i)
+    {
+        label cellI = addressing_[i];
+        label newCellI = morphMap.reverseCellMap()[cellI];
+        if (newCellI >= 0)
+        {
+            newAddressing[n] = newCellI;
+            n++;
+        }
+    }
+    newAddressing.setSize(n);
+
+    addressing_.transfer(newAddressing);
+
+    updateSet();
+}
+
+
+void cellZoneSet::writeDebug
+(
+    Ostream& os,
+    const primitiveMesh& mesh,
+    const label maxLen
+) const
+{
+    cellSet::writeDebug(os, mesh, maxLen);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/topoSets/cellZoneSet.H b/src/meshTools/sets/topoSets/cellZoneSet.H
new file mode 100644
index 0000000000000000000000000000000000000000..7ddb48d1146ab8d5e78003b49722ce6f18de0f76
--- /dev/null
+++ b/src/meshTools/sets/topoSets/cellZoneSet.H
@@ -0,0 +1,170 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::cellZoneSet
+
+Description
+    Like cellSet but updates cellZone when writing.
+
+SourceFiles
+    cellZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cellZoneSet_H
+#define cellZoneSet_H
+
+#include "cellSet.H"
+#include "boolList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class cellZoneSet Declaration
+\*---------------------------------------------------------------------------*/
+
+class cellZoneSet
+:
+    public cellSet
+{
+    // Private data
+
+        const polyMesh& mesh_;
+
+        labelList addressing_;
+
+public:
+
+    //- Runtime type information
+    TypeName("cellZoneSet");
+
+
+    // Constructors
+
+        //- Construct from objectRegistry and name
+        cellZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            readOption r=MUST_READ,
+            writeOption w=NO_WRITE
+        );
+
+        //- Construct from additional size of labelHashSet
+        cellZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const label,
+            writeOption w=NO_WRITE
+        );
+
+        //- Construct from existing set
+        cellZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const topoSet&,
+            writeOption w=NO_WRITE
+        );
+
+
+
+    // Destructor
+
+        virtual ~cellZoneSet();
+
+
+    // Member functions
+
+        const labelList& addressing() const
+        {
+            return addressing_;
+        }
+
+        labelList& addressing()
+        {
+            return addressing_;
+        }
+
+        //- Sort addressing and make cellSet part consistent with addressing
+        void updateSet();
+
+        //- Invert contents. (insert all members 0..maxLen-1 which were not in
+        //  set)
+        virtual void invert(const label maxLen);
+
+        //- Subset contents. Only elements present in both sets remain.
+        virtual void subset(const topoSet& set);
+
+        //- Add elements present in set.
+        virtual void addSet(const topoSet& set);
+
+        //- Delete elements present in set.
+        virtual void deleteSet(const topoSet& set);
+
+        //- Sync cellZoneSet across coupled patches.
+        virtual void sync(const polyMesh& mesh);
+
+        //- Write maxLen items with label and coordinates. 
+        virtual void writeDebug
+        (
+            Ostream& os,
+            const primitiveMesh&,
+            const label maxLen
+        ) const;
+
+        //- Write cellZone
+        virtual bool writeObject
+        (
+            IOstream::streamFormat,
+            IOstream::versionNumber,
+            IOstream::compressionType
+        ) const;
+
+        //- Update any stored data for new labels
+        virtual void updateMesh(const mapPolyMesh& morphMap);
+
+        //- Return max index+1.
+        virtual label maxSize(const polyMesh& mesh) const;
+
+
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/topoSets/faceSet.C b/src/meshTools/sets/topoSets/faceSet.C
index e4a84a14e9b8b57900ee5d3e4d700072a604684e..c5bd79297319bfcd907df167b77d526dd386a05e 100644
--- a/src/meshTools/sets/topoSets/faceSet.C
+++ b/src/meshTools/sets/topoSets/faceSet.C
@@ -43,6 +43,7 @@ defineTypeNameAndDebug(faceSet, 0);
 
 addToRunTimeSelectionTable(topoSet, faceSet, word);
 addToRunTimeSelectionTable(topoSet, faceSet, size);
+addToRunTimeSelectionTable(topoSet, faceSet, set);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -79,6 +80,18 @@ faceSet::faceSet
 {}
 
 
+faceSet::faceSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const topoSet& set,
+    writeOption w
+)
+:
+    topoSet(mesh, name, set, w)
+{}
+
+
 faceSet::faceSet
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/topoSets/faceSet.H b/src/meshTools/sets/topoSets/faceSet.H
index aefb3380f867cb59876d60bb35e2e26d44f1e731..11e7387b9723b1371f3600c265e601cd15c3727e 100644
--- a/src/meshTools/sets/topoSets/faceSet.H
+++ b/src/meshTools/sets/topoSets/faceSet.H
@@ -83,6 +83,15 @@ public:
             writeOption w=NO_WRITE
         );
 
+        //- Construct from existing set
+        faceSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const topoSet&,
+            writeOption w=NO_WRITE
+        );
+
         //- Construct from additional labelHashSet
         faceSet
         (
diff --git a/src/meshTools/sets/topoSets/faceZoneSet.C b/src/meshTools/sets/topoSets/faceZoneSet.C
new file mode 100644
index 0000000000000000000000000000000000000000..85069a83ae2d6de8e9282ff2ae0fab82fd537b27
--- /dev/null
+++ b/src/meshTools/sets/topoSets/faceZoneSet.C
@@ -0,0 +1,413 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "faceZoneSet.H"
+#include "mapPolyMesh.H"
+#include "polyMesh.H"
+#include "processorPolyPatch.H"
+#include "cyclicPolyPatch.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(faceZoneSet, 0);
+
+addToRunTimeSelectionTable(topoSet, faceZoneSet, word);
+addToRunTimeSelectionTable(topoSet, faceZoneSet, size);
+addToRunTimeSelectionTable(topoSet, faceZoneSet, set);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void faceZoneSet::updateSet()
+{
+    labelList order;
+    sortedOrder(addressing_, order);
+    inplaceReorder(order, addressing_);
+    inplaceReorder(order, flipMap_);
+
+    faceSet::clearStorage();
+    faceSet::resize(2*addressing_.size());
+    forAll(addressing_, i)
+    {
+        faceSet::insert(addressing_[i]);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+faceZoneSet::faceZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    readOption r,
+    writeOption w
+)
+:
+    faceSet(mesh, name, 1000),  // do not read faceSet
+    mesh_(mesh),
+    addressing_(0),
+    flipMap_(0)
+{
+    const faceZoneMesh& faceZones = mesh.faceZones();
+    label zoneID = faceZones.findZoneID(name);
+
+    if
+    (
+        (r == IOobject::MUST_READ)
+     || (r == IOobject::READ_IF_PRESENT && zoneID != -1)
+    )
+    {
+        const faceZone& fz = faceZones[zoneID];
+        addressing_ = fz;
+        flipMap_ = fz.flipMap();
+    }
+
+    updateSet();
+
+    check(mesh.nFaces());
+}
+
+
+faceZoneSet::faceZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const label size,
+    writeOption w
+)
+:
+    faceSet(mesh, name, size, w),
+    mesh_(mesh),
+    addressing_(0),
+    flipMap_(0)
+{
+    updateSet();
+}
+
+
+faceZoneSet::faceZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const topoSet& set,
+    writeOption w
+)
+:
+    faceSet(mesh, name, set.size(), w),
+    mesh_(mesh),
+    addressing_(refCast<const faceZoneSet>(set).addressing()),
+    flipMap_(refCast<const faceZoneSet>(set).flipMap())
+{
+    updateSet();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+faceZoneSet::~faceZoneSet()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void faceZoneSet::invert(const label maxLen)
+{
+    label n = 0;
+
+    for (label faceI = 0; faceI < maxLen; faceI++)
+    {
+        if (!found(faceI))
+        {
+            addressing_[n] = faceI;
+            flipMap_[n] = false;         //? or true?
+            n++;
+        }
+    }
+    addressing_.setSize(n);
+    flipMap_.setSize(n);
+    updateSet();
+}
+
+
+void faceZoneSet::subset(const topoSet& set)
+{
+    label nConflict = 0;
+
+    DynamicList<label> newAddressing(addressing_.size());
+    DynamicList<bool> newFlipMap(flipMap_.size());
+
+    Map<label> faceToIndex(addressing_.size());
+    forAll(addressing_, i)
+    {
+        faceToIndex.insert(addressing_[i], i);
+    }
+
+    const faceZoneSet& fSet = refCast<const faceZoneSet>(set);
+
+    forAll(fSet.addressing(), i)
+    {
+        label faceI = fSet.addressing()[i];
+
+        Map<label>::const_iterator iter = faceToIndex.find(faceI);
+
+        if (iter != faceToIndex.end())
+        {
+            label index = iter();
+
+            if (fSet.flipMap()[i] != flipMap_[index])
+            {
+                nConflict++;
+            }
+            newAddressing.append(faceI);
+            newFlipMap.append(flipMap_[index]);
+        }
+    }
+
+    if (nConflict > 0)
+    {
+        WarningIn(" faceZoneSet::subset(const topoSet&)")
+            << "subset : there are " << nConflict
+            << " faces with different orientation in faceZonesSets "
+            << name() << " and " << set.name() << endl;
+    }
+
+    addressing_.transfer(newAddressing);
+    flipMap_.transfer(newFlipMap);
+    updateSet();
+}
+
+
+void faceZoneSet::addSet(const topoSet& set)
+{
+    label nConflict = 0;
+
+    DynamicList<label> newAddressing(addressing_);
+    DynamicList<bool> newFlipMap(flipMap_);
+
+    Map<label> faceToIndex(addressing_.size());
+    forAll(addressing_, i)
+    {
+        faceToIndex.insert(addressing_[i], i);
+    }
+
+    const faceZoneSet& fSet = refCast<const faceZoneSet>(set);
+
+    forAll(fSet.addressing(), i)
+    {
+        label faceI = fSet.addressing()[i];
+
+        Map<label>::const_iterator iter = faceToIndex.find(faceI);
+
+        if (iter != faceToIndex.end())
+        {
+            label index = iter();
+
+            if (fSet.flipMap()[i] != flipMap_[index])
+            {
+                nConflict++;
+            }
+        }
+        else
+        {
+            newAddressing.append(faceI);
+            newFlipMap.append(fSet.flipMap()[i]);
+        }
+    }
+
+    if (nConflict > 0)
+    {
+        WarningIn("faceZoneSet::addSet(const topoSet&)")
+            << "addSet : there are " << nConflict
+            << " faces with different orientation in faceZonesSets "
+            << name() << " and " << set.name() << endl;
+    }
+
+    addressing_.transfer(newAddressing);
+    flipMap_.transfer(newFlipMap);
+    updateSet();
+}
+
+
+void faceZoneSet::deleteSet(const topoSet& set)
+{
+    label nConflict = 0;
+
+    DynamicList<label> newAddressing(addressing_.size());
+    DynamicList<bool> newFlipMap(flipMap_.size());
+
+    const faceZoneSet& fSet = refCast<const faceZoneSet>(set);
+
+    Map<label> faceToIndex(fSet.addressing().size());
+    forAll(fSet.addressing(), i)
+    {
+        faceToIndex.insert(fSet.addressing()[i], i);
+    }
+
+    forAll(addressing_, i)
+    {
+        label faceI = addressing_[i];
+
+        Map<label>::const_iterator iter = faceToIndex.find(faceI);
+
+        if (iter != faceToIndex.end())
+        {
+            label index = iter();
+
+            if (fSet.flipMap()[index] != flipMap_[i])
+            {
+                nConflict++;
+            }
+        }
+        else
+        {
+            // Not found in fSet so add
+            newAddressing.append(faceI);
+            newFlipMap.append(fSet.flipMap()[i]);
+        }
+    }
+
+    if (nConflict > 0)
+    {
+        WarningIn("faceZoneSet::deleteSet(const topoSet&)")
+            << "deleteSet : there are " << nConflict
+            << " faces with different orientation in faceZonesSets "
+            << name() << " and " << set.name() << endl;
+    }
+
+    addressing_.transfer(newAddressing);
+    flipMap_.transfer(newFlipMap);
+    updateSet();
+}
+
+
+void faceZoneSet::sync(const polyMesh& mesh)
+{}
+
+
+label faceZoneSet::maxSize(const polyMesh& mesh) const
+{
+    return mesh.nFaces();
+}
+
+
+//- Write using given format, version and compression
+bool faceZoneSet::writeObject
+(
+    IOstream::streamFormat s,
+    IOstream::versionNumber v,
+    IOstream::compressionType c
+) const
+{
+    // Write shadow faceSet
+    word oldTypeName = typeName;
+    const_cast<word&>(type()) = faceSet::typeName;
+    bool ok = faceSet::writeObject(s, v, c);
+    const_cast<word&>(type()) = oldTypeName;
+
+    // Modify faceZone
+    faceZoneMesh& faceZones = const_cast<polyMesh&>(mesh_).faceZones();
+    label zoneID = faceZones.findZoneID(name());
+
+    if (zoneID == -1)
+    {
+        zoneID = faceZones.size();
+
+        faceZones.setSize(zoneID+1);
+        faceZones.set
+        (
+            zoneID,
+            new faceZone
+            (
+                name(),
+                addressing_,
+                flipMap_,
+                zoneID,
+                faceZones
+            )
+        );
+    }
+    else
+    {
+        faceZones[zoneID].resetAddressing(addressing_, flipMap_);
+    }
+    faceZones.clearAddressing();
+
+    return ok && faceZones.write();
+}
+
+
+void faceZoneSet::updateMesh(const mapPolyMesh& morphMap)
+{
+    // faceZone
+    labelList newAddressing(addressing_.size());
+    boolList newFlipMap(flipMap_.size());
+
+    label n = 0;
+    forAll(addressing_, i)
+    {
+        label faceI = addressing_[i];
+        label newFaceI = morphMap.reverseFaceMap()[faceI];
+        if (newFaceI >= 0)
+        {
+            newAddressing[n] = newFaceI;
+            newFlipMap[n] = flipMap_[i];
+            n++;
+        }
+    }
+    newAddressing.setSize(n);
+    newFlipMap.setSize(n);
+
+    addressing_.transfer(newAddressing);
+    flipMap_.transfer(newFlipMap);
+
+    updateSet();
+}
+
+
+void faceZoneSet::writeDebug
+(
+    Ostream& os,
+    const primitiveMesh& mesh,
+    const label maxLen
+) const
+{
+    faceSet::writeDebug(os, mesh, maxLen);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/topoSets/faceZoneSet.H b/src/meshTools/sets/topoSets/faceZoneSet.H
new file mode 100644
index 0000000000000000000000000000000000000000..335bd9edead8f913a2593f975d7becb2d4a2a53e
--- /dev/null
+++ b/src/meshTools/sets/topoSets/faceZoneSet.H
@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::faceZoneSet
+
+Description
+    Like faceSet but updates faceZone when writing.
+
+SourceFiles
+    faceZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faceZoneSet_H
+#define faceZoneSet_H
+
+#include "faceSet.H"
+#include "boolList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class faceZoneSet Declaration
+\*---------------------------------------------------------------------------*/
+
+class faceZoneSet
+:
+    public faceSet
+{
+    // Private data
+
+        const polyMesh& mesh_;
+
+        labelList addressing_;
+
+        boolList flipMap_;
+
+   // Private Member Functions
+
+
+public:
+
+    //- Runtime type information
+    TypeName("faceZoneSet");
+
+
+    // Constructors
+
+        //- Construct from objectRegistry and name
+        faceZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            readOption r=MUST_READ,
+            writeOption w=NO_WRITE
+        );
+
+        //- Construct from additional size of labelHashSet
+        faceZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const label,
+            writeOption w=NO_WRITE
+        );
+
+        //- Construct from existing set
+        faceZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const topoSet&,
+            writeOption w=NO_WRITE
+        );
+
+
+
+    // Destructor
+
+        virtual ~faceZoneSet();
+
+
+    // Member functions
+
+        const labelList& addressing() const
+        {
+            return addressing_;
+        }
+
+        labelList& addressing()
+        {
+            return addressing_;
+        }
+
+
+        const boolList& flipMap() const
+        {
+            return flipMap_;
+        }
+
+        boolList& flipMap()
+        {
+            return flipMap_;
+        }
+
+
+        //- Sort addressing and make faceSet part consistent with addressing
+        void updateSet();
+
+        //- Invert contents. (insert all members 0..maxLen-1 which were not in
+        //  set)
+        virtual void invert(const label maxLen);
+
+        //- Subset contents. Only elements present in both sets remain.
+        virtual void subset(const topoSet& set);
+
+        //- Add elements present in set.
+        virtual void addSet(const topoSet& set);
+
+        //- Delete elements present in set.
+        virtual void deleteSet(const topoSet& set);
+
+        //- Sync faceZoneSet across coupled patches.
+        virtual void sync(const polyMesh& mesh);
+
+        //- Write maxLen items with label and coordinates. 
+        virtual void writeDebug
+        (
+            Ostream& os,
+            const primitiveMesh&,
+            const label maxLen
+        ) const;
+
+        //- Write faceZone
+        virtual bool writeObject
+        (
+            IOstream::streamFormat,
+            IOstream::versionNumber,
+            IOstream::compressionType
+        ) const;
+
+        //- Update any stored data for new labels
+        virtual void updateMesh(const mapPolyMesh& morphMap);
+
+        //- Return max index+1.
+        virtual label maxSize(const polyMesh& mesh) const;
+
+
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/topoSets/pointSet.C b/src/meshTools/sets/topoSets/pointSet.C
index 1e4884450cb88b1db2979bcc5c3ed3443cfbeb51..e370299fa1abf061bf475d5c7555042ec15bd24a 100644
--- a/src/meshTools/sets/topoSets/pointSet.C
+++ b/src/meshTools/sets/topoSets/pointSet.C
@@ -42,6 +42,7 @@ defineTypeNameAndDebug(pointSet, 0);
 
 addToRunTimeSelectionTable(topoSet, pointSet, word);
 addToRunTimeSelectionTable(topoSet, pointSet, size);
+addToRunTimeSelectionTable(topoSet, pointSet, set);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -78,6 +79,18 @@ pointSet::pointSet
 {}
 
 
+pointSet::pointSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const topoSet& set,
+    writeOption w
+)
+:
+    topoSet(mesh, name, set, w)
+{}
+
+
 pointSet::pointSet
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/topoSets/pointSet.H b/src/meshTools/sets/topoSets/pointSet.H
index b9c8ef93aceb6a2fcd7124f6da1cb2615d1348dc..05f4ef6e36b193ab16b8cbf26f46b0130ab1c1e7 100644
--- a/src/meshTools/sets/topoSets/pointSet.H
+++ b/src/meshTools/sets/topoSets/pointSet.H
@@ -83,6 +83,15 @@ public:
             writeOption w=NO_WRITE
         );
 
+        //- Construct from additional labelHashSet
+        pointSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const topoSet&,
+            writeOption w=NO_WRITE
+        );
+
         //- Construct from additional labelHashSet
         pointSet
         (
diff --git a/src/meshTools/sets/topoSets/pointZoneSet.C b/src/meshTools/sets/topoSets/pointZoneSet.C
new file mode 100644
index 0000000000000000000000000000000000000000..32f447d1447ed5723d589eb6ed642e9b73bc3a9d
--- /dev/null
+++ b/src/meshTools/sets/topoSets/pointZoneSet.C
@@ -0,0 +1,314 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "pointZoneSet.H"
+#include "mapPolyMesh.H"
+#include "polyMesh.H"
+#include "processorPolyPatch.H"
+#include "cyclicPolyPatch.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(pointZoneSet, 0);
+
+addToRunTimeSelectionTable(topoSet, pointZoneSet, word);
+addToRunTimeSelectionTable(topoSet, pointZoneSet, size);
+addToRunTimeSelectionTable(topoSet, pointZoneSet, set);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void pointZoneSet::updateSet()
+{
+    labelList order;
+    sortedOrder(addressing_, order);
+    inplaceReorder(order, addressing_);
+
+    pointSet::clearStorage();
+    pointSet::resize(2*addressing_.size());
+    forAll(addressing_, i)
+    {
+        pointSet::insert(addressing_[i]);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+pointZoneSet::pointZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    readOption r,
+    writeOption w
+)
+:
+    pointSet(mesh, name, 1000),  // do not read pointSet
+    mesh_(mesh),
+    addressing_(0)
+{
+    const pointZoneMesh& pointZones = mesh.pointZones();
+    label zoneID = pointZones.findZoneID(name);
+
+    if
+    (
+        (r == IOobject::MUST_READ)
+     || (r == IOobject::READ_IF_PRESENT && zoneID != -1)
+    )
+    {
+        const pointZone& fz = pointZones[zoneID];
+        addressing_ = fz;
+    }
+
+    updateSet();
+
+    check(mesh.nPoints());
+}
+
+
+pointZoneSet::pointZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const label size,
+    writeOption w
+)
+:
+    pointSet(mesh, name, size, w),
+    mesh_(mesh),
+    addressing_(0)
+{
+    updateSet();
+}
+
+
+pointZoneSet::pointZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const topoSet& set,
+    writeOption w
+)
+:
+    pointSet(mesh, name, set.size(), w),
+    mesh_(mesh),
+    addressing_(refCast<const pointZoneSet>(set).addressing())
+{
+    updateSet();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+pointZoneSet::~pointZoneSet()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void pointZoneSet::invert(const label maxLen)
+{
+    label n = 0;
+
+    for (label pointI = 0; pointI < maxLen; pointI++)
+    {
+        if (!found(pointI))
+        {
+            addressing_[n] = pointI;
+            n++;
+        }
+    }
+    addressing_.setSize(n);
+    updateSet();
+}
+
+
+void pointZoneSet::subset(const topoSet& set)
+{
+    DynamicList<label> newAddressing(addressing_.size());
+
+    const pointZoneSet& fSet = refCast<const pointZoneSet>(set);
+
+    forAll(fSet.addressing(), i)
+    {
+        label pointI = fSet.addressing()[i];
+
+        if (found(pointI))
+        {
+            newAddressing.append(pointI);
+        }
+    }
+
+    addressing_.transfer(newAddressing);
+    updateSet();
+}
+
+
+void pointZoneSet::addSet(const topoSet& set)
+{
+    DynamicList<label> newAddressing(addressing_);
+
+    const pointZoneSet& fSet = refCast<const pointZoneSet>(set);
+
+    forAll(fSet.addressing(), i)
+    {
+        label pointI = fSet.addressing()[i];
+
+        if (!found(pointI))
+        {
+            newAddressing.append(pointI);
+        }
+    }
+
+    addressing_.transfer(newAddressing);
+    updateSet();
+}
+
+
+void pointZoneSet::deleteSet(const topoSet& set)
+{
+    DynamicList<label> newAddressing(addressing_.size());
+
+    const pointZoneSet& fSet = refCast<const pointZoneSet>(set);
+
+    forAll(addressing_, i)
+    {
+        label pointI = addressing_[i];
+
+        if (!fSet.found(pointI))
+        {
+            // Not found in fSet so add
+            newAddressing.append(pointI);
+        }
+    }
+
+    addressing_.transfer(newAddressing);
+    updateSet();
+}
+
+
+void pointZoneSet::sync(const polyMesh& mesh)
+{}
+
+
+label pointZoneSet::maxSize(const polyMesh& mesh) const
+{
+    return mesh.nPoints();
+}
+
+
+//- Write using given format, version and compression
+bool pointZoneSet::writeObject
+(
+    IOstream::streamFormat s,
+    IOstream::versionNumber v,
+    IOstream::compressionType c
+) const
+{
+    // Write shadow pointSet
+    word oldTypeName = typeName;
+    const_cast<word&>(type()) = pointSet::typeName;
+    bool ok = pointSet::writeObject(s, v, c);
+    const_cast<word&>(type()) = oldTypeName;
+
+    // Modify pointZone
+    pointZoneMesh& pointZones = const_cast<polyMesh&>(mesh_).pointZones();
+    label zoneID = pointZones.findZoneID(name());
+
+    if (zoneID == -1)
+    {
+        zoneID = pointZones.size();
+
+        pointZones.setSize(zoneID+1);
+        pointZones.set
+        (
+            zoneID,
+            new pointZone
+            (
+                name(),
+                addressing_,
+                zoneID,
+                pointZones
+            )
+        );
+    }
+    else
+    {
+        pointZones[zoneID] = addressing_;
+    }
+    pointZones.clearAddressing();
+
+    return ok && pointZones.write();
+}
+
+
+void pointZoneSet::updateMesh(const mapPolyMesh& morphMap)
+{
+    // pointZone
+    labelList newAddressing(addressing_.size());
+
+    label n = 0;
+    forAll(addressing_, i)
+    {
+        label pointI = addressing_[i];
+        label newPointI = morphMap.reversePointMap()[pointI];
+        if (newPointI >= 0)
+        {
+            newAddressing[n] = newPointI;
+            n++;
+        }
+    }
+    newAddressing.setSize(n);
+
+    addressing_.transfer(newAddressing);
+
+    updateSet();
+}
+
+
+void pointZoneSet::writeDebug
+(
+    Ostream& os,
+    const primitiveMesh& mesh,
+    const label maxLen
+) const
+{
+    pointSet::writeDebug(os, mesh, maxLen);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/topoSets/pointZoneSet.H b/src/meshTools/sets/topoSets/pointZoneSet.H
new file mode 100644
index 0000000000000000000000000000000000000000..7f93886e0bc71b3b8eddc31231f262df5ab6fe29
--- /dev/null
+++ b/src/meshTools/sets/topoSets/pointZoneSet.H
@@ -0,0 +1,173 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::pointZoneSet
+
+Description
+    Like pointSet but updates pointZone when writing.
+
+SourceFiles
+    pointZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pointZoneSet_H
+#define pointZoneSet_H
+
+#include "pointSet.H"
+#include "boolList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class pointZoneSet Declaration
+\*---------------------------------------------------------------------------*/
+
+class pointZoneSet
+:
+    public pointSet
+{
+    // Private data
+
+        const polyMesh& mesh_;
+
+        labelList addressing_;
+
+   // Private Member Functions
+
+
+public:
+
+    //- Runtime type information
+    TypeName("pointZoneSet");
+
+
+    // Constructors
+
+        //- Construct from objectRegistry and name
+        pointZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            readOption r=MUST_READ,
+            writeOption w=NO_WRITE
+        );
+
+        //- Construct from additional size of labelHashSet
+        pointZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const label,
+            writeOption w=NO_WRITE
+        );
+
+        //- Construct from existing set
+        pointZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const topoSet&,
+            writeOption w=NO_WRITE
+        );
+
+
+
+    // Destructor
+
+        virtual ~pointZoneSet();
+
+
+    // Member functions
+
+        const labelList& addressing() const
+        {
+            return addressing_;
+        }
+
+        labelList& addressing()
+        {
+            return addressing_;
+        }
+
+        //- Sort addressing and make pointSet part consistent with addressing
+        void updateSet();
+
+        //- Invert contents. (insert all members 0..maxLen-1 which were not in
+        //  set)
+        virtual void invert(const label maxLen);
+
+        //- Subset contents. Only elements present in both sets remain.
+        virtual void subset(const topoSet& set);
+
+        //- Add elements present in set.
+        virtual void addSet(const topoSet& set);
+
+        //- Delete elements present in set.
+        virtual void deleteSet(const topoSet& set);
+
+        //- Sync pointZoneSet across coupled patches.
+        virtual void sync(const polyMesh& mesh);
+
+        //- Write maxLen items with label and coordinates. 
+        virtual void writeDebug
+        (
+            Ostream& os,
+            const primitiveMesh&,
+            const label maxLen
+        ) const;
+
+        //- Write pointZone
+        virtual bool writeObject
+        (
+            IOstream::streamFormat,
+            IOstream::versionNumber,
+            IOstream::compressionType
+        ) const;
+
+        //- Update any stored data for new labels
+        virtual void updateMesh(const mapPolyMesh& morphMap);
+
+        //- Return max index+1.
+        virtual label maxSize(const polyMesh& mesh) const;
+
+
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/topoSets/topoSet.C b/src/meshTools/sets/topoSets/topoSet.C
index 32f8d1ff989c6981b216fa5e83bb6fe23f269661..28bc14f469fe1c0facf4ec6b48aae8b6ee0d7eb5 100644
--- a/src/meshTools/sets/topoSets/topoSet.C
+++ b/src/meshTools/sets/topoSets/topoSet.C
@@ -39,6 +39,7 @@ namespace Foam
 defineTypeNameAndDebug(topoSet, 0);
 defineRunTimeSelectionTable(topoSet, word);
 defineRunTimeSelectionTable(topoSet, size);
+defineRunTimeSelectionTable(topoSet, set);
 
 
 // Construct named object from existing set.
@@ -103,6 +104,37 @@ autoPtr<topoSet> topoSet::New
 }
 
 
+// Construct named object from existing set.
+autoPtr<topoSet> topoSet::New
+(
+    const word& setType,
+    const polyMesh& mesh,
+    const word& name,
+    const topoSet& set,
+    writeOption w
+)
+{
+    setConstructorTable::iterator cstrIter =
+        setConstructorTablePtr_
+            ->find(setType);
+
+    if (cstrIter == setConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "topoSet::New(const word&, "
+            "const polyMesh&, const word&, const topoSet&, writeOption)"
+        )   << "Unknown set type " << setType
+            << endl << endl
+            << "Valid set types : " << endl
+            << setConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<topoSet>(cstrIter()(mesh, name, set, w));
+}
+
+
 Foam::fileName topoSet::topoSet::localPath
 (
     const polyMesh& mesh,
@@ -255,9 +287,9 @@ void topoSet::writeDebug
 ) const
 {
     // Bounding box of contents.
-    boundBox bb(pointField(coords, toc()));
+    boundBox bb(pointField(coords, toc()), true);
 
-    Pout<< "Set bounding box: min = "
+    os  << "Set bounding box: min = "
         << bb.min() << "    max = " << bb.max() << " meters. " << endl << endl;
 
     label n = 0;
@@ -538,18 +570,18 @@ void topoSet::writeDebug(Ostream& os, const label maxLen) const
 }
 
 
-void topoSet::writeDebug
-(
-    Ostream&,
-    const primitiveMesh&,
-    const label
-) const
-{
-    notImplemented
-    (
-        "topoSet::writeDebug(Ostream&, const primitiveMesh&, const label)"
-    );
-}
+//void topoSet::writeDebug
+//(
+//    Ostream&,
+//    const primitiveMesh&,
+//    const label
+//) const
+//{
+//    notImplemented
+//    (
+//        "topoSet::writeDebug(Ostream&, const primitiveMesh&, const label)"
+//    );
+//}
 
 
 bool topoSet::writeData(Ostream& os) const
@@ -564,13 +596,13 @@ void topoSet::updateMesh(const mapPolyMesh&)
 }
 
 
-//- Return max index+1.
-label topoSet::maxSize(const polyMesh&) const
-{
-    notImplemented("topoSet::maxSize(const polyMesh&)");
-
-    return -1;
-}
+////- Return max index+1.
+//label topoSet::maxSize(const polyMesh&) const
+//{
+//    notImplemented("topoSet::maxSize(const polyMesh&)");
+//
+//    return -1;
+//}
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
diff --git a/src/meshTools/sets/topoSets/topoSet.H b/src/meshTools/sets/topoSets/topoSet.H
index 51049ab9e25a3309dcfcaeb555355c59e39f1daa..124681305a6c93abf1fb25471738dce47ae061f0 100644
--- a/src/meshTools/sets/topoSets/topoSet.H
+++ b/src/meshTools/sets/topoSets/topoSet.H
@@ -55,7 +55,6 @@ namespace Foam
 class mapPolyMesh;
 class polyMesh;
 class primitiveMesh;
-//class directPolyTopoChange;
 
 /*---------------------------------------------------------------------------*\
                            Class topoSet Declaration
@@ -154,6 +153,21 @@ public:
             (mesh, name, size, w)
         );
 
+        // For the constructor as copy
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            topoSet,
+            set,
+            (
+                const polyMesh& mesh,
+                const word& name,
+                const topoSet& set,
+                writeOption w
+            ),
+            (mesh, name, set, w)
+        );
+
 
     // Constructors
 
@@ -208,7 +222,7 @@ public:
 
     // Selectors
 
-        //- Return a reference to the selected source
+        //- Return a pointer to a toposet read from file
         static autoPtr<topoSet> New
         (
             const word& setType,
@@ -218,7 +232,7 @@ public:
             writeOption w=NO_WRITE
         );
 
-        //- Return a reference to the selected source
+        //- Return a pointer to a new toposet of given size
         static autoPtr<topoSet> New
         (
             const word& setType,
@@ -228,6 +242,16 @@ public:
             writeOption w=NO_WRITE
         );
 
+        //- Return a pointer to a new toposet as copy of another toposet
+        static autoPtr<topoSet> New
+        (
+            const word& setType,
+            const polyMesh& mesh,
+            const word& name,
+            const topoSet& set,
+            writeOption w=NO_WRITE
+        );
+
 
     // Destructor
 
@@ -256,13 +280,13 @@ public:
         virtual void writeDebug(Ostream& os, const label maxLen) const;
 
         //- Like above but also writes mesh related quantity
-        //  (usually coordinate). Not implemented.
+        //  (usually coordinate).
         virtual void writeDebug
         (
             Ostream& os,
             const primitiveMesh&,
             const label maxLen
-        ) const;
+        ) const = 0;
 
         //- Write contents.
         virtual bool writeData(Ostream&) const;
@@ -271,7 +295,7 @@ public:
         virtual void updateMesh(const mapPolyMesh& morphMap);
 
         //- Return max allowable index (+1). Not implemented.
-        virtual label maxSize(const polyMesh& mesh) const;
+        virtual label maxSize(const polyMesh& mesh) const = 0;
 
 
 
diff --git a/tutorials/incompressible/pimpleDyMFoam/movingCone/0/cellMotionUx b/tutorials/incompressible/pimpleDyMFoam/movingCone/0/cellMotionUx
deleted file mode 100644
index 5ff626074066c0db5a318340109c0e95272165e4..0000000000000000000000000000000000000000
--- a/tutorials/incompressible/pimpleDyMFoam/movingCone/0/cellMotionUx
+++ /dev/null
@@ -1,60 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    object      motionU;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 1 -1 0 0 0 0];
-
-internalField   uniform 0;
-
-boundaryField
-{
-    movingWall
-    {
-        type            fixedValue;
-        value           uniform 1;
-    }
-    farFieldMoving
-    {
-        type            slip;
-    }
-    fixedWall
-    {
-        type            fixedValue;
-        value           uniform 0;
-    }
-    axis
-    {
-        type            symmetryPlane;
-    }
-    left
-    {
-        type            fixedValue;
-        value           uniform 0;
-    }
-    farField
-    {
-        type            slip;
-    }
-    back
-    {
-        type            wedge;
-    }
-    front
-    {
-        type            wedge;
-    }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleDyMFoam/movingCone/0/pointMotionUx b/tutorials/incompressible/pimpleDyMFoam/movingCone/0/pointMotionUx
index e0a68b5aff1f4d4d3c3ccd6ed8238c2e8e4b339c..ebad161cb509d65fc102fe1504e7fea9d927c5b8 100644
--- a/tutorials/incompressible/pimpleDyMFoam/movingCone/0/pointMotionUx
+++ b/tutorials/incompressible/pimpleDyMFoam/movingCone/0/pointMotionUx
@@ -27,7 +27,33 @@ boundaryField
     }
     farFieldMoving
     {
-        type            slip;
+        //type            slip;
+        type            surfaceSlipDisplacement;
+        geometry
+        {
+            top
+            {
+                type        searchablePlane;
+                planeType   pointAndNormal;
+                basePoint   (0 0.0025 0);
+               normalVector (0 1 0); 
+            }
+        };
+
+        // Find projection with surface:
+        //     fixedNormal : intersections along prespecified direction
+        //     pointNormal : intersections along current pointNormal of patch
+        //     nearest     : nearest point on surface
+        followMode fixedNormal;
+
+        // if fixedNormal : normal
+        projectDirection (0 1 0);
+
+        //- -1 or component to knock out before doing projection
+        //wedgePlane      -1;
+
+        //- Points that should remain fixed
+        //frozenPointsZone fixedPointsZone;
     }
     fixedWall
     {
@@ -45,7 +71,28 @@ boundaryField
     }
     farField
     {
-        type            slip;
+        //type            slip;
+        type            surfaceSlipDisplacement;
+        geometry
+        {
+            top
+            {
+                type        searchablePlane;
+                planeType   pointAndNormal;
+                basePoint   (0 0.0025 0);
+               normalVector (0 1 0); 
+            }
+        };
+
+        followMode fixedNormal;
+
+        projectDirection (0 1 0);
+
+        //- -1 or component to knock out before doing projection
+        wedgePlane      -1;
+
+        //- Points that should remain fixed
+        //frozenPointsZone fixedPointsZone;
     }
     back
     {
diff --git a/tutorials/incompressible/pimpleDyMFoam/movingCone/system/controlDict b/tutorials/incompressible/pimpleDyMFoam/movingCone/system/controlDict
index fac973c9f9e5336e3df6a40df659c4206c6a5131..ae0f289880a4414148128e12161575a585486fd8 100644
--- a/tutorials/incompressible/pimpleDyMFoam/movingCone/system/controlDict
+++ b/tutorials/incompressible/pimpleDyMFoam/movingCone/system/controlDict
@@ -17,6 +17,9 @@ FoamFile
 
 application     pimpleDyMFoam;
 
+// For surfaceSlip boundary conditions
+libs            ("libfvMotionSolvers.so");
+
 startFrom       startTime;
 
 startTime       0;