diff --git a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
index 4abdc6c607..b57e181a03 100644
--- a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
+++ b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
@@ -101,296 +101,13 @@ Description
 #include "syncTools.H"
 #include "ReadFields.H"
 #include "mappedWallPolyPatch.H"
+#include "fvMeshTools.H"
 #include "zeroGradientFvPatchFields.H"
 
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template<class GeoField>
-void addPatchFields(fvMesh& mesh, const word& patchFieldType)
-{
-    HashTable<const GeoField*> flds
-    (
-        mesh.objectRegistry::lookupClass<GeoField>()
-    );
-
-    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
-    {
-        const GeoField& fld = *iter();
-
-        typename GeoField::GeometricBoundaryField& bfld =
-            const_cast<typename GeoField::GeometricBoundaryField&>
-            (
-                fld.boundaryField()
-            );
-
-        label sz = bfld.size();
-        bfld.setSize(sz+1);
-        bfld.set
-        (
-            sz,
-            GeoField::PatchFieldType::New
-            (
-                patchFieldType,
-                mesh.boundary()[sz],
-                fld.dimensionedInternalField()
-            )
-        );
-    }
-}
-
-
-// Remove last patch field
-template<class GeoField>
-void trimPatchFields(fvMesh& mesh, const label nPatches)
-{
-    HashTable<const GeoField*> flds
-    (
-        mesh.objectRegistry::lookupClass<GeoField>()
-    );
-
-    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
-    {
-        const GeoField& fld = *iter();
-
-        const_cast<typename GeoField::GeometricBoundaryField&>
-        (
-            fld.boundaryField()
-        ).setSize(nPatches);
-    }
-}
-
-
-// Reorder patch field
-template<class GeoField>
-void reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
-{
-    HashTable<const GeoField*> flds
-    (
-        mesh.objectRegistry::lookupClass<GeoField>()
-    );
-
-    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
-    {
-        const GeoField& fld = *iter();
-
-        typename GeoField::GeometricBoundaryField& bfld =
-            const_cast<typename GeoField::GeometricBoundaryField&>
-            (
-                fld.boundaryField()
-            );
-
-        bfld.reorder(oldToNew);
-    }
-}
-
-
-// Adds patch if not yet there. Returns patchID.
-label addPatch(fvMesh& mesh, const polyPatch& patch)
-{
-    polyBoundaryMesh& polyPatches =
-        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
-
-    label patchI = polyPatches.findPatchID(patch.name());
-    if (patchI != -1)
-    {
-        if (polyPatches[patchI].type() == patch.type())
-        {
-            // Already there
-            return patchI;
-        }
-        else
-        {
-            FatalErrorIn("addPatch(fvMesh&, const polyPatch*)")
-                << "Already have patch " << patch.name()
-                << " but of type " << patch.type()
-                << exit(FatalError);
-        }
-    }
-
-
-    label insertPatchI = polyPatches.size();
-    label startFaceI = mesh.nFaces();
-
-    forAll(polyPatches, patchI)
-    {
-        const polyPatch& pp = polyPatches[patchI];
-
-        if (isA<processorPolyPatch>(pp))
-        {
-            insertPatchI = patchI;
-            startFaceI = pp.start();
-            break;
-        }
-    }
-
-
-    // Below is all quite a hack. Feel free to change once there is a better
-    // mechanism to insert and reorder patches.
-
-    // Clear local fields and e.g. polyMesh parallelInfo.
-    mesh.clearOut();
-
-    label sz = polyPatches.size();
-
-    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
-
-    // Add polyPatch at the end
-    polyPatches.setSize(sz+1);
-    polyPatches.set
-    (
-        sz,
-        patch.clone
-        (
-            polyPatches,
-            insertPatchI,   //index
-            0,              //size
-            startFaceI      //start
-        )
-    );
-    fvPatches.setSize(sz+1);
-    fvPatches.set
-    (
-        sz,
-        fvPatch::New
-        (
-            polyPatches[sz],  // point to newly added polyPatch
-            mesh.boundary()
-        )
-    );
-
-    addPatchFields<volScalarField>
-    (
-        mesh,
-        calculatedFvPatchField<scalar>::typeName
-    );
-    addPatchFields<volVectorField>
-    (
-        mesh,
-        calculatedFvPatchField<vector>::typeName
-    );
-    addPatchFields<volSphericalTensorField>
-    (
-        mesh,
-        calculatedFvPatchField<sphericalTensor>::typeName
-    );
-    addPatchFields<volSymmTensorField>
-    (
-        mesh,
-        calculatedFvPatchField<symmTensor>::typeName
-    );
-    addPatchFields<volTensorField>
-    (
-        mesh,
-        calculatedFvPatchField<tensor>::typeName
-    );
-
-    // Surface fields
-
-    addPatchFields<surfaceScalarField>
-    (
-        mesh,
-        calculatedFvPatchField<scalar>::typeName
-    );
-    addPatchFields<surfaceVectorField>
-    (
-        mesh,
-        calculatedFvPatchField<vector>::typeName
-    );
-    addPatchFields<surfaceSphericalTensorField>
-    (
-        mesh,
-        calculatedFvPatchField<sphericalTensor>::typeName
-    );
-    addPatchFields<surfaceSymmTensorField>
-    (
-        mesh,
-        calculatedFvPatchField<symmTensor>::typeName
-    );
-    addPatchFields<surfaceTensorField>
-    (
-        mesh,
-        calculatedFvPatchField<tensor>::typeName
-    );
-
-    // Create reordering list
-    // patches before insert position stay as is
-    labelList oldToNew(sz+1);
-    for (label i = 0; i < insertPatchI; i++)
-    {
-        oldToNew[i] = i;
-    }
-    // patches after insert position move one up
-    for (label i = insertPatchI; i < sz; i++)
-    {
-        oldToNew[i] = i+1;
-    }
-    // appended patch gets moved to insert position
-    oldToNew[sz] = insertPatchI;
-
-    // Shuffle into place
-    polyPatches.reorder(oldToNew);
-    fvPatches.reorder(oldToNew);
-
-    reorderPatchFields<volScalarField>(mesh, oldToNew);
-    reorderPatchFields<volVectorField>(mesh, oldToNew);
-    reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
-    reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
-    reorderPatchFields<volTensorField>(mesh, oldToNew);
-    reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
-    reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
-    reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
-    reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
-    reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
-
-    return insertPatchI;
-}
-
-
-// Reorder and delete patches.
-void reorderPatches
-(
-    fvMesh& mesh,
-    const labelList& oldToNew,
-    const label nNewPatches
-)
-{
-    polyBoundaryMesh& polyPatches =
-        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
-    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
-
-    // Shuffle into place
-    polyPatches.reorder(oldToNew);
-    fvPatches.reorder(oldToNew);
-
-    reorderPatchFields<volScalarField>(mesh, oldToNew);
-    reorderPatchFields<volVectorField>(mesh, oldToNew);
-    reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
-    reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
-    reorderPatchFields<volTensorField>(mesh, oldToNew);
-    reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
-    reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
-    reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
-    reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
-    reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
-
-    // Remove last.
-    polyPatches.setSize(nNewPatches);
-    fvPatches.setSize(nNewPatches);
-    trimPatchFields<volScalarField>(mesh, nNewPatches);
-    trimPatchFields<volVectorField>(mesh, nNewPatches);
-    trimPatchFields<volSphericalTensorField>(mesh, nNewPatches);
-    trimPatchFields<volSymmTensorField>(mesh, nNewPatches);
-    trimPatchFields<volTensorField>(mesh, nNewPatches);
-    trimPatchFields<surfaceScalarField>(mesh, nNewPatches);
-    trimPatchFields<surfaceVectorField>(mesh, nNewPatches);
-    trimPatchFields<surfaceSphericalTensorField>(mesh, nNewPatches);
-    trimPatchFields<surfaceSymmTensorField>(mesh, nNewPatches);
-    trimPatchFields<surfaceTensorField>(mesh, nNewPatches);
-}
-
-
 // Prepend prefix to selected patches.
 void renamePatches
 (
@@ -1193,8 +910,8 @@ void createAndWriteRegion
         }
     }
 
-    reorderPatches(newMesh(), oldToNew, nNewPatches);
-
+    //reorderPatches(newMesh(), oldToNew, nNewPatches);
+    fvMeshTools::reorderPatches(newMesh(), oldToNew, nNewPatches, true);
 
     // Rename shared patches with region name
     if (prefixRegion)
@@ -1360,7 +1077,15 @@ labelList addRegionPatches
             mesh.boundaryMesh()
         );
 
-        interfacePatches[interI] = addPatch(mesh, patch1);
+        //interfacePatches[interI] = addPatch(mesh, patch1);
+        interfacePatches[interI] = fvMeshTools::addPatch
+        (
+            mesh,
+            patch1,
+            dictionary(),   //optional per field value
+            calculatedFvPatchField<scalar>::typeName,
+            true            //validBoundary
+        );
 
         mappedWallPolyPatch patch2
         (
@@ -1374,7 +1099,15 @@ labelList addRegionPatches
             point::zero,        // offset
             mesh.boundaryMesh()
         );
-        addPatch(mesh, patch2);
+        //addPatch(mesh, patch2);
+        fvMeshTools::addPatch
+        (
+            mesh,
+            patch2,
+            dictionary(),   //optional per field value
+            calculatedFvPatchField<scalar>::typeName,
+            true            //validBoundary
+        );
 
         Info<< "For interface between region " << regionNames[e[0]]
             << " and " << regionNames[e[1]] << " added patches" << endl
diff --git a/applications/utilities/mesh/manipulation/topoSet/topoSetDict b/applications/utilities/mesh/manipulation/topoSet/topoSetDict
index 31fa348e97..00720d2382 100644
--- a/applications/utilities/mesh/manipulation/topoSet/topoSetDict
+++ b/applications/utilities/mesh/manipulation/topoSet/topoSetDict
@@ -226,6 +226,12 @@ FoamFile
 //                            // (regular expressions allowed)
 //    }
 //
+//    // All boundary faces
+//    source boundaryToFace;
+//    sourceInfo
+//    {
+//    }
+//
 //    // All faces of faceZone
 //    source zoneToFace;
 //    sourceInfo
@@ -359,6 +365,21 @@ FoamFile
 //        cellSet c0;       // name of cellSet of slave side
 //    }
 //
+//    // Select based on surface. Orientation from normals on surface
+//    {
+//        name    fz0;
+//        type    faceZoneSet;
+//        action  new;
+//        source  searchableSurfaceToFaceZone;
+//        sourceInfo
+//        {
+//            surface searchableSphere;
+//            centre  (0.05 0.05 0.005);
+//            radius  0.025;
+//        }
+//    }
+//
+//
 //
 // pointZoneSet
 // ~~~~~~~~~~~~
diff --git a/src/OpenFOAM/db/dictionary/dictionary.C b/src/OpenFOAM/db/dictionary/dictionary.C
index 295749ea4c..2a2ae85253 100644
--- a/src/OpenFOAM/db/dictionary/dictionary.C
+++ b/src/OpenFOAM/db/dictionary/dictionary.C
@@ -424,38 +424,87 @@ const Foam::entry* Foam::dictionary::lookupScopedEntryPtr
     bool patternMatch
 ) const
 {
-    string::size_type dotPos = keyword.find('.');
-
-    if (dotPos == string::npos)
+    if (keyword[0] == ':')
     {
-        return lookupEntryPtr(keyword, recursive, patternMatch);
+        // Go up to top level
+        const dictionary* dictPtr = this;
+        while (&dictPtr->parent_ != &dictionary::null)
+        {
+            dictPtr = &dictPtr->parent_;
+        }
+
+        // At top. Recurse to find entries
+        return dictPtr->lookupScopedEntryPtr
+        (
+            keyword.substr(1, keyword.size()-1),
+            false,
+            patternMatch
+        );
     }
     else
     {
-        if (dotPos == 0)
-        {
-            const dictionary* dictPtr = this;
-            while (&dictPtr->parent_ != &dictionary::null)
-            {
-                dictPtr = &dictPtr->parent_;
-            }
+        string::size_type dotPos = keyword.find('.');
 
-            // At top
-            return dictPtr->lookupScopedEntryPtr
-            (
-                keyword.substr(1, keyword.size()-1),
-                false,
-                patternMatch
-            );
+        if (dotPos == string::npos)
+        {
+            // Non-scoped lookup
+            return lookupEntryPtr(keyword, recursive, patternMatch);
         }
         else
         {
-            wordList entryNames(fileName(keyword).components('.'));
-
-            const entry* entPtr = lookupEntryPtr(entryNames[0], false, true);
-
-            for (int i=1; i<entryNames.size(); ++i)
+            if (dotPos == 0)
             {
+                // Starting with a '.'. Go up for every 2nd '.' found
+
+                const dictionary* dictPtr = this;
+
+                string::size_type begVar = dotPos + 1;
+                string::const_iterator iter = keyword.begin() + begVar;
+                string::size_type endVar = begVar;
+                while
+                (
+                    iter != keyword.end()
+                 && *iter == '.'
+                )
+                {
+                    ++iter;
+                    ++endVar;
+
+                    // Go to parent
+                    if (&dictPtr->parent_ == &dictionary::null)
+                    {
+                        FatalIOErrorIn
+                        (
+                            "dictionary::lookupScopedEntryPtr"
+                            "(const word&, bool, bool)",
+                            *this
+                        )   << "No parent of current dictionary"
+                            << " when searching for "
+                            << keyword.substr(begVar, keyword.size()-begVar)
+                            << exit(FatalIOError);
+                    }
+                    dictPtr = &dictPtr->parent_;
+                }
+
+                return dictPtr->lookupScopedEntryPtr
+                (
+                    keyword.substr(endVar),
+                    false,
+                    patternMatch
+                );
+            }
+            else
+            {
+                // Extract the first word
+                word firstWord = keyword.substr(0, dotPos);
+
+                const entry* entPtr = lookupScopedEntryPtr
+                (
+                    firstWord,
+                    false,          //recursive
+                    patternMatch
+                );
+
                 if (!entPtr)
                 {
                     FatalIOErrorIn
@@ -463,46 +512,27 @@ const Foam::entry* Foam::dictionary::lookupScopedEntryPtr
                         "dictionary::lookupScopedEntryPtr"
                         "(const word&, bool, bool)",
                         *this
-                    )   << "keyword " << keyword
+                    )   << "keyword " << firstWord
                         << " is undefined in dictionary "
                         << name() << endl
                         << "Valid keywords are " << keys()
                         << exit(FatalIOError);
                 }
-                if (!entPtr->isDict())
+
+                if (entPtr->isDict())
                 {
-                    FatalIOErrorIn
+                    return entPtr->dict().lookupScopedEntryPtr
                     (
-                        "dictionary::lookupScopedEntryPtr"
-                        "(const word&, bool, bool)",
-                        *this
-                    )   << "Entry " << entPtr->name()
-                        << " is not a dictionary so cannot lookup sub entry "
-                        << entryNames[i]
-                        << exit(FatalIOError);
+                        keyword.substr(dotPos, keyword.size()-dotPos),
+                        false,
+                        patternMatch
+                    );
+                }
+                else
+                {
+                    return NULL;
                 }
-
-                entPtr = entPtr->dict().lookupEntryPtr
-                (
-                    entryNames[i],
-                    false,
-                    true
-                );
             }
-
-            if (!entPtr)
-            {
-                FatalIOErrorIn
-                (
-                    "dictionary::lookupScopedEntryPtr"
-                    "(const word&, bool, bool)",
-                    *this
-                )   << "keyword " << keyword
-                    << " is not a valid scoped entry in dictionary "
-                    << name()
-                    << exit(FatalIOError);
-            }
-            return entPtr;
         }
     }
 }
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
index a54438f105..b4148405cf 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
@@ -1078,7 +1078,11 @@ void Foam::polyBoundaryMesh::updateMesh()
 }
 
 
-void Foam::polyBoundaryMesh::reorder(const labelUList& oldToNew)
+void Foam::polyBoundaryMesh::reorder
+(
+    const labelUList& oldToNew,
+    const bool validBoundary
+)
 {
     // Change order of patches
     polyPatchList::reorder(oldToNew);
@@ -1091,7 +1095,10 @@ void Foam::polyBoundaryMesh::reorder(const labelUList& oldToNew)
         patches[patchI].index() = patchI;
     }
 
-    updateMesh();
+    if (validBoundary)
+    {
+        updateMesh();
+    }
 }
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
index 34f0a6d409..ad937a86fd 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
@@ -219,9 +219,10 @@ public:
 
         //- Reorders patches. Ordering does not have to be done in
         //  ascending or descending order. Reordering has to be unique.
-        //  (is shuffle) Calls updateMesh() after reordering to recalculate
-        //  data.
-        void reorder(const labelUList&);
+        //  (is shuffle) If validBoundary calls updateMesh()
+        //  after reordering to recalculate data (so call needs to be parallel
+        //  sync in that case)
+        void reorder(const labelUList&, const bool validBoundary);
 
         //- writeData member function required by regIOobject
         bool writeData(Ostream&) const;
diff --git a/src/dynamicMesh/Make/files b/src/dynamicMesh/Make/files
index 21dcc72398..bb0ec3ddef 100644
--- a/src/dynamicMesh/Make/files
+++ b/src/dynamicMesh/Make/files
@@ -81,6 +81,8 @@ fvMeshDistribute/fvMeshDistribute.C
 polyMeshAdder/faceCoupleInfo.C
 polyMeshAdder/polyMeshAdder.C
 
+fvMeshTools/fvMeshTools.C
+
 motionSmoother/motionSmoother.C
 motionSmoother/motionSmootherCheck.C
 motionSmoother/polyMeshGeometry/polyMeshGeometry.C
diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
index 33058afaf4..7fe993522a 100644
--- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
+++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
@@ -31,7 +31,6 @@ License
 #include "processorFvsPatchField.H"
 #include "processorCyclicPolyPatch.H"
 #include "processorCyclicFvPatchField.H"
-#include "processorCyclicFvsPatchField.H"
 #include "polyTopoChange.H"
 #include "removeCells.H"
 #include "polyModifyFace.H"
@@ -40,6 +39,7 @@ License
 #include "surfaceFields.H"
 #include "syncTools.H"
 #include "CompactListList.H"
+#include "fvMeshTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -274,156 +274,6 @@ Foam::label Foam::fvMeshDistribute::findNonEmptyPatch() const
 }
 
 
-//// Appends processorPolyPatch. Returns patchID.
-//Foam::label Foam::fvMeshDistribute::addProcPatch
-//(
-//    const word& patchName,
-//    const label nbrProc
-//)
-//{
-//    // Clear local fields and e.g. polyMesh globalMeshData.
-//    mesh_.clearOut();
-//
-//
-//    polyBoundaryMesh& polyPatches =
-//        const_cast<polyBoundaryMesh&>(mesh_.boundaryMesh());
-//    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh_.boundary());
-//
-//    if (polyPatches.findPatchID(patchName) != -1)
-//    {
-//        FatalErrorIn
-//        (
-//            "fvMeshDistribute::addProcPatch(const word&, const label)"
-//        )
-//            << "Cannot create patch " << patchName << " since already exists."
-//            << nl
-//            << "Current patch names:" << polyPatches.names()
-//            << exit(FatalError);
-//    }
-//
-//
-//
-//    // Add the patch
-//    // ~~~~~~~~~~~~~
-//
-//    label sz = polyPatches.size();
-//
-//    // Add polyPatch
-//    polyPatches.setSize(sz+1);
-//    polyPatches.set
-//    (
-//        sz,
-//        new processorPolyPatch
-//        (
-//            patchName,
-//            0,              // size
-//            mesh_.nFaces(),
-//            sz,
-//            mesh_.boundaryMesh(),
-//            Pstream::myProcNo(),
-//            nbrProc
-//        )
-//    );
-//    fvPatches.setSize(sz+1);
-//    fvPatches.set
-//    (
-//        sz,
-//        fvPatch::New
-//        (
-//            polyPatches[sz],  // point to newly added polyPatch
-//            mesh_.boundary()
-//        )
-//    );
-//
-//    return sz;
-//}
-
-
-// Appends polyPatch. Returns patchID.
-Foam::label Foam::fvMeshDistribute::addPatch(polyPatch* patchPtr)
-{
-    // Clear local fields and e.g. polyMesh globalMeshData.
-    mesh_.clearOut();
-
-    polyBoundaryMesh& polyPatches =
-        const_cast<polyBoundaryMesh&>(mesh_.boundaryMesh());
-    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh_.boundary());
-
-    if (polyPatches.findPatchID(patchPtr->name()) != -1)
-    {
-        FatalErrorIn("fvMeshDistribute::addPatch(polyPatch*)")
-            << "Cannot create patch " << patchPtr->name()
-            << " since already exists." << nl
-            << "Current patch names:" << polyPatches.names()
-            << exit(FatalError);
-    }
-
-
-    // Add the patch
-    // ~~~~~~~~~~~~~
-
-    label sz = polyPatches.size();
-
-    // Add polyPatch
-    polyPatches.setSize(sz+1);
-    polyPatches.set(sz, patchPtr);
-    fvPatches.setSize(sz+1);
-    fvPatches.set
-    (
-        sz,
-        fvPatch::New
-        (
-            polyPatches[sz],  // point to newly added polyPatch
-            mesh_.boundary()
-        )
-    );
-
-    return sz;
-}
-
-
-// Deletes last patch
-void Foam::fvMeshDistribute::deleteTrailingPatch()
-{
-    // Clear local fields and e.g. polyMesh globalMeshData.
-    mesh_.clearOut();
-
-    polyBoundaryMesh& polyPatches =
-        const_cast<polyBoundaryMesh&>(mesh_.boundaryMesh());
-    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh_.boundary());
-
-    if (polyPatches.empty())
-    {
-        FatalErrorIn("fvMeshDistribute::deleteTrailingPatch(fvMesh&)")
-            << "No patches in mesh"
-            << abort(FatalError);
-    }
-
-    label sz = polyPatches.size();
-
-    label nFaces = polyPatches[sz-1].size();
-
-    // Remove last polyPatch
-    if (debug)
-    {
-        Pout<< "deleteTrailingPatch : Removing patch " << sz-1
-            << " : " << polyPatches[sz-1].name() << " size:" << nFaces << endl;
-    }
-
-    if (nFaces)
-    {
-        FatalErrorIn("fvMeshDistribute::deleteTrailingPatch()")
-            << "There are still " << nFaces << " faces in patch to be deleted "
-            << sz-1 << ' ' << polyPatches[sz-1].name()
-            << abort(FatalError);
-    }
-
-    // Remove actual patch
-    polyPatches.setSize(sz-1);
-    fvPatches.setSize(sz-1);
-}
-
-
 // Delete all processor patches. Move any processor faces into the last
 // non-processor patch.
 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::deleteProcPatches
@@ -469,25 +319,28 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::deleteProcPatches
 
 
     // Delete (now empty) processor patches.
-    forAllReverse(mesh_.boundaryMesh(), patchI)
     {
-        const polyPatch& pp = mesh_.boundaryMesh()[patchI];
-
-        if (isA<processorPolyPatch>(pp))
+        labelList oldToNew(identity(mesh_.boundaryMesh().size()));
+        label newI = 0;
+        // Non processor patches first
+        forAll(mesh_.boundaryMesh(), patchI)
         {
-            deleteTrailingPatch();
-            deleteTrailingPatchFields<volScalarField>();
-            deleteTrailingPatchFields<volVectorField>();
-            deleteTrailingPatchFields<volSphericalTensorField>();
-            deleteTrailingPatchFields<volSymmTensorField>();
-            deleteTrailingPatchFields<volTensorField>();
-
-            deleteTrailingPatchFields<surfaceScalarField>();
-            deleteTrailingPatchFields<surfaceVectorField>();
-            deleteTrailingPatchFields<surfaceSphericalTensorField>();
-            deleteTrailingPatchFields<surfaceSymmTensorField>();
-            deleteTrailingPatchFields<surfaceTensorField>();
+            if (!isA<processorPolyPatch>(mesh_.boundaryMesh()[patchI]))
+            {
+                oldToNew[patchI] = newI++;
+            }
         }
+        label nNonProcPatches = newI;
+
+        // Processor patches as last
+        forAll(mesh_.boundaryMesh(), patchI)
+        {
+            if (isA<processorPolyPatch>(mesh_.boundaryMesh()[patchI]))
+            {
+                oldToNew[patchI] = newI++;
+            }
+        }
+        fvMeshTools::reorderPatches(mesh_, oldToNew, nNonProcPatches, false);
     }
 
     return map;
@@ -1082,65 +935,28 @@ void Foam::fvMeshDistribute::addProcPatches
                       + "to"
                       + name(procI);
 
+                    processorPolyPatch pp
+                    (
+                        patchName,
+                        0,              // size
+                        mesh_.nFaces(),
+                        mesh_.boundaryMesh().size(),
+                        mesh_.boundaryMesh(),
+                        Pstream::myProcNo(),
+                        nbrProc[bFaceI]
+                    );
                     procPatchID[procI].insert
                     (
                         referPatchID[bFaceI],
-                        addPatch
+                        fvMeshTools::addPatch
                         (
-                            new processorPolyPatch
-                            (
-                                patchName,
-                                0,              // size
-                                mesh_.nFaces(),
-                                mesh_.boundaryMesh().size(),
-                                mesh_.boundaryMesh(),
-                                Pstream::myProcNo(),
-                                nbrProc[bFaceI]
-                            )
+                            mesh_,
+                            pp,
+                            dictionary(),
+                            processorFvPatchField<scalar>::typeName,
+                            false       // not parallel sync
                         )
                     );
-
-                    addPatchFields<volScalarField>
-                    (
-                        processorFvPatchField<scalar>::typeName
-                    );
-                    addPatchFields<volVectorField>
-                    (
-                        processorFvPatchField<vector>::typeName
-                    );
-                    addPatchFields<volSphericalTensorField>
-                    (
-                        processorFvPatchField<sphericalTensor>::typeName
-                    );
-                    addPatchFields<volSymmTensorField>
-                    (
-                        processorFvPatchField<symmTensor>::typeName
-                    );
-                    addPatchFields<volTensorField>
-                    (
-                        processorFvPatchField<tensor>::typeName
-                    );
-
-                    addPatchFields<surfaceScalarField>
-                    (
-                        processorFvPatchField<scalar>::typeName
-                    );
-                    addPatchFields<surfaceVectorField>
-                    (
-                        processorFvPatchField<vector>::typeName
-                    );
-                    addPatchFields<surfaceSphericalTensorField>
-                    (
-                        processorFvPatchField<sphericalTensor>::typeName
-                    );
-                    addPatchFields<surfaceSymmTensorField>
-                    (
-                        processorFvPatchField<symmTensor>::typeName
-                    );
-                    addPatchFields<surfaceTensorField>
-                    (
-                        processorFvPatchField<tensor>::typeName
-                    );
                 }
                 else
                 {
@@ -1158,66 +974,29 @@ void Foam::fvMeshDistribute::addProcPatches
                       + "through"
                       + cycName;
 
+                    processorCyclicPolyPatch pp
+                    (
+                        patchName,
+                        0,              // size
+                        mesh_.nFaces(),
+                        mesh_.boundaryMesh().size(),
+                        mesh_.boundaryMesh(),
+                        Pstream::myProcNo(),
+                        nbrProc[bFaceI],
+                        cycName
+                    );
                     procPatchID[procI].insert
                     (
                         referPatchID[bFaceI],
-                        addPatch
+                        fvMeshTools::addPatch
                         (
-                            new processorCyclicPolyPatch
-                            (
-                                patchName,
-                                0,              // size
-                                mesh_.nFaces(),
-                                mesh_.boundaryMesh().size(),
-                                mesh_.boundaryMesh(),
-                                Pstream::myProcNo(),
-                                nbrProc[bFaceI],
-                                cycName
-                            )
+                            mesh_,
+                            pp,
+                            dictionary(),   // optional per field patchField
+                            processorCyclicFvPatchField<scalar>::typeName,
+                            false           // not parallel sync
                         )
                     );
-
-                    addPatchFields<volScalarField>
-                    (
-                        processorCyclicFvPatchField<scalar>::typeName
-                    );
-                    addPatchFields<volVectorField>
-                    (
-                        processorCyclicFvPatchField<vector>::typeName
-                    );
-                    addPatchFields<volSphericalTensorField>
-                    (
-                        processorCyclicFvPatchField<sphericalTensor>::typeName
-                    );
-                    addPatchFields<volSymmTensorField>
-                    (
-                        processorCyclicFvPatchField<symmTensor>::typeName
-                    );
-                    addPatchFields<volTensorField>
-                    (
-                        processorCyclicFvPatchField<tensor>::typeName
-                    );
-
-                    addPatchFields<surfaceScalarField>
-                    (
-                        processorCyclicFvPatchField<scalar>::typeName
-                    );
-                    addPatchFields<surfaceVectorField>
-                    (
-                        processorCyclicFvPatchField<vector>::typeName
-                    );
-                    addPatchFields<surfaceSphericalTensorField>
-                    (
-                        processorCyclicFvPatchField<sphericalTensor>::typeName
-                    );
-                    addPatchFields<surfaceSymmTensorField>
-                    (
-                        processorCyclicFvPatchField<symmTensor>::typeName
-                    );
-                    addPatchFields<surfaceTensorField>
-                    (
-                        processorCyclicFvPatchField<tensor>::typeName
-                    );
                 }
             }
         }
@@ -2445,68 +2224,32 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
     // parallel comms. After this points and edges should again be consistent.
     mergeSharedPoints(constructPointMap);
 
-//    // Bit of hack: processorFvPatchField does not get reset since created
-//    // from nothing so explicitly reset.
-//    initPatchFields<volScalarField, processorFvPatchField<scalar> >
-//    (
-//        pTraits<scalar>::zero
-//    );
-//    initPatchFields<volVectorField, processorFvPatchField<vector> >
-//    (
-//        pTraits<vector>::zero
-//    );
-//    initPatchFields
-//    <
-//        volSphericalTensorField,
-//        processorFvPatchField<sphericalTensor>
-//    >
-//    (
-//        pTraits<sphericalTensor>::zero
-//    );
-//    initPatchFields<volSymmTensorField, processorFvPatchField<symmTensor> >
-//    (
-//        pTraits<symmTensor>::zero
-//    );
-//    initPatchFields<volTensorField, processorFvPatchField<tensor> >
-//    (
-//        pTraits<tensor>::zero
-//    );
-//    initPatchFields<surfaceScalarField, processorFvsPatchField<scalar> >
-//    (
-//        pTraits<scalar>::zero
-//    );
-//    initPatchFields<surfaceVectorField, processorFvsPatchField<vector> >
-//    (
-//        pTraits<vector>::zero
-//    );
-//    initPatchFields
-//    <
-//        surfaceSphericalTensorField,
-//        processorFvsPatchField<sphericalTensor>
-//    >
-//    (
-//        pTraits<sphericalTensor>::zero
-//    );
-//    initPatchFields
-//    <
-//        surfaceSymmTensorField,
-//        processorFvsPatchField<symmTensor>
-//    >
-//    (
-//        pTraits<symmTensor>::zero
-//    );
-//    initPatchFields<surfaceTensorField, processorFvsPatchField<tensor> >
-//    (
-//        pTraits<tensor>::zero
-//    );
-//XXXXX
     // Bit of hack: processorFvPatchField does not get reset since created
     // from nothing so explicitly reset.
-    correctBoundaryConditions<volScalarField>();
-    correctBoundaryConditions<volVectorField>();
-    correctBoundaryConditions<volSphericalTensorField>();
-    correctBoundaryConditions<volSymmTensorField>();
-    correctBoundaryConditions<volTensorField>();
+    initPatchFields<volScalarField, processorFvPatchField<scalar> >
+    (
+        pTraits<scalar>::zero
+    );
+    initPatchFields<volVectorField, processorFvPatchField<vector> >
+    (
+        pTraits<vector>::zero
+    );
+    initPatchFields
+    <
+        volSphericalTensorField,
+        processorFvPatchField<sphericalTensor>
+    >
+    (
+        pTraits<sphericalTensor>::zero
+    );
+    initPatchFields<volSymmTensorField, processorFvPatchField<symmTensor> >
+    (
+        pTraits<symmTensor>::zero
+    );
+    initPatchFields<volTensorField, processorFvPatchField<tensor> >
+    (
+        pTraits<tensor>::zero
+    );
 
     initPatchFields<surfaceScalarField, processorFvsPatchField<scalar> >
     (
@@ -2536,7 +2279,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
     (
         pTraits<tensor>::zero
     );
-//XXXXX
+
 
     mesh_.setInstance(mesh_.time().timeName());
 
diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
index eb771003a0..0aa64bf447 100644
--- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
+++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
@@ -102,20 +102,6 @@ class fvMeshDistribute
             //- Find patch to put exposed faces into.
             label findNonEmptyPatch() const;
 
-            //- Appends polyPatch. Returns patchID.
-            label addPatch(polyPatch*);
-
-            //- Add patch field
-            template<class GeoField>
-            void addPatchFields(const word& patchFieldType);
-
-            //- Deletes last patch.
-            void deleteTrailingPatch();
-
-            // Delete trailing patch fields
-            template<class GeoField>
-            void deleteTrailingPatchFields();
-
             //- Save boundary fields
             template <class T, class Mesh>
             void saveBoundaryFields
diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistributeTemplates.C b/src/dynamicMesh/fvMeshDistribute/fvMeshDistributeTemplates.C
index 00364a81ad..a09d538083 100644
--- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistributeTemplates.C
+++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistributeTemplates.C
@@ -55,65 +55,6 @@ void Foam::fvMeshDistribute::printFieldInfo(const fvMesh& mesh)
 }
 
 
-template<class GeoField>
-void Foam::fvMeshDistribute::addPatchFields(const word& patchFieldType)
-{
-    HashTable<const GeoField*> flds
-    (
-        mesh_.objectRegistry::lookupClass<GeoField>()
-    );
-
-    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
-    {
-        const GeoField& fld = *iter();
-
-        typename GeoField::GeometricBoundaryField& bfld =
-            const_cast<typename GeoField::GeometricBoundaryField&>
-            (
-                fld.boundaryField()
-            );
-
-        label sz = bfld.size();
-        bfld.setSize(sz + 1);
-        bfld.set
-        (
-            sz,
-            GeoField::PatchFieldType::New
-            (
-                patchFieldType,
-                mesh_.boundary()[sz],
-                fld.dimensionedInternalField()
-            )
-        );
-    }
-}
-
-
-// Delete trailing patch fields
-template<class GeoField>
-void Foam::fvMeshDistribute::deleteTrailingPatchFields()
-{
-    HashTable<const GeoField*> flds
-    (
-        mesh_.objectRegistry::lookupClass<GeoField>()
-    );
-
-    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
-    {
-        const GeoField& fld = *iter();
-
-        typename GeoField::GeometricBoundaryField& bfld =
-            const_cast<typename GeoField::GeometricBoundaryField&>
-            (
-                fld.boundaryField()
-            );
-
-        // Shrink patchFields
-        bfld.setSize(bfld.size() - 1);
-    }
-}
-
-
 // Save whole boundary field
 template <class T, class Mesh>
 void Foam::fvMeshDistribute::saveBoundaryFields
diff --git a/src/dynamicMesh/fvMeshTools/fvMeshTools.C b/src/dynamicMesh/fvMeshTools/fvMeshTools.C
new file mode 100644
index 0000000000..08d3cdc014
--- /dev/null
+++ b/src/dynamicMesh/fvMeshTools/fvMeshTools.C
@@ -0,0 +1,351 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvMeshTools.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Adds patch if not yet there. Returns patchID.
+Foam::label Foam::fvMeshTools::addPatch
+(
+    fvMesh& mesh,
+    const polyPatch& patch,
+    const dictionary& patchFieldDict,
+    const word& defaultPatchFieldType,
+    const bool validBoundary
+)
+{
+    polyBoundaryMesh& polyPatches =
+        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
+
+    label patchI = polyPatches.findPatchID(patch.name());
+    if (patchI != -1)
+    {
+        // Already there
+        return patchI;
+    }
+
+
+    label insertPatchI = polyPatches.size();
+    label startFaceI = mesh.nFaces();
+
+    forAll(polyPatches, patchI)
+    {
+        const polyPatch& pp = polyPatches[patchI];
+
+        if (isA<processorPolyPatch>(pp))
+        {
+            insertPatchI = patchI;
+            startFaceI = pp.start();
+            break;
+        }
+    }
+
+
+    // Below is all quite a hack. Feel free to change once there is a better
+    // mechanism to insert and reorder patches.
+
+    // Clear local fields and e.g. polyMesh parallelInfo.
+    mesh.clearOut();
+
+    label sz = polyPatches.size();
+
+    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
+
+    // Add polyPatch at the end
+    polyPatches.setSize(sz+1);
+    polyPatches.set
+    (
+        sz,
+        patch.clone
+        (
+            polyPatches,
+            insertPatchI,   //index
+            0,              //size
+            startFaceI      //start
+        )
+    );
+    fvPatches.setSize(sz+1);
+    fvPatches.set
+    (
+        sz,
+        fvPatch::New
+        (
+            polyPatches[sz],  // point to newly added polyPatch
+            mesh.boundary()
+        )
+    );
+
+    addPatchFields<volScalarField>
+    (
+        mesh,
+        patchFieldDict,
+        defaultPatchFieldType,
+        pTraits<scalar>::zero
+    );
+    addPatchFields<volVectorField>
+    (
+        mesh,
+        patchFieldDict,
+        defaultPatchFieldType,
+        pTraits<vector>::zero
+    );
+    addPatchFields<volSphericalTensorField>
+    (
+        mesh,
+        patchFieldDict,
+        defaultPatchFieldType,
+        pTraits<sphericalTensor>::zero
+    );
+    addPatchFields<volSymmTensorField>
+    (
+        mesh,
+        patchFieldDict,
+        defaultPatchFieldType,
+        pTraits<symmTensor>::zero
+    );
+    addPatchFields<volTensorField>
+    (
+        mesh,
+        patchFieldDict,
+        defaultPatchFieldType,
+        pTraits<tensor>::zero
+    );
+
+    // Surface fields
+
+    addPatchFields<surfaceScalarField>
+    (
+        mesh,
+        patchFieldDict,
+        defaultPatchFieldType,
+        pTraits<scalar>::zero
+    );
+    addPatchFields<surfaceVectorField>
+    (
+        mesh,
+        patchFieldDict,
+        defaultPatchFieldType,
+        pTraits<vector>::zero
+    );
+    addPatchFields<surfaceSphericalTensorField>
+    (
+        mesh,
+        patchFieldDict,
+        defaultPatchFieldType,
+        pTraits<sphericalTensor>::zero
+    );
+    addPatchFields<surfaceSymmTensorField>
+    (
+        mesh,
+        patchFieldDict,
+        defaultPatchFieldType,
+        pTraits<symmTensor>::zero
+    );
+    addPatchFields<surfaceTensorField>
+    (
+        mesh,
+        patchFieldDict,
+        defaultPatchFieldType,
+        pTraits<tensor>::zero
+    );
+
+    // Create reordering list
+    // patches before insert position stay as is
+    labelList oldToNew(sz+1);
+    for (label i = 0; i < insertPatchI; i++)
+    {
+        oldToNew[i] = i;
+    }
+    // patches after insert position move one up
+    for (label i = insertPatchI; i < sz; i++)
+    {
+        oldToNew[i] = i+1;
+    }
+    // appended patch gets moved to insert position
+    oldToNew[sz] = insertPatchI;
+
+    // Shuffle into place
+    polyPatches.reorder(oldToNew, validBoundary);
+    fvPatches.reorder(oldToNew);
+
+    reorderPatchFields<volScalarField>(mesh, oldToNew);
+    reorderPatchFields<volVectorField>(mesh, oldToNew);
+    reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
+    reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
+    reorderPatchFields<volTensorField>(mesh, oldToNew);
+    reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
+    reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
+    reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
+    reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
+    reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
+
+    return insertPatchI;
+}
+
+
+void Foam::fvMeshTools::setPatchFields
+(
+    fvMesh& mesh,
+    const label patchI,
+    const dictionary& patchFieldDict
+)
+{
+    setPatchFields<volScalarField>(mesh, patchI, patchFieldDict);
+    setPatchFields<volVectorField>(mesh, patchI, patchFieldDict);
+    setPatchFields<volSphericalTensorField>(mesh, patchI, patchFieldDict);
+    setPatchFields<volSymmTensorField>(mesh, patchI, patchFieldDict);
+    setPatchFields<volTensorField>(mesh, patchI, patchFieldDict);
+    setPatchFields<surfaceScalarField>(mesh, patchI, patchFieldDict);
+    setPatchFields<surfaceVectorField>(mesh, patchI, patchFieldDict);
+    setPatchFields<surfaceSphericalTensorField>
+    (
+        mesh,
+        patchI,
+        patchFieldDict
+    );
+    setPatchFields<surfaceSymmTensorField>(mesh, patchI, patchFieldDict);
+    setPatchFields<surfaceTensorField>(mesh, patchI, patchFieldDict);
+}
+
+
+void Foam::fvMeshTools::zeroPatchFields(fvMesh& mesh, const label patchI)
+{
+    setPatchFields<volScalarField>(mesh, patchI, pTraits<scalar>::zero);
+    setPatchFields<volVectorField>(mesh, patchI, pTraits<vector>::zero);
+    setPatchFields<volSphericalTensorField>
+    (
+        mesh,
+        patchI,
+        pTraits<sphericalTensor>::zero
+    );
+    setPatchFields<volSymmTensorField>
+    (
+        mesh,
+        patchI,
+        pTraits<symmTensor>::zero
+    );
+    setPatchFields<volTensorField>(mesh, patchI, pTraits<tensor>::zero);
+    setPatchFields<surfaceScalarField>(mesh, patchI, pTraits<scalar>::zero);
+    setPatchFields<surfaceVectorField>(mesh, patchI, pTraits<vector>::zero);
+    setPatchFields<surfaceSphericalTensorField>
+    (
+        mesh,
+        patchI,
+        pTraits<sphericalTensor>::zero
+    );
+    setPatchFields<surfaceSymmTensorField>
+    (
+        mesh,
+        patchI,
+        pTraits<symmTensor>::zero
+    );
+    setPatchFields<surfaceTensorField>(mesh, patchI, pTraits<tensor>::zero);
+}
+
+
+// Deletes last patch
+void Foam::fvMeshTools::trimPatches(fvMesh& mesh, const label nPatches)
+{
+    // Clear local fields and e.g. polyMesh globalMeshData.
+    mesh.clearOut();
+
+    polyBoundaryMesh& polyPatches =
+        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
+    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
+
+    if (polyPatches.empty())
+    {
+        FatalErrorIn("fvMeshTools::trimPatches(fvMesh&, const label)")
+            << "No patches in mesh"
+            << abort(FatalError);
+    }
+
+    label nFaces = 0;
+    for (label patchI = nPatches; patchI < polyPatches.size(); patchI++)
+    {
+        nFaces += polyPatches[patchI].size();
+    }
+    reduce(nFaces, sumOp<label>());
+
+    if (nFaces)
+    {
+        FatalErrorIn("fvMeshTools::trimPatches(fvMesh&, const label)")
+            << "There are still " << nFaces
+            << " faces in " << polyPatches.size()-nPatches
+            << " patches to be deleted" << abort(FatalError);
+    }
+
+    // Remove actual patches
+    polyPatches.setSize(nPatches);
+    fvPatches.setSize(nPatches);
+
+    trimPatchFields<volScalarField>(mesh, nPatches);
+    trimPatchFields<volVectorField>(mesh, nPatches);
+    trimPatchFields<volSphericalTensorField>(mesh, nPatches);
+    trimPatchFields<volSymmTensorField>(mesh, nPatches);
+    trimPatchFields<volTensorField>(mesh, nPatches);
+
+    trimPatchFields<surfaceScalarField>(mesh, nPatches);
+    trimPatchFields<surfaceVectorField>(mesh, nPatches);
+    trimPatchFields<surfaceSphericalTensorField>(mesh, nPatches);
+    trimPatchFields<surfaceSymmTensorField>(mesh, nPatches);
+    trimPatchFields<surfaceTensorField>(mesh, nPatches);
+}
+
+
+void Foam::fvMeshTools::reorderPatches
+(
+    fvMesh& mesh,
+    const labelList& oldToNew,
+    const label nNewPatches,
+    const bool validBoundary
+)
+{
+    polyBoundaryMesh& polyPatches =
+        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
+    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
+
+    // Shuffle into place
+    polyPatches.reorder(oldToNew, validBoundary);
+    fvPatches.reorder(oldToNew);
+
+    reorderPatchFields<volScalarField>(mesh, oldToNew);
+    reorderPatchFields<volVectorField>(mesh, oldToNew);
+    reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
+    reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
+    reorderPatchFields<volTensorField>(mesh, oldToNew);
+    reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
+    reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
+    reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
+    reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
+    reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
+
+    // Remove last.
+    trimPatches(mesh, nNewPatches);
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/fvMeshTools/fvMeshTools.H b/src/dynamicMesh/fvMeshTools/fvMeshTools.H
new file mode 100644
index 0000000000..a5d443fffd
--- /dev/null
+++ b/src/dynamicMesh/fvMeshTools/fvMeshTools.H
@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fvMeshTools
+
+Description
+    A collection of tools for operating on an fvMesh.
+
+SourceFiles
+    fvMeshTools.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fvMeshTools_H
+#define fvMeshTools_H
+
+#include "fvMesh.H"
+//#include "HashSet.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class fvMeshTools Declaration
+\*---------------------------------------------------------------------------*/
+
+class fvMeshTools
+{
+    template<class GeoField>
+    static void addPatchFields
+    (
+        fvMesh&,
+        const dictionary& patchFieldDict,
+        const word& defaultPatchFieldType,
+        const typename GeoField::value_type& defaultPatchValue
+    );
+
+    //- Set patchFields according to dictionary
+    template<class GeoField>
+    static void setPatchFields
+    (
+        fvMesh& mesh,
+        const label patchI,
+        const dictionary& patchFieldDict
+    );
+
+    //- Set patchFields to value
+    template<class GeoField>
+    static void setPatchFields
+    (
+        fvMesh& mesh,
+        const label patchI,
+        const typename GeoField::value_type& value
+    );
+
+    // Remove last patch fields
+    template<class GeoField>
+    static void trimPatchFields(fvMesh&, const label nPatches);
+
+    template<class GeoField>
+    static void reorderPatchFields(fvMesh&, const labelList& oldToNew);
+
+    // Remove trialing patches
+    static void trimPatches(fvMesh&, const label nPatches);
+
+
+public:
+
+    //- Add patch. Supply per field the new patchField per field as a
+    //  subdictionary or a default type. If validBoundary call is parallel
+    //  synced and all add the same patch with same settings
+    static label addPatch
+    (
+        fvMesh& mesh,
+        const polyPatch& patch,
+        const dictionary& patchFieldDict,
+        const word& defaultPatchFieldType,
+        const bool validBoundary
+    );
+
+    //- Change patchField on registered fields according to dictionary
+    static void setPatchFields
+    (
+        fvMesh& mesh,
+        const label patchI,
+        const dictionary& patchFieldDict
+    );
+
+    //- Change patchField to zero on registered fields
+    static void zeroPatchFields(fvMesh& mesh, const label patchI);
+
+    // -Reorder and remove trailing patches. If validBoundary call is parallel
+    //  synced and all add the same patch with same settings
+    static void reorderPatches
+    (
+        fvMesh&,
+        const labelList& oldToNew,
+        const label nPatches,
+        const bool validBoundary
+    );
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "fvMeshToolsTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/fvMeshTools/fvMeshToolsTemplates.C b/src/dynamicMesh/fvMeshTools/fvMeshToolsTemplates.C
new file mode 100644
index 0000000000..3e9ca9d088
--- /dev/null
+++ b/src/dynamicMesh/fvMeshTools/fvMeshToolsTemplates.C
@@ -0,0 +1,208 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvMeshTools.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class GeoField>
+void Foam::fvMeshTools::addPatchFields
+(
+    fvMesh& mesh,
+    const dictionary& patchFieldDict,
+    const word& defaultPatchFieldType,
+    const typename GeoField::value_type& defaultPatchValue
+)
+{
+    HashTable<const GeoField*> flds
+    (
+        mesh.objectRegistry::lookupClass<GeoField>()
+    );
+
+    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
+    {
+        const GeoField& fld = *iter();
+
+        typename GeoField::GeometricBoundaryField& bfld =
+            const_cast<typename GeoField::GeometricBoundaryField&>
+            (
+                fld.boundaryField()
+            );
+
+        label sz = bfld.size();
+        bfld.setSize(sz+1);
+
+        if (patchFieldDict.found(fld.name()))
+        {
+            bfld.set
+            (
+                sz,
+                GeoField::PatchFieldType::New
+                (
+                    mesh.boundary()[sz],
+                    fld.dimensionedInternalField(),
+                    patchFieldDict.subDict(fld.name())
+                )
+            );
+        }
+        else
+        {
+            bfld.set
+            (
+                sz,
+                GeoField::PatchFieldType::New
+                (
+                    defaultPatchFieldType,
+                    mesh.boundary()[sz],
+                    fld.dimensionedInternalField()
+                )
+            );
+            bfld[sz] == defaultPatchValue;
+        }
+    }
+}
+
+
+template<class GeoField>
+void Foam::fvMeshTools::setPatchFields
+(
+    fvMesh& mesh,
+    const label patchI,
+    const dictionary& patchFieldDict
+)
+{
+    HashTable<const GeoField*> flds
+    (
+        mesh.objectRegistry::lookupClass<GeoField>()
+    );
+
+    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
+    {
+        const GeoField& fld = *iter();
+
+        typename GeoField::GeometricBoundaryField& bfld =
+            const_cast<typename GeoField::GeometricBoundaryField&>
+            (
+                fld.boundaryField()
+            );
+
+        if (patchFieldDict.found(fld.name()))
+        {
+            bfld.set
+            (
+                patchI,
+                GeoField::PatchFieldType::New
+                (
+                    mesh.boundary()[patchI],
+                    fld.dimensionedInternalField(),
+                    patchFieldDict.subDict(fld.name())
+                )
+            );
+        }
+    }
+}
+
+
+
+
+template<class GeoField>
+void Foam::fvMeshTools::setPatchFields
+(
+    fvMesh& mesh,
+    const label patchI,
+    const typename GeoField::value_type& value
+)
+{
+    HashTable<const GeoField*> flds
+    (
+        mesh.objectRegistry::lookupClass<GeoField>()
+    );
+
+    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
+    {
+        const GeoField& fld = *iter();
+
+        typename GeoField::GeometricBoundaryField& bfld =
+            const_cast<typename GeoField::GeometricBoundaryField&>
+            (
+                fld.boundaryField()
+            );
+
+        bfld[patchI] == value;
+    }
+}
+
+
+// Remove last patch field
+template<class GeoField>
+void Foam::fvMeshTools::trimPatchFields(fvMesh& mesh, const label nPatches)
+{
+    HashTable<const GeoField*> flds
+    (
+        mesh.objectRegistry::lookupClass<GeoField>()
+    );
+
+    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
+    {
+        const GeoField& fld = *iter();
+
+        const_cast<typename GeoField::GeometricBoundaryField&>
+        (
+            fld.boundaryField()
+        ).setSize(nPatches);
+    }
+}
+
+
+// Reorder patch field
+template<class GeoField>
+void Foam::fvMeshTools::reorderPatchFields
+(
+    fvMesh& mesh,
+    const labelList& oldToNew
+)
+{
+    HashTable<const GeoField*> flds
+    (
+        mesh.objectRegistry::lookupClass<GeoField>()
+    );
+
+    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
+    {
+        const GeoField& fld = *iter();
+
+        typename GeoField::GeometricBoundaryField& bfld =
+            const_cast<typename GeoField::GeometricBoundaryField&>
+            (
+                fld.boundaryField()
+            );
+        bfld.reorder(oldToNew);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H
index bc3624548e..9e8269df93 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,8 +66,8 @@ Description
 
 
     - added faces get same patchID as face they are extruded from
-    - 'side' faces (i.e. on the edge of pp) get the patchID of the
-    other patch they are connected to.
+    - 'side' faces (i.e. on the edge of pp) get the patchID/zoneID of the
+    other patch/zone they are connected to (hopefully only 1)
 
 
     E.g. 3 boundary faces on patches a,b. b gets extruded, a doesn't.
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index 2621adbc23..b35aff1ab4 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -128,6 +128,7 @@ faceZoneSources = sets/faceZoneSources
 $(faceZoneSources)/faceZoneToFaceZone/faceZoneToFaceZone.C
 $(faceZoneSources)/setsToFaceZone/setsToFaceZone.C
 $(faceZoneSources)/setToFaceZone/setToFaceZone.C
+$(faceZoneSources)/searchableSurfaceToFaceZone/searchableSurfaceToFaceZone.C
 
 cellZoneSources = sets/cellZoneSources
 $(cellZoneSources)/setToCellZone/setToCellZone.C
diff --git a/src/meshTools/sets/faceZoneSources/searchableSurfaceToFaceZone/searchableSurfaceToFaceZone.C b/src/meshTools/sets/faceZoneSources/searchableSurfaceToFaceZone/searchableSurfaceToFaceZone.C
new file mode 100644
index 0000000000..b91ec33b82
--- /dev/null
+++ b/src/meshTools/sets/faceZoneSources/searchableSurfaceToFaceZone/searchableSurfaceToFaceZone.C
@@ -0,0 +1,211 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "searchableSurfaceToFaceZone.H"
+#include "polyMesh.H"
+#include "faceZoneSet.H"
+#include "searchableSurface.H"
+#include "syncTools.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(searchableSurfaceToFaceZone, 0);
+    addToRunTimeSelectionTable
+    (
+        topoSetSource,
+        searchableSurfaceToFaceZone,
+        word
+    );
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::searchableSurfaceToFaceZone::usage_
+(
+    searchableSurfaceToFaceZone::typeName,
+    "\n    Usage: searchableSurfaceToFaceZone surface\n\n"
+    "    Select all faces whose cell-cell centre vector intersects the surface "
+    "\n"
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from dictionary
+Foam::searchableSurfaceToFaceZone::searchableSurfaceToFaceZone
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    topoSetSource(mesh),
+    surfacePtr_
+    (
+        searchableSurface::New
+        (
+            word(dict.lookup("surface")),
+            mesh.objectRegistry::db(),
+            dict
+        )
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceToFaceZone::~searchableSurfaceToFaceZone()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::searchableSurfaceToFaceZone::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if (!isA<faceZoneSet>(set))
+    {
+        WarningIn
+        (
+            "searchableSurfaceToFaceZone::applyToSet"
+            "(const topoSetSource::setAction"
+            ", topoSet"
+        )   << "Operation only allowed on a faceZoneSet." << endl;
+    }
+    else
+    {
+        faceZoneSet& fzSet = refCast<faceZoneSet>(set);
+
+        // Get cell-cell centre vectors
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        pointField start(mesh_.nFaces());
+        pointField end(mesh_.nFaces());
+
+        const pointField& cc = mesh_.cellCentres();
+
+        // Internal faces
+        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+        {
+            start[faceI] = cc[mesh_.faceOwner()[faceI]];
+            end[faceI] = cc[mesh_.faceNeighbour()[faceI]];
+        }
+
+        // Boundary faces
+        vectorField nbrCellCentres;
+        syncTools::swapBoundaryCellList(mesh_, cc, nbrCellCentres);
+
+        const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
+        forAll(pbm, patchI)
+        {
+            const polyPatch& pp = pbm[patchI];
+
+            if (pp.coupled())
+            {
+                forAll(pp, i)
+                {
+                    label faceI = pp.start()+i;
+                    start[faceI] = cc[mesh_.faceOwner()[faceI]];
+                    end[faceI] = nbrCellCentres[faceI-mesh_.nInternalFaces()];
+                }
+            }
+            else
+            {
+                forAll(pp, i)
+                {
+                    label faceI = pp.start()+i;
+                    start[faceI] = cc[mesh_.faceOwner()[faceI]];
+                    end[faceI] = mesh_.faceCentres()[faceI];
+                }
+            }
+        }
+
+
+        // Do all intersection tests
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        List<pointIndexHit> hits;
+        surfacePtr_().findLine(start, end, hits);
+        pointField normals;
+        surfacePtr_().getNormal(hits, normals);
+
+
+        // Select intersected faces
+        // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+        if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+        {
+            Info<< "    Adding all faces from surface "
+                << surfacePtr_().name() << " ..." << endl;
+
+            DynamicList<label> newAddressing(fzSet.addressing());
+            DynamicList<bool> newFlipMap(fzSet.flipMap());
+
+            forAll(hits, faceI)
+            {
+                if (hits[faceI].hit() && !fzSet.found(faceI))
+                {
+                    newAddressing.append(faceI);
+                    vector d = end[faceI]-start[faceI];
+                    newFlipMap.append((normals[faceI] & d) < 0);
+                }
+            }
+
+            fzSet.addressing().transfer(newAddressing);
+            fzSet.flipMap().transfer(newFlipMap);
+            fzSet.updateSet();
+        }
+        else if (action == topoSetSource::DELETE)
+        {
+            Info<< "    Removing all faces from surface "
+                << surfacePtr_().name() << " ..." << endl;
+
+            // Start off empty
+            DynamicList<label> newAddressing(fzSet.addressing().size());
+            DynamicList<bool> newFlipMap(fzSet.flipMap().size());
+
+            forAll(fzSet.addressing(), i)
+            {
+                if (!hits[fzSet.addressing()[i]].hit())
+                {
+                    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/searchableSurfaceToFaceZone/searchableSurfaceToFaceZone.H b/src/meshTools/sets/faceZoneSources/searchableSurfaceToFaceZone/searchableSurfaceToFaceZone.H
new file mode 100644
index 0000000000..7269a3d38e
--- /dev/null
+++ b/src/meshTools/sets/faceZoneSources/searchableSurfaceToFaceZone/searchableSurfaceToFaceZone.H
@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::searchableSurfaceToFaceZone
+
+Description
+    A topoSetSource to select faces based on intersection (of cell-cell
+    vector) with a surface.
+
+SourceFiles
+    searchableSurfaceToFaceZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchableSurfaceToFaceZone_H
+#define searchableSurfaceToFaceZone_H
+
+#include "topoSetSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class searchableSurface;
+
+/*---------------------------------------------------------------------------*\
+                   Class searchableSurfaceToFaceZone Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchableSurfaceToFaceZone
+:
+    public topoSetSource
+{
+    // Private data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- Surface
+        autoPtr<searchableSurface> surfacePtr_;
+
+public:
+
+    //- Runtime type information
+    TypeName("searchableSurfaceToFaceZone");
+
+    // Constructors
+
+        //- Construct from dictionary
+        searchableSurfaceToFaceZone
+        (
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~searchableSurfaceToFaceZone();
+
+
+    // Member Functions
+
+        virtual sourceType setType() const
+        {
+            return FACEZONESOURCE;
+        }
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet&
+        ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/controlDict b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/controlDict
index 5d0174164b..97f299f8a5 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/controlDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/controlDict
@@ -15,6 +15,12 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+libs
+(
+    "libcompressibleTurbulenceModel.so"
+    "libcompressibleRASModels.so"
+);
+
 application     chtMultiRegionFoam;
 
 startFrom       latestTime;