diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
index b1bcf33d82..be4fffa8df 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
@@ -132,456 +132,12 @@ Notes:
 #include "surfaceFields.H"
 #include "pointFields.H"
 //#include "ReadFields.H"
-
+#include "fvMeshTools.H"
 
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-//template<class GeoField>
-//void addPatchFields(const 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();
-//
-//        for (label i = 0; i < sz; i++)
-//        {
-//            bfld.set
-//            (
-//                i,
-//                bfld.clone(GeoField::PatchFieldType::New
-//                (
-//                    patchFieldType,
-//                    fld.mesh().boundary()[sz],
-//                    fld.dimensionedInternalField()
-//                )
-//            );
-//
-//
-//
-//        Pout<< "fld:" << fld.name() << " had " << sz << " patches." << endl;
-//        Pout<< "fld before:" << fld << endl;
-//        Pout<< "adding on patch:" << fld.mesh().boundary()[sz].name() << endl;
-//
-//        bfld.setSize(sz+1);
-//        bfld.set
-//        (
-//            sz,
-//            GeoField::PatchFieldType::New
-//            (
-//                patchFieldType,
-//                fld.mesh().boundary()[sz],
-//                fld.dimensionedInternalField()
-//            )
-//        );
-//
-//        bfld[sz].operator=(pTraits<typename GeoField::value_type>::zero);
-//
-//        Pout<< "fld:" << fld.name() << " now " << bfld.size() << " patches."
-//            << endl;
-//
-//        const typename GeoField::PatchFieldType& pfld = bfld[sz];
-//        Pout<< "pfld:" << pfld << endl;
-//
-//
-//        Pout<< "fld value:" << fld << endl;
-//    }
-//}
-
-
-// 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);
-    }
-}
-
-
-//void addCalculatedPatchFields(const fvMesh& mesh)
-//{
-//    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
-//    );
-//
-//    // Point fields
-//
-//    addPatchFields<pointScalarField>
-//    (
-//        mesh,
-//        calculatedFvPatchField<scalar>::typeName
-//    );
-//    addPatchFields<pointVectorField>
-//    (
-//        mesh,
-//        calculatedFvPatchField<vector>::typeName
-//    );
-//}
-//
-//
-//void addAllPatchFields(fvMesh& mesh, const label insertPatchI)
-//{
-//    polyBoundaryMesh& polyPatches =
-//        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
-//    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
-//
-//    label sz = polyPatches.size();
-//
-//    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);
-//    for (label i = 0; i < insertPatchI; i++)
-//    {
-//        oldToNew[i] = i;
-//    }
-//    // patches after insert position move one up
-//    for (label i = insertPatchI; i < sz-1; i++)
-//    {
-//        oldToNew[i] = i+1;
-//    }
-//    // appended patch gets moved to insert position
-//    oldToNew[sz-1] = 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);
-//}
-
-
-//// Adds patch if not yet there. Returns patchID.
-//template<class PatchType>
-//label addPatch(fvMesh& mesh, const word& patchName, const dictionary& dict)
-//{
-//    polyBoundaryMesh& polyPatches =
-//        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
-//
-//    label patchI = polyPatches.findPatchID(patchName);
-//    if (patchI != -1)
-//    {
-//        if (isA<PatchType>(polyPatches[patchI]))
-//        {
-//            // Already there
-//            return patchI;
-//        }
-//        else
-//        {
-//            FatalErrorIn("addPatch<PatchType>(fvMesh&, const word&)")
-//                << "Already have patch " << patchName
-//                << " but of type " << PatchType::typeName
-//                << 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;
-//        }
-//    }
-//
-//    dictionary patchDict(dict);
-//    patchDict.set("type", PatchType::typeName);
-//    patchDict.set("nFaces", 0);
-//    patchDict.set("startFace", startFaceI);
-//
-//
-//    // 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,
-//        polyPatch::New
-//        (
-//            patchName,
-//            patchDict,
-//            insertPatchI,
-//            polyPatches
-//        )
-//    );
-//    fvPatches.setSize(sz+1);
-//    fvPatches.set
-//    (
-//        sz,
-//        fvPatch::New
-//        (
-//            polyPatches[sz],  // point to newly added polyPatch
-//            mesh.boundary()
-//        )
-//    );
-//
-//    addAllPatchFields(mesh, insertPatchI);
-//
-//    return insertPatchI;
-//}
-//
-//
-//template<class PatchType>
-//label addPatch(fvMesh& mesh, const word& patchName)
-//{
-//Pout<< "addPatch:" << patchName << endl;
-//
-//    polyBoundaryMesh& polyPatches =
-//        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
-//
-//    label patchI = polyPatches.findPatchID(patchName);
-//    if (patchI != -1)
-//    {
-//        if (isA<PatchType>(polyPatches[patchI]))
-//        {
-//            // Already there
-//            return patchI;
-//        }
-//        else
-//        {
-//            FatalErrorIn("addPatch<PatchType>(fvMesh&, const word&)")
-//                << "Already have patch " << patchName
-//                << " but of type " << PatchType::typeName
-//                << 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,
-//        polyPatch::New
-//        (
-//            PatchType::typeName,
-//            patchName,
-//            0,              // size
-//            startFaceI,
-//            insertPatchI,
-//            polyPatches
-//        )
-//    );
-//    fvPatches.setSize(sz+1);
-//    fvPatches.set
-//    (
-//        sz,
-//        fvPatch::New
-//        (
-//            polyPatches[sz],  // point to newly added polyPatch
-//            mesh.boundary()
-//        )
-//    );
-//
-//    addAllPatchFields(mesh, insertPatchI);
-//
-//    return insertPatchI;
-//}
-
-
 label findPatchID(const List<polyPatch*>& newPatches, const word& name)
 {
     forAll(newPatches, i)
@@ -712,53 +268,6 @@ label addPatch
 }
 
 
-// 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);
-    reorderPatchFields<pointScalarField>(mesh, oldToNew);
-    reorderPatchFields<pointVectorField>(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);
-    trimPatchFields<pointScalarField>(mesh, nNewPatches);
-    trimPatchFields<pointVectorField>(mesh, nNewPatches);
-}
-
-
 // Remove zero-sized patches
 void deleteEmptyPatches(fvMesh& mesh)
 {
@@ -837,7 +346,7 @@ void deleteEmptyPatches(fvMesh& mesh)
         }
     }
 
-    reorderPatches(mesh, oldToNew, usedI);
+    fvMeshTools::reorderPatches(mesh, oldToNew, usedI, true);
 }
 
 
@@ -1805,18 +1314,26 @@ int main(int argc, char *argv[])
         mappedPatchBase::sampleModeNames_[dict.lookup("sampleMode")];
 
     const Switch oneD(dict.lookup("oneD"));
+    Switch oneDNonManifoldEdges(false);
+    word oneDPatchType(emptyPolyPatch::typeName);
+    if (oneD)
+    {
+        oneDNonManifoldEdges = dict.lookupOrDefault("nonManifold", false);
+        dict.lookup("oneDPolyPatchType") >> oneDPatchType;
+    }
+
     const Switch adaptMesh(dict.lookup("adaptMesh"));
 
     if (hasZones)
     {
-        Pout<< "Extruding zones " << zoneNames
+        Info<< "Extruding zones " << zoneNames
             << " on mesh " << regionName
             << " into shell mesh " << shellRegionName
             << endl;
     }
     else
     {
-        Pout<< "Extruding faceSets " << zoneNames
+        Info<< "Extruding faceSets " << zoneNames
             << " on mesh " << regionName
             << " into shell mesh " << shellRegionName
             << endl;
@@ -1832,6 +1349,26 @@ int main(int argc, char *argv[])
     }
 
 
+    if (oneD)
+    {
+        if (oneDNonManifoldEdges)
+        {
+            Info<< "Extruding as 1D columns with sides in patch type "
+                << oneDPatchType
+                << " and connected points (except on non-manifold areas)."
+                << endl;
+        }
+        else
+        {
+            Info<< "Extruding as 1D columns with sides in patch type "
+                << oneDPatchType
+                << " and duplicated points (overlapping volumes)."
+                << endl;
+        }
+    }
+
+
+
 
     //// Read objects in time directory
     //IOobjectList objects(mesh, runTime.timeName());
@@ -1894,7 +1431,7 @@ int main(int argc, char *argv[])
     {
         meshInstance = oldInstance;
     }
-    Pout<< "Writing meshes to " << meshInstance << nl << endl;
+    Info<< "Writing meshes to " << meshInstance << nl << endl;
 
 
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
@@ -2130,7 +1667,7 @@ int main(int argc, char *argv[])
     const labelListList& edgeFaces = extrudePatch.edgeFaces();
 
 
-    Pout<< "extrudePatch :"
+    Info<< "extrudePatch :"
         << " faces:" << extrudePatch.size()
         << " points:" << extrudePatch.nPoints()
         << " edges:" << extrudePatch.nEdges()
@@ -2325,7 +1862,7 @@ int main(int argc, char *argv[])
     (
         mesh,
         zoneNames,
-        (oneD ? dict.lookup("oneDPolyPatchType") : word::null),
+        (oneD ? oneDPatchType : word::null),
 
         regionPatches,
         zoneSidePatch
@@ -2416,10 +1953,18 @@ int main(int argc, char *argv[])
             {
                 ePatches[i] = zoneSidePatch[zoneID[eFaces[i]]];
             }
-            //- Set nonManifoldEdge[edgeI] for non-manifold edges only
-            //  The other option is to have non-manifold edges everywhere
-            //  and generate space overlapping columns of cells.
-            if (eFaces.size() != 2)
+
+            if (oneDNonManifoldEdges)
+            {
+                //- Set nonManifoldEdge[edgeI] for non-manifold edges only
+                //  The other option is to have non-manifold edges everywhere
+                //  and generate space overlapping columns of cells.
+                if (eFaces.size() != 2)
+                {
+                    nonManifoldEdge[edgeI] = 1;
+                }
+            }
+            else
             {
                 nonManifoldEdge[edgeI] = 1;
             }
@@ -2834,7 +2379,7 @@ int main(int argc, char *argv[])
         "point to patch point addressing";
 
 
-    Pout<< "Writing mesh " << regionMesh.name()
+    Info<< "Writing mesh " << regionMesh.name()
         << " to " << regionMesh.facesInstance() << nl
         << endl;
 
@@ -3013,7 +2558,7 @@ int main(int argc, char *argv[])
         // Remove any unused patches
         deleteEmptyPatches(mesh);
 
-        Pout<< "Writing mesh " << mesh.name()
+        Info<< "Writing mesh " << mesh.name()
             << " to " << mesh.facesInstance() << nl
             << endl;
 
diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMeshDict b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMeshDict
index a768615f85..ffa2db8b19 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMeshDict
+++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMeshDict
@@ -40,11 +40,20 @@ adaptMesh       true;
 // Sample mode for inter-region communication
 sampleMode      nearestPatchFace;
 
-// Extrude 1D-columns of cells?
-oneD            false;
 
-// If oneD is true. Specify which boundary is wanted between the layers
-//oneDPolyPatchType empty; //wedge
+// 1 D extrusion
+// ~~~~~~~~~~~~~
+
+    // Extrude 1D-columns of cells? This by default duplicates points so can
+    // have overlapping columns (i.e. non space filling)
+    oneD            false;
+
+    //- If oneD: specify which boundary is wanted between the layers
+    //oneDPolyPatchType empty; //wedge
+
+    //- If oneD: specify whether to duplicate points (i.e. disconnect 1D
+    //  columns) or only on non-manifold extrusion areas. Default is false.
+    //  nonManifold true;
 
 
 //- Extrusion model to use. The only logical choice is linearNormal?
diff --git a/applications/utilities/mesh/manipulation/createBaffles/Make/files b/applications/utilities/mesh/manipulation/createBaffles/Make/files
index 4d4c0ea538..734f2bce76 100644
--- a/applications/utilities/mesh/manipulation/createBaffles/Make/files
+++ b/applications/utilities/mesh/manipulation/createBaffles/Make/files
@@ -1,3 +1,7 @@
+faceSelection/faceSelection.C
+faceSelection/faceZoneSelection.C
+faceSelection/searchableSurfaceSelection.C
+
 createBaffles.C
 
 EXE = $(FOAM_APPBIN)/createBaffles
diff --git a/applications/utilities/mesh/manipulation/createBaffles/Make/options b/applications/utilities/mesh/manipulation/createBaffles/Make/options
index f7e0c60fed..3c8db28129 100644
--- a/applications/utilities/mesh/manipulation/createBaffles/Make/options
+++ b/applications/utilities/mesh/manipulation/createBaffles/Make/options
@@ -1,4 +1,5 @@
 EXE_INC = \
+    -IfaceSelection \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
diff --git a/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C b/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C
index 6fb57aa027..c2e2a06ebb 100644
--- a/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C
+++ b/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C
@@ -26,25 +26,26 @@ Description
     mergeOrSplitBaffles.
 
     Note: if any coupled patch face is selected for baffling the opposite
-    member has to be selected for baffling as well. Note that this
-    is the same as repatching. This was added only for convenience so
-    you don't have to filter coupled boundary out of your set.
+    member has to be selected for baffling as well.
+
+    - if the patch already exists will not override it nor its fields
+    - if the patch does not exist it will be created together with 'calculated'
+      patchfields unless the field is mentioned in the patchFields section.
 
 \*---------------------------------------------------------------------------*/
 
-#include "syncTools.H"
 #include "argList.H"
 #include "Time.H"
-#include "faceSet.H"
 #include "polyTopoChange.H"
 #include "polyModifyFace.H"
 #include "polyAddFace.H"
 #include "ReadFields.H"
 #include "volFields.H"
 #include "surfaceFields.H"
-#include "ZoneIDs.H"
 #include "fvMeshMapper.H"
-#include "SetPatchFields.H"
+#include "faceSelection.H"
+
+#include "fvMeshTools.H"
 
 using namespace Foam;
 
@@ -107,21 +108,6 @@ void modifyOrAddFace
 }
 
 
-label findPatchID(const polyMesh& mesh, const word& name)
-{
-    const label patchI = mesh.boundaryMesh().findPatchID(name);
-
-    if (patchI == -1)
-    {
-        FatalErrorIn("findPatchID(const polyMesh&, const word&)")
-            << "Cannot find patch " << name << endl
-            << "Valid patches are " << mesh.boundaryMesh().names()
-            << exit(FatalError);
-    }
-    return patchI;
-}
-
-
 // Main program:
 
 int main(int argc, char *argv[])
@@ -129,102 +115,84 @@ int main(int argc, char *argv[])
     argList::addNote
     (
         "Makes internal faces into boundary faces.\n"
-        "Does not duplicate points, unlike mergeOrSplitBaffles."
+        "Does not duplicate points."
     );
 
     #include "addOverwriteOption.H"
-    #include "addRegionOption.H"
-
-    argList::validArgs.append("faceZone");
-    argList::validArgs.append("(masterPatch slavePatch)");
     argList::addOption
     (
-        "additionalPatches",
-        "((master2 slave2) .. (masterN slaveN))"
+        "dict",
+        "file",
+        "specify alternative dictionary for the createBaffles description"
     );
-    argList::addBoolOption
-    (
-        "internalFacesOnly",
-        "do not convert boundary faces"
-    );
-    argList::addBoolOption
-    (
-        "updateFields",
-        "update fields to include new patches:"
-        " NOTE: updated field values may need to be edited"
-    );
-
+    #include "addRegionOption.H"
     #include "setRootCase.H"
     #include "createTime.H"
     runTime.functionObjects().off();
     #include "createNamedMesh.H"
 
+
+    const bool overwrite = args.optionFound("overwrite");
+
     const word oldInstance = mesh.pointsInstance();
 
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-    const faceZoneMesh& faceZones = mesh.faceZones();
-
-    // Faces to baffle
-    faceZoneID zoneID(args.additionalArgs()[0], faceZones);
-
-    Info<< "Converting faces on zone " << zoneID.name()
-        << " into baffles." << nl << endl;
-
-    if (zoneID.index() == -1)
-    {
-        FatalErrorIn(args.executable()) << "Cannot find faceZone "
-            << zoneID.name() << endl
-            << "Valid zones are " << faceZones.names()
-            << exit(FatalError);
-    }
-
-    const faceZone& fZone = faceZones[zoneID.index()];
-
-    Info<< "Found " << returnReduce(fZone.size(), sumOp<label>())
-        << " faces on zone " << zoneID.name() << nl << endl;
-
-    // Make sure patches and zoneFaces are synchronised across couples
-    patches.checkParallelSync(true);
-    fZone.checkParallelSync(true);
-
-    // Patches to put baffles into
-    DynamicList<label> newMasterPatches(1);
-    DynamicList<label> newSlavePatches(1);
-
-    const Pair<word> patchNames(IStringStream(args.additionalArgs()[1])());
-    newMasterPatches.append(findPatchID(mesh, patchNames[0]));
-    newSlavePatches.append(findPatchID(mesh, patchNames[1]));
-    Info<< "Using master patch " << patchNames[0]
-        << " at index " << newMasterPatches[0] << endl;
-    Info<< "Using slave patch " << patchNames[1]
-        << " at index " << newSlavePatches[0] << endl;
-
-
-    // Additional patches
-    if (args.optionFound("additionalPatches"))
-    {
-        const List<Pair<word> > patchNames
+    const word dictName
+    (
+        args.optionLookupOrDefault<word>
         (
-            args.optionLookup("additionalPatches")()
+            "dict",
+            "createBafflesDict"
+        )
+    );
+
+    Switch internalFacesOnly(false);
+
+    Switch noFields(false);
+
+    PtrList<faceSelection> selectors;
+    {
+        Info<< "Reading baffle criteria from " << dictName << nl << endl;
+        IOdictionary dict
+        (
+            IOobject
+            (
+                dictName,
+                mesh.time().system(),
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            )
         );
 
-        newMasterPatches.reserve(patchNames.size() + 1);
-        newSlavePatches.reserve(patchNames.size() + 1);
-        forAll(patchNames, i)
+        dict.lookup("internalFacesOnly") >> internalFacesOnly;
+        noFields = dict.lookupOrDefault("noFields", false);
+
+        const dictionary& selectionsDict = dict.subDict("baffles");
+
+        label n = 0;
+        forAllConstIter(dictionary, selectionsDict, iter)
         {
-            newMasterPatches.append(findPatchID(mesh, patchNames[i][0]));
-            newSlavePatches.append(findPatchID(mesh, patchNames[i][1]));
-            Info<< "Using additional patches " << patchNames[i]
-                << " at indices " << newMasterPatches.last()
-                << " and " << newSlavePatches.last()
-                << endl;
+            if (iter().isDict())
+            {
+                n++;
+            }
+        }
+        selectors.setSize(n);
+        n = 0;
+        forAllConstIter(dictionary, selectionsDict, iter)
+        {
+            if (iter().isDict())
+            {
+                selectors.set
+                (
+                    n++,
+                    faceSelection::New(iter().keyword(), mesh, iter().dict())
+                );
+            }
         }
     }
 
 
-    const bool overwrite = args.optionFound("overwrite");
-    const bool internalFacesOnly = args.optionFound("internalFacesOnly");
-
     if (internalFacesOnly)
     {
         Info<< "Not converting faces on non-coupled patches." << nl << endl;
@@ -237,41 +205,230 @@ int main(int argc, char *argv[])
     // Read vol fields.
     Info<< "Reading geometric fields" << nl << endl;
     PtrList<volScalarField> vsFlds;
-    ReadFields(mesh, objects, vsFlds);
+    if (!noFields) ReadFields(mesh, objects, vsFlds);
 
     PtrList<volVectorField> vvFlds;
-    ReadFields(mesh, objects, vvFlds);
+    if (!noFields) ReadFields(mesh, objects, vvFlds);
 
     PtrList<volSphericalTensorField> vstFlds;
-    ReadFields(mesh, objects, vstFlds);
+    if (!noFields) ReadFields(mesh, objects, vstFlds);
 
     PtrList<volSymmTensorField> vsymtFlds;
-    ReadFields(mesh, objects, vsymtFlds);
+    if (!noFields) ReadFields(mesh, objects, vsymtFlds);
 
     PtrList<volTensorField> vtFlds;
-    ReadFields(mesh, objects, vtFlds);
+    if (!noFields) ReadFields(mesh, objects, vtFlds);
 
     // Read surface fields.
 
     PtrList<surfaceScalarField> ssFlds;
-    ReadFields(mesh, objects, ssFlds);
+    if (!noFields) ReadFields(mesh, objects, ssFlds);
 
     PtrList<surfaceVectorField> svFlds;
-    ReadFields(mesh, objects, svFlds);
+    if (!noFields) ReadFields(mesh, objects, svFlds);
 
     PtrList<surfaceSphericalTensorField> sstFlds;
-    ReadFields(mesh, objects, sstFlds);
+    if (!noFields) ReadFields(mesh, objects, sstFlds);
 
     PtrList<surfaceSymmTensorField> ssymtFlds;
-    ReadFields(mesh, objects, ssymtFlds);
+    if (!noFields) ReadFields(mesh, objects, ssymtFlds);
 
     PtrList<surfaceTensorField> stFlds;
-    ReadFields(mesh, objects, stFlds);
+    if (!noFields) ReadFields(mesh, objects, stFlds);
+
+
+
+
+    // Creating (if necessary) baffles
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    forAll(selectors, selectorI)
+    {
+        const word& name = selectors[selectorI].name();
+
+        if (mesh.faceZones().findZoneID(name) == -1)
+        {
+            mesh.faceZones().clearAddressing();
+            label sz = mesh.faceZones().size();
+
+            labelList addr(0);
+            boolList flip(0);
+            mesh.faceZones().setSize(sz+1);
+            mesh.faceZones().set
+            (
+                sz,
+                new faceZone(name, addr, flip, sz, mesh.faceZones())
+            );
+        }
+    }
+
+
+    // Select faces
+    // ~~~~~~~~~~~~
+
+    //- Per face zoneID it is in and flip status.
+    labelList faceToZoneID(mesh.nFaces(), -1);
+    boolList faceToFlip(mesh.nFaces(), false);
+    forAll(selectors, selectorI)
+    {
+        const word& name = selectors[selectorI].name();
+        label zoneID = mesh.faceZones().findZoneID(name);
+
+        selectors[selectorI].select(zoneID, faceToZoneID, faceToFlip);
+    }
+
+    // Add faces to faceZones
+    labelList nFaces(mesh.faceZones().size(), 0);
+    forAll(faceToZoneID, faceI)
+    {
+        label zoneID = faceToZoneID[faceI];
+        if (zoneID != -1)
+        {
+            nFaces[zoneID]++;
+        }
+    }
+
+    forAll(selectors, selectorI)
+    {
+        const word& name = selectors[selectorI].name();
+        label zoneID = mesh.faceZones().findZoneID(name);
+
+        label& n = nFaces[zoneID];
+        labelList addr(n);
+        boolList flip(n);
+        n = 0;
+        forAll(faceToZoneID, faceI)
+        {
+            label zone = faceToZoneID[faceI];
+            if (zone == zoneID)
+            {
+                addr[n] = faceI;
+                flip[n] = faceToFlip[faceI];
+                n++;
+            }
+        }
+
+        Info<< "Created zone " << name
+            << " at index " << zoneID
+            << " with " << n << " faces" << endl;
+
+        mesh.faceZones().set
+        (
+            zoneID,
+            new faceZone(name, addr, flip, zoneID, mesh.faceZones())
+        );
+    }
+
+
+
+    // Count patches to add
+    // ~~~~~~~~~~~~~~~~~~~~
+    HashSet<word> bafflePatches;
+    {
+        forAll(selectors, selectorI)
+        {
+            const dictionary& patchSources
+            (
+                selectors[selectorI].dict().subDict("patches")
+            );
+            forAllConstIter(dictionary, patchSources, iter)
+            {
+                //const word& patchName = iter().keyword();
+                const word patchName(iter().dict()["name"]);
+                bafflePatches.insert(patchName);
+            }
+        }
+    }
+
+
+
+    // Create baffles
+    // ~~~~~~~~~~~~~~
+    // Is done in multiple steps
+    // - create patches with 'calculated' patchFields
+    // - move faces into these patches
+    // - change the patchFields to the wanted type
+    // This order is done so e.g. fixedJump works:
+    // - you cannot create patchfields at the same time as patches since
+    //   they do an evaluate upon construction
+    // - you want to create the patchField only after you have faces
+    //   so you don't get the 'create-from-nothing' mapping problem.
+
+
+    // Pass 1: add patches
+    // ~~~~~~~~~~~~~~~~~~~
+
+    //HashSet<word> addedPatches;
+    {
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+        forAll(selectors, selectorI)
+        {
+            const dictionary& patchSources
+            (
+                selectors[selectorI].dict().subDict("patches")
+            );
+            forAllConstIter(dictionary, patchSources, iter)
+            {
+                //const word& patchName = iter().keyword();
+                const word patchName(iter().dict()["name"]);
+
+                label destPatchI = pbm.findPatchID(patchName);
+
+                if (destPatchI == -1)
+                {
+                    dictionary patchDict = iter().dict();
+                    patchDict.set("nFaces", 0);
+                    patchDict.set("startFace", 0);
+
+                    Info<< "Adding new patch " << patchName
+                        << " from " << patchDict << endl;
+
+                    autoPtr<polyPatch> ppPtr
+                    (
+                        polyPatch::New
+                        (
+                            patchName,
+                            patchDict,
+                            0,
+                            pbm
+                        )
+                    );
+
+                    // Add patch, create calculated everywhere
+                    fvMeshTools::addPatch
+                    (
+                        mesh,
+                        ppPtr(),
+                        dictionary(),   // do not set specialised patchFields
+                        calculatedFvPatchField<scalar>::typeName,
+                        true            // parallel sync'ed addition
+                    );
+
+                    //addedPatches.insert(patchName);
+                }
+                else
+                {
+                    Info<< "Patch '" << patchName << "' already exists.  Only "
+                        << "moving patch faces - type will remain the same"
+                        << endl;
+                }
+            }
+        }
+    }
+
+
+
+    // Make sure patches and zoneFaces are synchronised across couples
+    mesh.boundaryMesh().checkParallelSync(true);
+    mesh.faceZones().checkParallelSync(true);
+
 
 
     // Mesh change container
     polyTopoChange meshMod(mesh);
 
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
 
     // Do the actual changes. Note:
     // - loop in incrementing face order (not necessary if faceZone ordered).
@@ -284,161 +441,197 @@ int main(int argc, char *argv[])
     PackedBoolList modifiedFace(mesh.nFaces());
     label nModified = 0;
 
-    forAll(newMasterPatches, i)
+    forAll(selectors, selectorI)
     {
-        // Pass 1. Do selected side of zone
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        const word& name = selectors[selectorI].name();
+        label zoneID = mesh.faceZones().findZoneID(name);
+        const faceZone& fZone = mesh.faceZones()[zoneID];
 
-        for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+        const dictionary& patchSources
+        (
+            selectors[selectorI].dict().subDict("patches")
+        );
+
+        DynamicList<label> newMasterPatches(patchSources.size());
+        DynamicList<label> newSlavePatches(patchSources.size());
+
+        bool master = true;
+
+        forAllConstIter(dictionary, patchSources, iter)
         {
-            label zoneFaceI = fZone.whichFace(faceI);
-
-            if (zoneFaceI != -1)
+            //const word& patchName = iter().keyword();
+            const word patchName(iter().dict()["name"]);
+            label patchI = pbm.findPatchID(patchName);
+            if (master)
             {
-                if (!fZone.flipMap()[zoneFaceI])
-                {
-                    // Use owner side of face
-                    modifyOrAddFace
-                    (
-                        meshMod,
-                        mesh.faces()[faceI],    // modified face
-                        faceI,                  // label of face
-                        mesh.faceOwner()[faceI],// owner
-                        false,                  // face flip
-                        newMasterPatches[i],    // patch for face
-                        zoneID.index(),         // zone for face
-                        false,                  // face flip in zone
-                        modifiedFace            // modify or add status
-                    );
-                }
-                else
-                {
-                    // Use neighbour side of face
-                    modifyOrAddFace
-                    (
-                        meshMod,
-                        mesh.faces()[faceI].reverseFace(),  // modified face
-                        faceI,                      // label of face
-                        mesh.faceNeighbour()[faceI],// owner
-                        true,                       // face flip
-                        newMasterPatches[i],        // patch for face
-                        zoneID.index(),             // zone for face
-                        true,                       // face flip in zone
-                        modifiedFace                // modify or add status
-                    );
-                }
-
-                nModified++;
+                newMasterPatches.append(patchI);
             }
+            else
+            {
+                newSlavePatches.append(patchI);
+            }
+            master = !master;
         }
 
 
-        // Pass 2. Do other side of zone
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+        forAll(newMasterPatches, i)
         {
-            label zoneFaceI = fZone.whichFace(faceI);
+            // Pass 1. Do selected side of zone
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-            if (zoneFaceI != -1)
+            for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
             {
-                if (!fZone.flipMap()[zoneFaceI])
+                label zoneFaceI = fZone.whichFace(faceI);
+
+                if (zoneFaceI != -1)
                 {
-                    // Use neighbour side of face
-                    modifyOrAddFace
-                    (
-                        meshMod,
-                        mesh.faces()[faceI].reverseFace(),  // modified face
-                        faceI,                              // label of face
-                        mesh.faceNeighbour()[faceI],        // owner
-                        true,                               // face flip
-                        newSlavePatches[i],                 // patch for face
-                        zoneID.index(),                     // zone for face
-                        true,                               // face flip in zone
-                        modifiedFace                        // modify or add
-                    );
-                }
-                else
-                {
-                    // Use owner side of face
-                    modifyOrAddFace
-                    (
-                        meshMod,
-                        mesh.faces()[faceI],    // modified face
-                        faceI,                  // label of face
-                        mesh.faceOwner()[faceI],// owner
-                        false,                  // face flip
-                        newSlavePatches[i],     // patch for face
-                        zoneID.index(),         // zone for face
-                        false,                  // face flip in zone
-                        modifiedFace            // modify or add status
-                    );
-                }
-            }
-        }
-
-
-        // Modify any boundary faces
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~
-
-        // Normal boundary:
-        // - move to new patch. Might already be back-to-back baffle
-        // you want to add cyclic to. Do warn though.
-        //
-        // Processor boundary:
-        // - do not move to cyclic
-        // - add normal patches though.
-
-        // For warning once per patch.
-        labelHashSet patchWarned;
-
-        forAll(patches, patchI)
-        {
-            const polyPatch& pp = patches[patchI];
-
-            label newPatchI = newMasterPatches[i];
-
-            if (pp.coupled() && patches[newPatchI].coupled())
-            {
-                // Do not allow coupled faces to be moved to different coupled
-                // patches.
-            }
-            else if (pp.coupled() || !internalFacesOnly)
-            {
-                forAll(pp, i)
-                {
-                    label faceI = pp.start()+i;
-
-                    label zoneFaceI = fZone.whichFace(faceI);
-
-                    if (zoneFaceI != -1)
+                    if (!fZone.flipMap()[zoneFaceI])
                     {
-                        if (patchWarned.insert(patchI))
-                        {
-                            WarningIn(args.executable())
-                                << "Found boundary face (in patch " << pp.name()
-                                << ") in faceZone " << fZone.name()
-                                << " to convert to baffle patch "
-                                << patches[newPatchI].name()
-                                << endl
-                                << "    Run with -internalFacesOnly option"
-                                << " if you don't wish to convert"
-                                << " boundary faces." << endl;
-                        }
-
+                        // Use owner side of face
                         modifyOrAddFace
                         (
                             meshMod,
-                            mesh.faces()[faceI],        // modified face
+                            mesh.faces()[faceI],    // modified face
+                            faceI,                  // label of face
+                            mesh.faceOwner()[faceI],// owner
+                            false,                  // face flip
+                            newMasterPatches[i],    // patch for face
+                            fZone.index(),          // zone for face
+                            false,                  // face flip in zone
+                            modifiedFace            // modify or add status
+                        );
+                    }
+                    else
+                    {
+                        // Use neighbour side of face
+                        modifyOrAddFace
+                        (
+                            meshMod,
+                            mesh.faces()[faceI].reverseFace(),  // modified face
                             faceI,                      // label of face
-                            mesh.faceOwner()[faceI],    // owner
-                            false,                      // face flip
-                            newPatchI,                  // patch for face
-                            zoneID.index(),             // zone for face
-                            fZone.flipMap()[zoneFaceI], // face flip in zone
+                            mesh.faceNeighbour()[faceI],// owner
+                            true,                       // face flip
+                            newMasterPatches[i],        // patch for face
+                            fZone.index(),              // zone for face
+                            true,                       // face flip in zone
                             modifiedFace                // modify or add status
                         );
-                        nModified++;
+                    }
+
+                    nModified++;
+                }
+            }
+
+
+            // Pass 2. Do other side of zone
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+            {
+                label zoneFaceI = fZone.whichFace(faceI);
+
+                if (zoneFaceI != -1)
+                {
+                    if (!fZone.flipMap()[zoneFaceI])
+                    {
+                        // Use neighbour side of face
+                        modifyOrAddFace
+                        (
+                            meshMod,
+                            mesh.faces()[faceI].reverseFace(),  // modified face
+                            faceI,                          // label of face
+                            mesh.faceNeighbour()[faceI],    // owner
+                            true,                           // face flip
+                            newSlavePatches[i],             // patch for face
+                            fZone.index(),                  // zone for face
+                            true,                           // face flip in zone
+                            modifiedFace                    // modify or add
+                        );
+                    }
+                    else
+                    {
+                        // Use owner side of face
+                        modifyOrAddFace
+                        (
+                            meshMod,
+                            mesh.faces()[faceI],    // modified face
+                            faceI,                  // label of face
+                            mesh.faceOwner()[faceI],// owner
+                            false,                  // face flip
+                            newSlavePatches[i],     // patch for face
+                            fZone.index(),          // zone for face
+                            false,                  // face flip in zone
+                            modifiedFace            // modify or add status
+                        );
+                    }
+                }
+            }
+
+
+            // Modify any boundary faces
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            // Normal boundary:
+            // - move to new patch. Might already be back-to-back baffle
+            // you want to add cyclic to. Do warn though.
+            //
+            // Processor boundary:
+            // - do not move to cyclic
+            // - add normal patches though.
+
+            // For warning once per patch.
+            labelHashSet patchWarned;
+
+            forAll(pbm, patchI)
+            {
+                const polyPatch& pp = pbm[patchI];
+
+                label newPatchI = newMasterPatches[i];
+
+                if (pp.coupled() && pbm[newPatchI].coupled())
+                {
+                    // Do not allow coupled faces to be moved to different
+                    // coupled patches.
+                }
+                else if (pp.coupled() || !internalFacesOnly)
+                {
+                    forAll(pp, i)
+                    {
+                        label faceI = pp.start()+i;
+
+                        label zoneFaceI = fZone.whichFace(faceI);
+
+                        if (zoneFaceI != -1)
+                        {
+                            if (patchWarned.insert(patchI))
+                            {
+                                WarningIn(args.executable())
+                                    << "Found boundary face (in patch "
+                                    << pp.name()
+                                    << ") in faceZone " << fZone.name()
+                                    << " to convert to baffle patch "
+                                    << pbm[newPatchI].name()
+                                    << endl
+                                    << "    Run with -internalFacesOnly option"
+                                    << " if you don't wish to convert"
+                                    << " boundary faces." << endl;
+                            }
+
+                            modifyOrAddFace
+                            (
+                                meshMod,
+                                mesh.faces()[faceI],        // modified face
+                                faceI,                      // label of face
+                                mesh.faceOwner()[faceI],    // owner
+                                false,                      // face flip
+                                newPatchI,                  // patch for face
+                                fZone.index(),              // zone for face
+                                fZone.flipMap()[zoneFaceI], // face flip in zone
+                                modifiedFace                // modify or add
+                            );
+                            nModified++;
+                        }
                     }
                 }
             }
@@ -447,7 +640,8 @@ int main(int argc, char *argv[])
 
 
     Info<< "Converted " << returnReduce(nModified, sumOp<label>())
-        << " faces into boundary faces on patches " << patchNames << nl << endl;
+        << " faces into boundary faces in patches "
+        << bafflePatches.sortedToc() << nl << endl;
 
     if (!overwrite)
     {
@@ -460,14 +654,20 @@ int main(int argc, char *argv[])
     // Update fields
     mesh.updateMesh(map);
 
+
+
     // Correct boundary faces mapped-out-of-nothing.
+    // This is just a hack to correct the value field.
     {
         fvMeshMapper mapper(mesh, map);
         bool hasWarned = false;
-        forAll(newMasterPatches, i)
+
+        forAllConstIter(HashSet<word>, bafflePatches, iter)
         {
-            label patchI = newMasterPatches[i];
+            label patchI = mesh.boundaryMesh().findPatchID(iter.key());
+
             const fvPatchMapper& pm = mapper.boundaryMap()[patchI];
+
             if (pm.sizeBeforeMapping() == 0)
             {
                 if (!hasWarned)
@@ -478,40 +678,48 @@ int main(int argc, char *argv[])
                         << "You might have to edit these fields." << endl;
                 }
 
-                SetPatchFields(vsFlds, patchI, pTraits<scalar>::zero);
-                SetPatchFields(vvFlds, patchI, pTraits<vector>::zero);
-                SetPatchFields(vstFlds, patchI, pTraits<sphericalTensor>::zero);
-                SetPatchFields(vsymtFlds, patchI, pTraits<symmTensor>::zero);
-                SetPatchFields(vtFlds, patchI, pTraits<tensor>::zero);
-
-                SetPatchFields(ssFlds, patchI, pTraits<scalar>::zero);
-                SetPatchFields(svFlds, patchI, pTraits<vector>::zero);
-                SetPatchFields(sstFlds, patchI, pTraits<sphericalTensor>::zero);
-                SetPatchFields(ssymtFlds, patchI, pTraits<symmTensor>::zero);
-                SetPatchFields(stFlds, patchI, pTraits<tensor>::zero);
-            }
-        }
-        forAll(newSlavePatches, i)
-        {
-            label patchI = newSlavePatches[i];
-            const fvPatchMapper& pm = mapper.boundaryMap()[patchI];
-            if (pm.sizeBeforeMapping() == 0)
-            {
-                SetPatchFields(vsFlds, patchI, pTraits<scalar>::zero);
-                SetPatchFields(vvFlds, patchI, pTraits<vector>::zero);
-                SetPatchFields(vstFlds, patchI, pTraits<sphericalTensor>::zero);
-                SetPatchFields(vsymtFlds, patchI, pTraits<symmTensor>::zero);
-                SetPatchFields(vtFlds, patchI, pTraits<tensor>::zero);
-
-                SetPatchFields(ssFlds, patchI, pTraits<scalar>::zero);
-                SetPatchFields(svFlds, patchI, pTraits<vector>::zero);
-                SetPatchFields(sstFlds, patchI, pTraits<sphericalTensor>::zero);
-                SetPatchFields(ssymtFlds, patchI, pTraits<symmTensor>::zero);
-                SetPatchFields(stFlds, patchI, pTraits<tensor>::zero);
+                fvMeshTools::zeroPatchFields(mesh, patchI);
             }
         }
     }
 
+
+    // Pass 2: change patchFields
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    {
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+        forAll(selectors, selectorI)
+        {
+            const dictionary& patchSources
+            (
+                selectors[selectorI].dict().subDict("patches")
+            );
+            forAllConstIter(dictionary, patchSources, iter)
+            {
+                //const word& patchName = iter().keyword();
+                const word patchName(iter().dict()["name"]);
+                label patchI = pbm.findPatchID(patchName);
+
+                if (iter().dict().found("patchFields"))
+                {
+                    const dictionary& patchFieldsDict = iter().dict().subDict
+                    (
+                        "patchFields"
+                    );
+
+                    fvMeshTools::setPatchFields
+                    (
+                        mesh,
+                        patchI,
+                        patchFieldsDict
+                    );
+                }
+            }
+        }
+    }
+
+
     // Move mesh (since morphing might not do this)
     if (map().hasMotionPoints())
     {
@@ -522,6 +730,7 @@ int main(int argc, char *argv[])
     {
         mesh.setInstance(oldInstance);
     }
+
     Info<< "Writing mesh to " << runTime.timeName() << endl;
 
     mesh.write();
diff --git a/applications/utilities/mesh/manipulation/createBaffles/createBafflesDict b/applications/utilities/mesh/manipulation/createBaffles/createBafflesDict
new file mode 100644
index 0000000000..99c03b4cfc
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/createBaffles/createBafflesDict
@@ -0,0 +1,189 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      createBafflesDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Whether to convert internal faces only (so leave boundary faces intact).
+// This is only relevant if your face selection type can pick up boundary
+// faces.
+internalFacesOnly true;
+
+// Optionally do not read/convert/write any fields.
+//noFields true;
+
+
+// Baffles to create.
+baffles
+{
+    baffleFaces
+    {
+        //- Use predefined faceZone to select faces and orientation.
+        type        faceZone;
+        zoneName    baffleFaces;
+
+
+        //- Optional flip
+        //flip        false;
+
+        patches
+        {
+            master
+            {
+                //- Master side patch
+                name            baffles;
+                type            wall;
+
+                //- Optional override of added patchfields. If not specified
+                //  any added patchfields are of type calculated.
+                patchFields
+                {
+                    epsilon
+                    {
+                        type            epsilonWallFunction;
+                        Cmu             0.09;
+                        kappa           0.41;
+                        E               9.8;
+                        value           uniform 0;
+                    }
+                    k
+                    {
+                        type            kqRWallFunction;
+                        value           uniform 0;
+                    }
+                    nut
+                    {
+                        type            nutkWallFunction;
+                        Cmu             0.09;
+                        kappa           0.41;
+                        E               9.8;
+                        value           uniform 0;
+                    }
+                    nuTilda
+                    {
+                        type            zeroGradient;
+                    }
+                    p
+                    {
+                        type            zeroGradient;
+                    }
+                    U
+                    {
+                        type            fixedValue;
+                        value           uniform (0 0 0);
+                    }
+                }
+            }
+            slave
+            {
+                //- Slave side patch
+                name            baffles;
+                type            wall;
+
+                patchFields
+                {
+                    epsilon
+                    {
+                        type            epsilonWallFunction;
+                        Cmu             0.09;
+                        kappa           0.41;
+                        E               9.8;
+                        value           uniform 0;
+                    }
+                    k
+                    {
+                        type            kqRWallFunction;
+                        value           uniform 0;
+                    }
+                    nut
+                    {
+                        type            nutkWallFunction;
+                        Cmu             0.09;
+                        kappa           0.41;
+                        E               9.8;
+                        value           uniform 0;
+                    }
+                    nuTilda
+                    {
+                        type            zeroGradient;
+                    }
+                    p
+                    {
+                        type            zeroGradient;
+                    }
+                    U
+                    {
+                        type            fixedValue;
+                        value           uniform (0 0 0);
+                    }
+                }
+            }
+        }
+    }
+
+
+
+    cyclicFaces
+    {
+        //- Select faces and orientation through a searchableSurface
+        type        searchableSurface;
+        surface     searchablePlate;
+        origin      (0.099 -0.006 0.004);
+        span        (0 0.012 0.012);
+
+        patches
+        {
+            master
+            {
+                //- Master side patch
+
+                name            fan_half0;
+                type            cyclic;
+                neighbourPatch  fan_half1;
+
+                patchFields
+                {
+                    p
+                    {
+                        type            fan;
+                        patchType       cyclic;
+                        jump            uniform 0;
+                        value           uniform 0;
+                        jumpTable       polynomial 1((100 0));
+                    }
+                }
+            }
+            slave
+            {
+                //- Slave side patch
+
+                name            fan_half1;
+                type            cyclic;
+                neighbourPatch  fan_half0;
+
+                patchFields
+                {
+                    p
+                    {
+                        type            fan;
+                        patchType       cyclic;
+                        value           uniform 0;
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/createBaffles/faceSelection/faceSelection.C b/applications/utilities/mesh/manipulation/createBaffles/faceSelection/faceSelection.C
new file mode 100644
index 0000000000..57a6b554c7
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/createBaffles/faceSelection/faceSelection.C
@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "faceSelection.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(faceSelection, 0);
+    defineRunTimeSelectionTable(faceSelection, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::faceSelection::faceSelection
+(
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    name_(name),
+    mesh_(mesh),
+    dict_(dict),
+    flip_(dict.lookupOrDefault("flip", false))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::faceSelection::~faceSelection()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::faceSelection> Foam::faceSelection::New
+(
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+{
+    const word sampleType(dict.lookup("type"));
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(sampleType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "faceSelection::New"
+            "(const word&, const fvMesh&, const dictionary&)"
+        )   << "Unknown faceSelection type "
+            << sampleType << nl << nl
+            << "Valid faceSelection types : " << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<faceSelection>(cstrIter()(name, mesh, dict));
+}
+
+
+void Foam::faceSelection::select
+(
+    const label zoneID,
+    labelList& faceToZoneID,
+    boolList& faceToFlip
+) const
+{
+    if (flip_)
+    {
+        forAll(faceToZoneID, faceI)
+        {
+            if (faceToZoneID[faceI] == zoneID)
+            {
+                faceToFlip[faceI] = !faceToFlip[faceI];
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/createBaffles/faceSelection/faceSelection.H b/applications/utilities/mesh/manipulation/createBaffles/faceSelection/faceSelection.H
new file mode 100644
index 0000000000..7be333f030
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/createBaffles/faceSelection/faceSelection.H
@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::faceSelection
+
+Description
+    Face selection method for createBaffles
+
+SourceFiles
+    faceSelection.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faceSelection_H
+#define faceSelection_H
+
+#include "dictionary.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
+#include "boolList.H"
+#include "labelList.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class fvMesh;
+
+/*---------------------------------------------------------------------------*\
+                         Class faceSelection Declaration
+\*---------------------------------------------------------------------------*/
+
+class faceSelection
+{
+protected:
+
+    // Protected data
+
+        //- Name
+        const word name_;
+
+        //- Reference to mesh
+        const fvMesh& mesh_;
+
+        //- Input dictionary
+        const dictionary dict_;
+
+        //- Switch direction?
+        const Switch flip_;
+
+public:
+
+    //- Runtime type information
+    TypeName("faceSelection");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            faceSelection,
+            dictionary,
+            (
+                const word& name,
+                const fvMesh& mesh,
+                const dictionary& dict
+            ),
+            (name, mesh, dict)
+        );
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        faceSelection
+        (
+            const word& name,
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Clone
+        autoPtr<faceSelection> clone() const
+        {
+            notImplemented("autoPtr<faceSelection> clone() const");
+            return autoPtr<faceSelection>(NULL);
+        }
+
+
+    // Selectors
+
+        //- Return a reference to the selected faceSelection
+        static autoPtr<faceSelection> New
+        (
+            const word& name,
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~faceSelection();
+
+
+    // Member Functions
+
+        const word& name() const
+        {
+            return name_;
+        }
+
+        const dictionary& dict() const
+        {
+            return dict_;
+        }
+
+        virtual void select(const label, labelList&, boolList&) const = 0;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/createBaffles/faceSelection/faceZoneSelection.C b/applications/utilities/mesh/manipulation/createBaffles/faceSelection/faceZoneSelection.C
new file mode 100644
index 0000000000..bdc6d84e67
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/createBaffles/faceSelection/faceZoneSelection.C
@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "faceZoneSelection.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace faceSelections
+{
+    defineTypeNameAndDebug(faceZoneSelection, 0);
+    addToRunTimeSelectionTable(faceSelection, faceZoneSelection, dictionary);
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::faceSelections::faceZoneSelection::faceZoneSelection
+(
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    faceSelection(name, mesh, dict),
+    zoneName_(dict_.lookup("zoneName"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::faceSelections::faceZoneSelection::~faceZoneSelection()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::faceSelections::faceZoneSelection::select
+(
+    const label zoneID,
+    labelList& faceToZoneID,
+    boolList& faceToFlip
+) const
+{
+    label readID = mesh_.faceZones().findZoneID(zoneName_);
+
+    if (readID == -1)
+    {
+        FatalErrorIn
+        (
+            "faceSelections::faceZoneSelection::select(labelList&) const"
+        )   << "Cannot find faceZone " << zoneName_ << nl << "Valid zones are "
+            << mesh_.faceZones().names()
+            << exit(FatalError);
+    }
+
+    const faceZone& fZone = mesh_.faceZones()[readID];
+
+    forAll(fZone, i)
+    {
+        label faceI = fZone[i];
+
+        if (faceToZoneID[faceI] == -1)
+        {
+            faceToZoneID[faceI] = zoneID;
+            faceToFlip[faceI] = fZone.flipMap()[i];
+        }
+        else if (faceToZoneID[faceI] != zoneID)
+        {
+            FatalErrorIn
+            (
+                "faceSelections::faceZoneSelection::select(labelList&) const"
+            )   << "Face " << faceI << " already in faceZone "
+                << faceToZoneID[faceI]
+                << exit(FatalError);
+        }
+    }
+
+    faceSelection::select(zoneID, faceToZoneID, faceToFlip);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/createBaffles/faceSelection/faceZoneSelection.H b/applications/utilities/mesh/manipulation/createBaffles/faceSelection/faceZoneSelection.H
new file mode 100644
index 0000000000..1a2b15d366
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/createBaffles/faceSelection/faceZoneSelection.H
@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::faceSelections::faceZoneSelection
+
+Description
+    Deselect cells not reachable from 'inside' points
+
+SourceFiles
+    faceZoneSelection.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faceZoneSelection_H
+#define faceZoneSelection_H
+
+#include "faceSelection.H"
+#include "pointField.H"
+#include "boolList.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class regionSplit;
+
+namespace faceSelections
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class faceZoneSelection Declaration
+\*---------------------------------------------------------------------------*/
+
+class faceZoneSelection
+:
+    public faceSelection
+{
+    // Private data
+
+        //- Name of faceZone
+        const word zoneName_;
+
+
+   // Private Member Functions
+
+public:
+
+    //- Runtime type information
+    TypeName("faceZone");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        faceZoneSelection
+        (
+            const word& name,
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Clone
+        autoPtr<faceSelection> clone() const
+        {
+            notImplemented("autoPtr<faceSelection> clone() const");
+            return autoPtr<faceSelection>(NULL);
+        }
+
+
+    //- Destructor
+    virtual ~faceZoneSelection();
+
+
+    // Member Functions
+
+        //- Apply this selector
+        virtual void select(const label zoneID, labelList&, boolList&) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace faceSelections
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/createBaffles/faceSelection/searchableSurfaceSelection.C b/applications/utilities/mesh/manipulation/createBaffles/faceSelection/searchableSurfaceSelection.C
new file mode 100644
index 0000000000..1dbaaafaf8
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/createBaffles/faceSelection/searchableSurfaceSelection.C
@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "searchableSurfaceSelection.H"
+#include "addToRunTimeSelectionTable.H"
+#include "syncTools.H"
+#include "searchableSurface.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace faceSelections
+{
+    defineTypeNameAndDebug(searchableSurfaceSelection, 0);
+    addToRunTimeSelectionTable
+    (
+        faceSelection,
+        searchableSurfaceSelection,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::faceSelections::searchableSurfaceSelection::searchableSurfaceSelection
+(
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    faceSelection(name, mesh, dict),
+    surfacePtr_
+    (
+        searchableSurface::New
+        (
+            word(dict.lookup("surface")),
+            mesh.objectRegistry::db(),
+            dict
+        )
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::faceSelections::searchableSurfaceSelection::~searchableSurfaceSelection()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::faceSelections::searchableSurfaceSelection::select
+(
+    const label zoneID,
+    labelList& faceToZoneID,
+    boolList& faceToFlip
+) const
+{
+    // Get cell-cell centre vectors
+
+    pointField start(mesh_.nFaces());
+    pointField end(mesh_.nFaces());
+
+    // Internal faces
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        start[faceI] = mesh_.cellCentres()[mesh_.faceOwner()[faceI]];
+        end[faceI] = mesh_.cellCentres()[mesh_.faceNeighbour()[faceI]];
+    }
+
+    // Boundary faces
+    vectorField neighbourCellCentres;
+    syncTools::swapBoundaryCellList
+    (
+        mesh_,
+        mesh_.cellCentres(),
+        neighbourCellCentres
+    );
+
+    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] = mesh_.cellCentres()[mesh_.faceOwner()[faceI]];
+                end[faceI] = neighbourCellCentres[faceI-mesh_.nInternalFaces()];
+            }
+        }
+        else
+        {
+            forAll(pp, i)
+            {
+                label faceI = pp.start()+i;
+                start[faceI] = mesh_.cellCentres()[mesh_.faceOwner()[faceI]];
+                end[faceI] = mesh_.faceCentres()[faceI];
+            }
+        }
+    }
+
+    List<pointIndexHit> hits;
+    surfacePtr_().findLine(start, end, hits);
+    pointField normals;
+    surfacePtr_().getNormal(hits, normals);
+
+    //- Note: do not select boundary faces.
+
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        if (hits[faceI].hit())
+        {
+            faceToZoneID[faceI] = zoneID;
+            vector d = end[faceI]-start[faceI];
+            faceToFlip[faceI] = ((normals[faceI] & d) < 0);
+        }
+    }
+    forAll(pbm, patchI)
+    {
+        const polyPatch& pp = pbm[patchI];
+
+        if (pp.coupled())
+        {
+            forAll(pp, i)
+            {
+                label faceI = pp.start()+i;
+                if (hits[faceI].hit())
+                {
+                    faceToZoneID[faceI] = zoneID;
+                    vector d = end[faceI]-start[faceI];
+                    faceToFlip[faceI] = ((normals[faceI] & d) < 0);
+                }
+            }
+        }
+    }
+
+    faceSelection::select(zoneID, faceToZoneID, faceToFlip);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/createBaffles/faceSelection/searchableSurfaceSelection.H b/applications/utilities/mesh/manipulation/createBaffles/faceSelection/searchableSurfaceSelection.H
new file mode 100644
index 0000000000..f983132552
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/createBaffles/faceSelection/searchableSurfaceSelection.H
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::faceSelections::searchableSurfaceSelection
+
+Description
+    Selects all (internal or coupled) faces intersecting the searchableSurface.
+
+SourceFiles
+    searchableSurfaceSelection.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchableSurfaceSelection_H
+#define searchableSurfaceSelection_H
+
+#include "faceSelection.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class searchableSurface;
+
+namespace faceSelections
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class searchableSurfaceSelection Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchableSurfaceSelection
+:
+    public faceSelection
+{
+    // Private data
+
+        autoPtr<searchableSurface> surfacePtr_;
+
+public:
+
+    //- Runtime type information
+    TypeName("searchableSurface");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        searchableSurfaceSelection
+        (
+            const word& name,
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Clone
+        autoPtr<faceSelection> clone() const
+        {
+            notImplemented("autoPtr<faceSelection> clone() const");
+            return autoPtr<faceSelection>(NULL);
+        }
+
+
+    //- Destructor
+    virtual ~searchableSurfaceSelection();
+
+
+    // Member Functions
+
+        virtual void select(const label zoneID, labelList&, boolList&) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace faceSelections
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
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/etc/controlDict b/etc/controlDict
index a34d4a1977..a47741c94c 100644
--- a/etc/controlDict
+++ b/etc/controlDict
@@ -630,7 +630,7 @@ DebugSwitches
     muSgsSpalartAllmarasWallFunction 0;
     multiDirRefinement  0;
     multiHoleInjector   0;
-    multiLevel          1;
+    multiLevel          0;
     multivariateSelection 0;
     mutRoughWallFunction 0;
     mutSpalartAllmarasStandardRoughWallFunction  0;
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..e4794b634b 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,29 @@ 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(),   // optional per field patchField
+                            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 +975,30 @@ 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 +2226,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 +2281,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 239671d58c..ae9d512529 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..682a96422e
--- /dev/null
+++ b/src/dynamicMesh/fvMeshTools/fvMeshTools.C
@@ -0,0 +1,355 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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;
+    }
+
+
+    // Append at end unless there are processor patches
+    label insertPatchI = polyPatches.size();
+    label startFaceI = mesh.nFaces();
+
+    if (!isA<processorPolyPatch>(patch))
+    {
+        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..0e7e132d9c
--- /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"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+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. Inserts patch before all processor patches.
+    //  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/finiteVolume/Make/files b/src/finiteVolume/Make/files
index b982734e8b..5e097730fd 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -342,6 +342,7 @@ $(limitedGradSchemes)/cellMDLimitedGrad/cellMDLimitedGrads.C
 snGradSchemes = finiteVolume/snGradSchemes
 $(snGradSchemes)/snGradScheme/snGradSchemes.C
 $(snGradSchemes)/correctedSnGrad/correctedSnGrads.C
+$(snGradSchemes)/faceCorrectedSnGrad/faceCorrectedSnGrads.C
 $(snGradSchemes)/limitedSnGrad/limitedSnGrads.C
 $(snGradSchemes)/uncorrectedSnGrad/uncorrectedSnGrads.C
 $(snGradSchemes)/orthogonalSnGrad/orthogonalSnGrads.C
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.C b/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.C
new file mode 100644
index 0000000000..6b04bb2efb
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.C
@@ -0,0 +1,167 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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 "faceCorrectedSnGrad.H"
+#include "volPointInterpolation.H"
+#include "triangle.H"
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::fv::faceCorrectedSnGrad<Type>::~faceCorrectedSnGrad()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
+Foam::fv::faceCorrectedSnGrad<Type>::fullGradCorrection
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    const fvMesh& mesh = this->mesh();
+
+    GeometricField<Type, pointPatchField, pointMesh> pvf
+    (
+        volPointInterpolation::New(mesh).interpolate(vf)
+    );
+
+    // construct GeometricField<Type, fvsPatchField, surfaceMesh>
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsfCorr
+    (
+        new GeometricField<Type, fvsPatchField, surfaceMesh>
+        (
+            IOobject
+            (
+                "snGradCorr("+vf.name()+')',
+                vf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            vf.dimensions()*mesh.nonOrthDeltaCoeffs().dimensions()
+        )
+    );
+
+    Field<Type>& sfCorr = tsfCorr().internalField();
+
+    const pointField& points = mesh.points();
+    const faceList& faces = mesh.faces();
+    const vectorField& Sf = mesh.Sf().internalField();
+    const vectorField& C = mesh.C().internalField();
+    const scalarField& magSf = mesh.magSf().internalField();
+    const labelList& owner = mesh.owner();
+    const labelList& neighbour = mesh.neighbour();
+
+    forAll(sfCorr, facei)
+    {
+        typename outerProduct<vector, Type>::type fgrad
+        (
+            outerProduct<vector, Type>::type::zero
+        );
+
+        const face& fi = faces[facei];
+
+        vector nf(Sf[facei]/magSf[facei]);
+
+        for (label pi=0; pi<fi.size(); pi++)
+        {
+            // Next point index
+            label pj = (pi+1)%fi.size();
+
+            // Edge normal in plane of face
+            vector edgen(nf^(points[fi[pj]] - points[fi[pi]]));
+
+            // Edge centre field value
+            Type pvfe(0.5*(pvf[fi[pj]] + pvf[fi[pi]]));
+
+            // Integrate face gradient
+            fgrad += edgen*pvfe;
+        }
+
+        // Finalize face-gradient by dividing by face area
+        fgrad /= magSf[facei];
+
+        // Calculate correction vector
+        vector dCorr(C[neighbour[facei]] - C[owner[facei]]);
+        dCorr /= (nf & dCorr);
+
+        // if (mag(dCorr) > 2) dCorr *= 2/mag(dCorr);
+
+        sfCorr[facei] = dCorr&fgrad;
+    }
+
+    tsfCorr().boundaryField() = pTraits<Type>::zero;
+
+    return tsfCorr;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
+Foam::fv::faceCorrectedSnGrad<Type>::correction
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    const fvMesh& mesh = this->mesh();
+
+    // construct GeometricField<Type, fvsPatchField, surfaceMesh>
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tssf
+    (
+        new GeometricField<Type, fvsPatchField, surfaceMesh>
+        (
+            IOobject
+            (
+                "snGradCorr("+vf.name()+')',
+                vf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            vf.dimensions()*mesh.nonOrthDeltaCoeffs().dimensions()
+        )
+    );
+    GeometricField<Type, fvsPatchField, surfaceMesh>& ssf = tssf();
+
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    {
+        ssf.replace
+        (
+            cmpt,
+            faceCorrectedSnGrad<typename pTraits<Type>::cmptType>(mesh)
+           .fullGradCorrection(vf.component(cmpt))
+        );
+    }
+
+    return tssf;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.H b/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.H
new file mode 100644
index 0000000000..a13dd62d18
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.H
@@ -0,0 +1,157 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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::fv::faceCorrectedSnGrad
+
+Description
+    Simple central-difference snGrad scheme with non-orthogonal correction.
+
+SourceFiles
+    faceCorrectedSnGrad.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faceCorrectedSnGrad_H
+#define faceCorrectedSnGrad_H
+
+#include "snGradScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class faceCorrectedSnGrad Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class faceCorrectedSnGrad
+:
+    public snGradScheme<Type>
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise assignment
+        void operator=(const faceCorrectedSnGrad&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("faceCorrected");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        faceCorrectedSnGrad(const fvMesh& mesh)
+        :
+            snGradScheme<Type>(mesh)
+        {}
+
+
+        //- Construct from mesh and data stream
+        faceCorrectedSnGrad(const fvMesh& mesh, Istream&)
+        :
+            snGradScheme<Type>(mesh)
+        {}
+
+
+    //- Destructor
+    virtual ~faceCorrectedSnGrad();
+
+
+    // Member Functions
+
+        //- Return the interpolation weighting factors for the given field
+        virtual tmp<surfaceScalarField> deltaCoeffs
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const
+        {
+            return this->mesh().nonOrthDeltaCoeffs();
+        }
+
+        //- Return true if this scheme uses an explicit correction
+        virtual bool corrected() const
+        {
+            return true;
+        }
+
+        //- Return the explicit correction to the faceCorrectedSnGrad
+        //  for the given field using the gradient of the field
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+        fullGradCorrection
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
+        //- Return the explicit correction to the faceCorrectedSnGrad
+        //  for the given field using the gradients of the field components
+        virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+        correction(const GeometricField<Type, fvPatchField, volMesh>&) const;
+};
+
+
+// * * * * * * * * Template Member Function Specialisations  * * * * * * * * //
+
+template<>
+tmp<surfaceScalarField> faceCorrectedSnGrad<scalar>::correction
+(
+    const volScalarField& vsf
+) const;
+
+
+template<>
+tmp<surfaceVectorField> faceCorrectedSnGrad<vector>::correction
+(
+    const volVectorField& vvf
+) const;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "faceCorrectedSnGrad.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/createBaffles/SetPatchFields.C b/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrads.C
similarity index 65%
rename from applications/utilities/mesh/manipulation/createBaffles/SetPatchFields.C
rename to src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrads.C
index df7c4dccb8..ba6d83542e 100644
--- a/applications/utilities/mesh/manipulation/createBaffles/SetPatchFields.C
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrads.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,24 +23,39 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "SetPatchFields.H"
+#include "faceCorrectedSnGrad.H"
+#include "fvMesh.H"
 
-// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template<class GeoField>
-void Foam::SetPatchFields
-(
-    PtrList<GeoField>& fields,
-    const label patchI,
-    const typename GeoField::value_type& initVal
-)
+namespace Foam
 {
-    forAll(fields, i)
-    {
-        typename GeoField::PatchFieldType& pfld =
-            fields[i].boundaryField()[patchI];
-        pfld == initVal;
-    }
+namespace fv
+{
+    makeSnGradScheme(faceCorrectedSnGrad)
+}
+}
+
+
+template<>
+Foam::tmp<Foam::surfaceScalarField>
+Foam::fv::faceCorrectedSnGrad<Foam::scalar>::correction
+(
+    const volScalarField& vsf
+) const
+{
+    return fullGradCorrection(vsf);
+}
+
+
+template<>
+Foam::tmp<Foam::surfaceVectorField>
+Foam::fv::faceCorrectedSnGrad<Foam::vector>::correction
+(
+    const volVectorField& vvf
+) const
+{
+    return fullGradCorrection(vvf);
 }
 
 
diff --git a/src/finiteVolume/fvMesh/fvMesh.H b/src/finiteVolume/fvMesh/fvMesh.H
index 97c1b0db58..7cb9d97a00 100644
--- a/src/finiteVolume/fvMesh/fvMesh.H
+++ b/src/finiteVolume/fvMesh/fvMesh.H
@@ -312,6 +312,9 @@ public:
             //- Return face centres as surfaceVectorField
             const surfaceVectorField& Cf() const;
 
+            //- Return face deltas as surfaceVectorField
+            tmp<surfaceVectorField> delta() const;
+
 
         // Edit
 
diff --git a/src/finiteVolume/fvMesh/fvMeshGeometry.C b/src/finiteVolume/fvMesh/fvMeshGeometry.C
index 5abdb45d11..ddd395a722 100644
--- a/src/finiteVolume/fvMesh/fvMeshGeometry.C
+++ b/src/finiteVolume/fvMesh/fvMeshGeometry.C
@@ -372,6 +372,53 @@ const surfaceVectorField& fvMesh::Cf() const
 }
 
 
+tmp<surfaceVectorField> fvMesh::delta() const
+{
+    if (debug)
+    {
+        Info<< "void fvMesh::delta() : "
+            << "calculating face deltas"
+            << endl;
+    }
+
+    tmp<surfaceVectorField> tdelta
+    (
+        new surfaceVectorField
+        (
+            IOobject
+            (
+                "delta",
+                pointsInstance(),
+                meshSubDir,
+                *this,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            *this,
+            dimLength
+        )
+    );
+    surfaceVectorField& delta = tdelta();
+
+    const volVectorField& C = this->C();
+    const labelUList& owner = this->owner();
+    const labelUList& neighbour = this->neighbour();
+
+    forAll(owner, facei)
+    {
+        delta[facei] = C[neighbour[facei]] - C[owner[facei]];
+    }
+
+    forAll(delta.boundaryField(), patchi)
+    {
+        delta.boundaryField()[patchi] = boundary()[patchi].delta();
+    }
+
+    return tdelta;
+}
+
+
 const surfaceScalarField& fvMesh::phi() const
 {
     if (!phiPtr_)
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/applications/utilities/mesh/manipulation/createBaffles/SetPatchFields.H b/src/meshTools/sets/faceZoneSources/searchableSurfaceToFaceZone/searchableSurfaceToFaceZone.H
similarity index 51%
rename from applications/utilities/mesh/manipulation/createBaffles/SetPatchFields.H
rename to src/meshTools/sets/faceZoneSources/searchableSurfaceToFaceZone/searchableSurfaceToFaceZone.H
index cf8a13d480..7269a3d38e 100644
--- a/applications/utilities/mesh/manipulation/createBaffles/SetPatchFields.H
+++ b/src/meshTools/sets/faceZoneSources/searchableSurfaceToFaceZone/searchableSurfaceToFaceZone.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) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -21,43 +21,84 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Global
-    Foam::SetPatchFields
+Class
+    Foam::searchableSurfaceToFaceZone
 
 Description
-    Helper routine to initialise a patch field to a constant value
+    A topoSetSource to select faces based on intersection (of cell-cell
+    vector) with a surface.
 
 SourceFiles
-    SetPatchFields.C
+    searchableSurfaceToFaceZone.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef SetPatchFields_H
-#define SetPatchFields_H
+#ifndef searchableSurfaceToFaceZone_H
+#define searchableSurfaceToFaceZone_H
 
-#include "PtrList.H"
+#include "topoSetSource.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-//- Helper routine to read fields
-template<class GeoField>
-void SetPatchFields
-(
-    PtrList<GeoField>& fields,
-    const label patchI,
-    const typename GeoField::value_type& initVal
-);
+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
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-#   include "SetPatchFields.C"
-#endif
+} // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/thermophysicalModels/basic/basicThermo/basicThermo.C b/src/thermophysicalModels/basic/basicThermo/basicThermo.C
index 84e30e31a4..e2fc874b5e 100644
--- a/src/thermophysicalModels/basic/basicThermo/basicThermo.C
+++ b/src/thermophysicalModels/basic/basicThermo/basicThermo.C
@@ -112,7 +112,7 @@ Foam::basicThermo::basicThermo
     (
         IOobject
         (
-            phasePropertyName("alpha"),
+            phasePropertyName("thermo:alpha"),
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ,
@@ -167,7 +167,7 @@ Foam::basicThermo::basicThermo
     (
         IOobject
         (
-            phasePropertyName("alpha"),
+            phasePropertyName("thermo:alpha"),
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ,
@@ -237,7 +237,7 @@ const Foam::basicThermo& Foam::basicThermo::lookupThermo
 
 void Foam::basicThermo::validate
 (
-    const word& app,
+    const string& app,
     const word& a
 ) const
 {
@@ -252,7 +252,7 @@ void Foam::basicThermo::validate
 
 void Foam::basicThermo::validate
 (
-    const word& app,
+    const string& app,
     const word& a,
     const word& b
 ) const
@@ -275,7 +275,7 @@ void Foam::basicThermo::validate
 
 void Foam::basicThermo::validate
 (
-    const word& app,
+    const string& app,
     const word& a,
     const word& b,
     const word& c
@@ -301,7 +301,7 @@ void Foam::basicThermo::validate
 
 void Foam::basicThermo::validate
 (
-    const word& app,
+    const string& app,
     const word& a,
     const word& b,
     const word& c,
diff --git a/src/thermophysicalModels/basic/basicThermo/basicThermo.H b/src/thermophysicalModels/basic/basicThermo/basicThermo.H
index 659b6f35ca..e652c30ee6 100644
--- a/src/thermophysicalModels/basic/basicThermo/basicThermo.H
+++ b/src/thermophysicalModels/basic/basicThermo/basicThermo.H
@@ -185,7 +185,7 @@ public:
         //  with energy forms supported by the application
         void validate
         (
-            const word& app,
+            const string& app,
             const word&
         ) const;
 
@@ -193,7 +193,7 @@ public:
         //  with energy forms supported by the application
         void validate
         (
-            const word& app,
+            const string& app,
             const word&,
             const word&
         ) const;
@@ -202,7 +202,7 @@ public:
         //  with energy forms supported by the application
         void validate
         (
-            const word& app,
+            const string& app,
             const word&,
             const word&,
             const word&
@@ -212,7 +212,7 @@ public:
         //  with energy forms supported by the application
         void validate
         (
-            const word& app,
+            const string& app,
             const word&,
             const word&,
             const word&,
@@ -263,6 +263,14 @@ public:
             //- Enthalpy/Internal energy [J/kg]
             virtual const volScalarField& he() const = 0;
 
+            //- Enthalpy/Internal energy
+            //  for given pressure and temperature [J/kg]
+            virtual tmp<volScalarField> he
+            (
+                const volScalarField& p,
+                const volScalarField& T
+            ) const = 0;
+
             //- Enthalpy/Internal energy for cell-set [J/kg]
             virtual tmp<scalarField> he
             (
diff --git a/src/thermophysicalModels/basic/heThermo/heThermo.C b/src/thermophysicalModels/basic/heThermo/heThermo.C
index 8efaf0f5d5..29634e2769 100644
--- a/src/thermophysicalModels/basic/heThermo/heThermo.C
+++ b/src/thermophysicalModels/basic/heThermo/heThermo.C
@@ -174,7 +174,10 @@ Foam::heThermo<BasicThermo, MixtureType>::heThermo
     (
         IOobject
         (
-            BasicThermo::phasePropertyName(MixtureType::thermoType::heName()),
+            BasicThermo::phasePropertyName
+            (
+                MixtureType::thermoType::heName()
+            ),
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ,
@@ -205,7 +208,10 @@ Foam::heThermo<BasicThermo, MixtureType>::heThermo
     (
         IOobject
         (
-            BasicThermo::phasePropertyName(MixtureType::thermoType::heName()),
+            BasicThermo::phasePropertyName
+            (
+                MixtureType::thermoType::heName()
+            ),
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ,
@@ -230,6 +236,60 @@ Foam::heThermo<BasicThermo, MixtureType>::~heThermo()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class BasicThermo, class MixtureType>
+Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::he
+(
+    const volScalarField& p,
+    const volScalarField& T
+) const
+{
+    const fvMesh& mesh = this->T_.mesh();
+
+    tmp<volScalarField> the
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "he",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            he_.dimensions()
+        )
+    );
+
+    volScalarField& he = the();
+    scalarField& heCells = he.internalField();
+    const scalarField& pCells = p.internalField();
+    const scalarField& TCells = T.internalField();
+
+    forAll(heCells, celli)
+    {
+        heCells[celli] =
+            this->cellMixture(celli).HE(pCells[celli], TCells[celli]);
+    }
+
+    forAll(he.boundaryField(), patchi)
+    {
+        scalarField& hep = he.boundaryField()[patchi];
+        const scalarField& pp = p.boundaryField()[patchi];
+        const scalarField& Tp = T.boundaryField()[patchi];
+
+        forAll(hep, facei)
+        {
+            hep[facei] =
+                this->patchFaceMixture(patchi, facei).HE(pp[facei], Tp[facei]);
+        }
+    }
+
+    return the;
+}
+
+
 template<class BasicThermo, class MixtureType>
 Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::he
 (
diff --git a/src/thermophysicalModels/basic/heThermo/heThermo.H b/src/thermophysicalModels/basic/heThermo/heThermo.H
index ef9919a99c..f77d017e13 100644
--- a/src/thermophysicalModels/basic/heThermo/heThermo.H
+++ b/src/thermophysicalModels/basic/heThermo/heThermo.H
@@ -161,6 +161,14 @@ public:
 
         // Fields derived from thermodynamic state variables
 
+            //- Enthalpy/Internal energy
+            //  for given pressure and temperature [J/kg]
+            virtual tmp<volScalarField> he
+            (
+                const volScalarField& p,
+                const volScalarField& T
+            ) const;
+
             //- Enthalpy/Internal energy for cell-set [J/kg]
             virtual tmp<scalarField> he
             (
diff --git a/src/thermophysicalModels/basic/psiThermo/psiThermo.C b/src/thermophysicalModels/basic/psiThermo/psiThermo.C
index fd4bc962c9..dd69890313 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermo.C
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermo.C
@@ -44,7 +44,7 @@ Foam::psiThermo::psiThermo(const fvMesh& mesh, const word& phaseName)
     (
         IOobject
         (
-            phasePropertyName("psi"),
+            phasePropertyName("thermo:psi"),
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ,
@@ -58,7 +58,7 @@ Foam::psiThermo::psiThermo(const fvMesh& mesh, const word& phaseName)
     (
         IOobject
         (
-            phasePropertyName("mu"),
+            phasePropertyName("thermo:mu"),
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ,
diff --git a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
index 4aa7835cac..d34186eea9 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
@@ -43,7 +43,7 @@ Foam::rhoThermo::rhoThermo(const fvMesh& mesh, const word& phaseName)
     (
         IOobject
         (
-            phasePropertyName("rhoThermo"),
+            phasePropertyName("thermo:rho"),
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ,
@@ -57,7 +57,7 @@ Foam::rhoThermo::rhoThermo(const fvMesh& mesh, const word& phaseName)
     (
         IOobject
         (
-            phasePropertyName("psi"),
+            phasePropertyName("thermo:psi"),
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ,
@@ -71,7 +71,7 @@ Foam::rhoThermo::rhoThermo(const fvMesh& mesh, const word& phaseName)
     (
         IOobject
         (
-            phasePropertyName("mu"),
+            phasePropertyName("thermo:mu"),
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ,
@@ -95,7 +95,7 @@ Foam::rhoThermo::rhoThermo
     (
         IOobject
         (
-            phasePropertyName("rhoThermo"),
+            phasePropertyName("thermo:rho"),
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ,
@@ -109,7 +109,7 @@ Foam::rhoThermo::rhoThermo
     (
         IOobject
         (
-            phasePropertyName("psi"),
+            phasePropertyName("thermo:psi"),
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ,
@@ -123,7 +123,7 @@ Foam::rhoThermo::rhoThermo
     (
         IOobject
         (
-            phasePropertyName("mu"),
+            phasePropertyName("thermo:mu"),
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ,
diff --git a/src/thermophysicalModels/basic/rhoThermo/rhoThermos.C b/src/thermophysicalModels/basic/rhoThermo/rhoThermos.C
index fcfc3437dd..beaa43b735 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermos.C
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermos.C
@@ -30,6 +30,7 @@ License
 #include "perfectGas.H"
 #include "incompressiblePerfectGas.H"
 #include "rhoConst.H"
+#include "perfectFluid.H"
 #include "hConstThermo.H"
 #include "janafThermo.H"
 #include "sensibleEnthalpy.H"
@@ -101,6 +102,18 @@ makeThermo
     specie
 );
 
+makeThermo
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleEnthalpy,
+    hConstThermo,
+    perfectFluid,
+    specie
+);
+
 makeThermo
 (
     rhoThermo,
@@ -200,6 +213,18 @@ makeThermo
     specie
 );
 
+makeThermo
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleInternalEnergy,
+    hConstThermo,
+    perfectFluid,
+    specie
+);
+
 makeThermo
 (
     rhoThermo,
diff --git a/src/thermophysicalModels/solidThermo/solidThermo/solidThermo.C b/src/thermophysicalModels/solidThermo/solidThermo/solidThermo.C
index 102414e3fb..2e12bcf213 100644
--- a/src/thermophysicalModels/solidThermo/solidThermo/solidThermo.C
+++ b/src/thermophysicalModels/solidThermo/solidThermo/solidThermo.C
@@ -50,7 +50,7 @@ Foam::solidThermo::solidThermo
     (
         IOobject
         (
-            "rhoThermo",
+            phasePropertyName("thermo:rho"),
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ,
@@ -74,7 +74,7 @@ Foam::solidThermo::solidThermo
     (
         IOobject
         (
-            "rhoThermo",
+            phasePropertyName("thermo:rho"),
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ,
diff --git a/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluid.C b/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluid.C
new file mode 100644
index 0000000000..844f293e2d
--- /dev/null
+++ b/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluid.C
@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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 "perfectFluid.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Specie>
+Foam::perfectFluid<Specie>::perfectFluid(Istream& is)
+:
+    Specie(is),
+    rho0_(readScalar(is))
+{
+    is.check("perfectFluid<Specie>::perfectFluid(Istream& is)");
+}
+
+
+template<class Specie>
+Foam::perfectFluid<Specie>::perfectFluid(const dictionary& dict)
+:
+    Specie(dict),
+    rho0_(readScalar(dict.subDict("equationOfState").lookup("rho0")))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Specie>
+void Foam::perfectFluid<Specie>::write(Ostream& os) const
+{
+    Specie::write(os);
+
+    dictionary dict("equationOfState");
+    dict.add("rho0", rho0_);
+
+    os  << indent << dict.dictName() << dict;
+}
+
+
+// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
+
+template<class Specie>
+Foam::Ostream& Foam::operator<<(Ostream& os, const perfectFluid<Specie>& pf)
+{
+    os  << static_cast<const Specie&>(pf)
+        << token::SPACE << pf.rho0_;
+
+    os.check("Ostream& operator<<(Ostream&, const perfectFluid<Specie>&)");
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluid.H b/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluid.H
new file mode 100644
index 0000000000..76cf80c4a8
--- /dev/null
+++ b/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluid.H
@@ -0,0 +1,223 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::perfectFluid
+
+Description
+    Perfect gas equation of state.
+
+SourceFiles
+    perfectFluidI.H
+    perfectFluid.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef perfectFluid_H
+#define perfectFluid_H
+
+#include "autoPtr.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+template<class Specie> class perfectFluid;
+
+template<class Specie>
+inline perfectFluid<Specie> operator+
+(
+    const perfectFluid<Specie>&,
+    const perfectFluid<Specie>&
+);
+
+template<class Specie>
+inline perfectFluid<Specie> operator-
+(
+    const perfectFluid<Specie>&,
+    const perfectFluid<Specie>&
+);
+
+template<class Specie>
+inline perfectFluid<Specie> operator*
+(
+    const scalar,
+    const perfectFluid<Specie>&
+);
+
+template<class Specie>
+inline perfectFluid<Specie> operator==
+(
+    const perfectFluid<Specie>&,
+    const perfectFluid<Specie>&
+);
+
+template<class Specie>
+Ostream& operator<<
+(
+    Ostream&,
+    const perfectFluid<Specie>&
+);
+
+
+/*---------------------------------------------------------------------------*\
+                           Class perfectFluid Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Specie>
+class perfectFluid
+:
+    public Specie
+{
+    // Private data
+
+        //- The reference density
+        scalar rho0_;
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        inline perfectFluid(const Specie& sp, const scalar rho0);
+
+        //- Construct from Istream
+        perfectFluid(Istream&);
+
+        //- Construct from dictionary
+        perfectFluid(const dictionary& dict);
+
+        //- Construct as named copy
+        inline perfectFluid(const word& name, const perfectFluid&);
+
+        //- Construct and return a clone
+        inline autoPtr<perfectFluid> clone() const;
+
+        // Selector from Istream
+        inline static autoPtr<perfectFluid> New(Istream& is);
+
+        // Selector from dictionary
+        inline static autoPtr<perfectFluid> New(const dictionary& dict);
+
+
+    // Member functions
+
+        //- Return the instantiated type name
+        static word typeName()
+        {
+            return "perfectFluid<" + word(Specie::typeName_()) + '>';
+        }
+
+
+        // Fundamental properties
+
+            //- Is the equation of state is incompressible i.e. rho != f(p)
+            static const bool incompressible = false;
+
+            //- Is the equation of state is isochoric i.e. rho = const
+            static const bool isochoric = false;
+
+            //- Return density [kg/m^3]
+            inline scalar rho(scalar p, scalar T) const;
+
+            //- Return compressibility rho/p [s^2/m^2]
+            inline scalar psi(scalar p, scalar T) const;
+
+            //- Return compression factor []
+            inline scalar Z(scalar p, scalar T) const;
+
+            //- Return (cp - cv) [J/(kmol K]
+            inline scalar cpMcv(scalar p, scalar T) const;
+
+
+        // IO
+
+            //- Write to Ostream
+            void write(Ostream& os) const;
+
+
+    // Member operators
+
+        inline void operator+=(const perfectFluid&);
+        inline void operator-=(const perfectFluid&);
+
+        inline void operator*=(const scalar);
+
+
+    // Friend operators
+
+        friend perfectFluid operator+ <Specie>
+        (
+            const perfectFluid&,
+            const perfectFluid&
+        );
+
+        friend perfectFluid operator- <Specie>
+        (
+            const perfectFluid&,
+            const perfectFluid&
+        );
+
+        friend perfectFluid operator* <Specie>
+        (
+            const scalar s,
+            const perfectFluid&
+        );
+
+        friend perfectFluid operator== <Specie>
+        (
+            const perfectFluid&,
+            const perfectFluid&
+        );
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<< <Specie>
+        (
+            Ostream&,
+            const perfectFluid&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "perfectFluidI.H"
+
+#ifdef NoRepository
+#   include "perfectFluid.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluidI.H b/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluidI.H
new file mode 100644
index 0000000000..21fc092fe2
--- /dev/null
+++ b/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluidI.H
@@ -0,0 +1,220 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "perfectFluid.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Specie>
+inline Foam::perfectFluid<Specie>::perfectFluid
+(
+    const Specie& sp,
+    const scalar rho0
+)
+:
+    Specie(sp),
+    rho0_(rho0)
+{}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Specie>
+inline Foam::perfectFluid<Specie>::perfectFluid
+(
+    const word& name,
+    const perfectFluid<Specie>& pf
+)
+:
+    Specie(name, pf),
+    rho0_(pf.rho0_)
+{}
+
+
+template<class Specie>
+inline Foam::autoPtr<Foam::perfectFluid<Specie> >
+Foam::perfectFluid<Specie>::clone() const
+{
+    return autoPtr<perfectFluid<Specie> >(new perfectFluid<Specie>(*this));
+}
+
+
+template<class Specie>
+inline Foam::autoPtr<Foam::perfectFluid<Specie> >
+Foam::perfectFluid<Specie>::New(Istream& is)
+{
+    return autoPtr<perfectFluid<Specie> >(new perfectFluid<Specie>(is));
+}
+
+
+template<class Specie>
+inline Foam::autoPtr<Foam::perfectFluid<Specie> >
+Foam::perfectFluid<Specie>::New
+(
+    const dictionary& dict
+)
+{
+    return autoPtr<perfectFluid<Specie> >(new perfectFluid<Specie>(dict));
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Specie>
+inline Foam::scalar Foam::perfectFluid<Specie>::rho(scalar p, scalar T) const
+{
+    return rho0_ + p/(this->R()*T);
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::perfectFluid<Specie>::psi(scalar, scalar T) const
+{
+    return 1.0/(this->R()*T);
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::perfectFluid<Specie>::Z(scalar, scalar) const
+{
+    return 1.0;
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::perfectFluid<Specie>::cpMcv(scalar, scalar) const
+{
+    return this->RR;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class Specie>
+inline void Foam::perfectFluid<Specie>::operator+=
+(
+    const perfectFluid<Specie>& pf
+)
+{
+    scalar molr1 = this->nMoles();
+
+    Specie::operator+=(pf);
+
+    molr1 /= this->nMoles();
+    scalar molr2 = pf.nMoles()/this->nMoles();
+
+    rho0_ = molr1*rho0_ + molr2*pf.rho0_;
+}
+
+
+template<class Specie>
+inline void Foam::perfectFluid<Specie>::operator-=
+(
+    const perfectFluid<Specie>& pf
+)
+{
+    scalar molr1 = this->nMoles();
+
+    Specie::operator-=(pf);
+
+    molr1 /= this->nMoles();
+    scalar molr2 = pf.nMoles()/this->nMoles();
+
+    rho0_ = molr1*rho0_ - molr2*pf.rho0_;
+}
+
+
+template<class Specie>
+inline void Foam::perfectFluid<Specie>::operator*=(const scalar s)
+{
+    Specie::operator*=(s);
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+template<class Specie>
+inline Foam::perfectFluid<Specie> Foam::operator+
+(
+    const perfectFluid<Specie>& pf1,
+    const perfectFluid<Specie>& pf2
+)
+{
+    scalar nMoles = pf1.nMoles() + pf2.nMoles();
+    scalar molr1 = pf1.nMoles()/nMoles;
+    scalar molr2 = pf2.nMoles()/nMoles;
+
+    return rhoConst<Specie>
+    (
+        static_cast<const Specie&>(pf1)
+      + static_cast<const Specie&>(pf2),
+        molr1*pf1.rho0_ + molr2*pf2.rho0_
+    );
+}
+
+
+template<class Specie>
+inline Foam::perfectFluid<Specie> Foam::operator-
+(
+    const perfectFluid<Specie>& pf1,
+    const perfectFluid<Specie>& pf2
+)
+{
+    scalar nMoles = pf1.nMoles() + pf2.nMoles();
+    scalar molr1 = pf1.nMoles()/nMoles;
+    scalar molr2 = pf2.nMoles()/nMoles;
+
+    return rhoConst<Specie>
+    (
+        static_cast<const Specie&>(pf1)
+      - static_cast<const Specie&>(pf2),
+        molr1*pf1.rho0_ - molr2*pf2.rho0_
+    );
+}
+
+
+template<class Specie>
+inline Foam::perfectFluid<Specie> Foam::operator*
+(
+    const scalar s,
+    const perfectFluid<Specie>& pf
+)
+{
+    return perfectFluid<Specie>(s*static_cast<const Specie&>(pf), pf.rho0_);
+}
+
+
+template<class Specie>
+inline Foam::perfectFluid<Specie> Foam::operator==
+(
+    const perfectFluid<Specie>& pf1,
+    const perfectFluid<Specie>& pf2
+)
+{
+    return pf2 - pf1;
+}
+
+
+// ************************************************************************* //
diff --git a/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
index ddde795545..321c07443b 100644
--- a/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
+++ b/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
@@ -97,6 +97,19 @@ tmp<volScalarField> kOmegaSST::F3() const
 }
 
 
+tmp<volScalarField> kOmegaSST::F23() const
+{
+    tmp<volScalarField> f23(F2());
+
+    if (F3_)
+    {
+        f23() *= F3();
+    }
+
+    return f23;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 kOmegaSST::kOmegaSST
@@ -228,6 +241,15 @@ kOmegaSST::kOmegaSST
             10.0
         )
     ),
+    F3_
+    (
+        Switch::lookupOrAddToDict
+        (
+            "F3",
+            coeffDict_,
+            false
+        )
+    ),
 
     y_(mesh_),
 
@@ -289,7 +311,7 @@ kOmegaSST::kOmegaSST
       / max
         (
             a1_*omega_,
-            b1_*F2()*F3()*sqrt(2.0)*mag(symm(fvc::grad(U_)))
+            b1_*F23()*sqrt(2.0)*mag(symm(fvc::grad(U_)))
         )
     );
     mut_.correctBoundaryConditions();
@@ -370,6 +392,7 @@ bool kOmegaSST::read()
         a1_.readIfPresent(coeffDict());
         b1_.readIfPresent(coeffDict());
         c1_.readIfPresent(coeffDict());
+        F3_.readIfPresent("F3", coeffDict());
 
         return true;
     }
@@ -470,7 +493,7 @@ void kOmegaSST::correct()
 
 
     // Re-calculate viscosity
-    mut_ = a1_*rho_*k_/max(a1_*omega_, b1_*F2()*F3()*sqrt(S2));
+    mut_ = a1_*rho_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2));
     mut_.correctBoundaryConditions();
 
     // Re-calculate thermal diffusivity
diff --git a/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.H b/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.H
index 6cc40f6b68..222242a74c 100644
--- a/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.H
+++ b/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.H
@@ -38,7 +38,7 @@ Description
         Nov. 2001
     \endverbatim
 
-    with the addition of the F3 term for rough walls from
+    with the addition of the optional F3 term for rough walls from
     \verbatim
         Hellsten, A.
         "Some Improvements in Menter’s k-omega-SST turbulence model"
@@ -80,6 +80,7 @@ Description
             a1          0.31;
             b1          1.0;
             c1          10.0;
+            F3          no;
         }
     \endverbatim
 
@@ -138,6 +139,8 @@ protected:
             dimensionedScalar b1_;
             dimensionedScalar c1_;
 
+            Switch F3_;
+
 
         //- Wall distance
         //  Note: different to wall distance in parent RASModel
@@ -156,6 +159,7 @@ protected:
         tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
         tmp<volScalarField> F2() const;
         tmp<volScalarField> F3() const;
+        tmp<volScalarField> F23() const;
 
         tmp<volScalarField> blend
         (
diff --git a/src/turbulenceModels/derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchField.C b/src/turbulenceModels/derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchField.C
index 07cc0c9481..afd43395ba 100644
--- a/src/turbulenceModels/derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchField.C
+++ b/src/turbulenceModels/derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -111,14 +111,10 @@ Foam::porousBafflePressureFvPatchField<Type>::porousBafflePressureFvPatchField
 template<class Type>
 void Foam::porousBafflePressureFvPatchField<Type>::write(Ostream& os) const
 {
-
     fixedJumpFvPatchField<Type>::write(os);
-
     os.writeKeyword("D") << D_ << token::END_STATEMENT << nl;
     os.writeKeyword("I") << I_ << token::END_STATEMENT << nl;
     os.writeKeyword("length") << length_ << token::END_STATEMENT << nl;
-
-    this->writeEntry("value", os);
 }
 
 
diff --git a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
index cff0d48146..290a4c73a8 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
+++ b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
@@ -98,6 +98,19 @@ tmp<volScalarField> kOmegaSST::F3() const
 }
 
 
+tmp<volScalarField> kOmegaSST::F23() const
+{
+    tmp<volScalarField> f23(F2());
+
+    if (F3_)
+    {
+        f23() *= F3();
+    }
+
+    return f23;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 kOmegaSST::kOmegaSST
@@ -219,6 +232,15 @@ kOmegaSST::kOmegaSST
             10.0
         )
     ),
+    F3_
+    (
+        Switch::lookupOrAddToDict
+        (
+            "F3",
+            coeffDict_,
+            false
+        )
+    ),
 
     y_(mesh_),
 
@@ -268,7 +290,7 @@ kOmegaSST::kOmegaSST
       / max
         (
             a1_*omega_,
-            b1_*F2()*F3()*sqrt(2.0)*mag(symm(fvc::grad(U_)))
+            b1_*F23()*sqrt(2.0)*mag(symm(fvc::grad(U_)))
         )
     );
     nut_.correctBoundaryConditions();
@@ -362,6 +384,7 @@ bool kOmegaSST::read()
         a1_.readIfPresent(coeffDict());
         b1_.readIfPresent(coeffDict());
         c1_.readIfPresent(coeffDict());
+        F3_.readIfPresent("F3", coeffDict());
 
         return true;
     }
@@ -439,7 +462,7 @@ void kOmegaSST::correct()
 
 
     // Re-calculate viscosity
-    nut_ = a1_*k_/max(a1_*omega_, b1_*F2()*F3()*sqrt(S2));
+    nut_ = a1_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2));
     nut_.correctBoundaryConditions();
 }
 
diff --git a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H
index 88feb733bc..07d16d2732 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H
+++ b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H
@@ -36,7 +36,7 @@ Description
         Nov. 2001.
     \endverbatim
 
-    with the addition of the F3 term for rough walls from
+    with the addition of the optional F3 term for rough walls from
     \verbatim
         Hellsten, A.
         "Some Improvements in Menter’s k-omega-SST turbulence model"
@@ -77,6 +77,7 @@ Description
             a1          0.31;
             b1          1.0;
             c1          10.0;
+            F3          no;
         }
     \endverbatim
 
@@ -132,6 +133,9 @@ protected:
             dimensionedScalar b1_;
             dimensionedScalar c1_;
 
+            Switch F3_;
+
+
         //- Wall distance field
         //  Note: different to wall distance in parent RASModel
         wallDist y_;
@@ -148,6 +152,7 @@ protected:
         tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
         tmp<volScalarField> F2() const;
         tmp<volScalarField> F3() const;
+        tmp<volScalarField> F23() const;
 
         tmp<volScalarField> blend
         (
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/T b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/T
index 8fbb73eff0..92fbe4aaac 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/T
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/T
@@ -44,10 +44,6 @@ boundaryField
     {
         type            empty;
     }
-    "baffle1Wall.*"
-    {
-        type           calculated;
-    }
 }
 
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/U b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/U
index 1cee83b6cf..76239f12dd 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/U
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/U
@@ -46,10 +46,6 @@ boundaryField
     {
         type            empty;
     }
-    "baffle1Wall.*"
-    {
-        type            calculated;
-    }
 }
 
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/alphat b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/alphat
index 977203512a..3bc99dd808 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/alphat
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/alphat
@@ -45,10 +45,6 @@ boundaryField
     {
         type            empty;
     }
-    "baffle1Wall.*"
-    {
-        type            calculated;
-    }
 }
 
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/epsilon b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/epsilon
index 21a6d42fa8..69f9ac6803 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/epsilon
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/epsilon
@@ -44,10 +44,6 @@ boundaryField
     {
         type            empty;
     }
-    "baffle1Wall.*"
-    {
-         type           calculated;
-    }
 }
 
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/k b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/k
index 5f46125669..9ecfe3d9f8 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/k
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/k
@@ -44,10 +44,6 @@ boundaryField
     {
         type            empty;
     }
-    "baffle1Wall.*"
-    {
-        type            calculated;
-    }
 }
 
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/mut b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/mut
index 41310b7ced..f83a004214 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/mut
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/mut
@@ -45,10 +45,6 @@ boundaryField
     {
         type            empty;
     }
-    "baffle1Wall.*"
-    {
-        type           calculated;
-    }
 }
 
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/p b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/p
index 8bbeb86853..d20b276591 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/p
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/p
@@ -45,10 +45,6 @@ boundaryField
     {
         type            empty;
     }
-    "baffle1Wall.*"
-    {
-        type           calculated;
-    }
 }
 
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/p_rgh b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/p_rgh
index 8cb7d58cdf..26fa39b255 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/p_rgh
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/p_rgh
@@ -45,10 +45,6 @@ boundaryField
     {
         type            empty;
     }
-    "baffle1Wall.*"
-    {
-        type           calculated;
-    }
 }
 
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/Allrun b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/Allrun
index ed131900d3..e83f5f5977 100755
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/Allrun
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/Allrun
@@ -13,15 +13,11 @@ unset FOAM_SETNAN
 unset FOAM_SIGFPE
 
 # Create first baffle
-runApplication createBaffles baffleFaces '(baffle1Wall_0 baffle1Wall_1)' -overwrite
+runApplication createBaffles -overwrite
 
 # Create region
 runApplication extrudeToRegionMesh -overwrite
 
-# Set the BC's for the baffle
-runApplication changeDictionary -dict system/changeDictionaryDict.baffle
-rm log.changeDictionary
-
 # Set Bc's for the region baffle
 runApplication changeDictionary -dict system/changeDictionaryDict.baffleRegion -literalRE
 rm log.changeDictionary
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/polyMesh/blockMeshDict b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/polyMesh/blockMeshDict
index b1def0a3a2..4b273990da 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/polyMesh/blockMeshDict
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/polyMesh/blockMeshDict
@@ -84,28 +84,6 @@ boundary
             (4 5 6 7)
         );
     }
-
-    baffle1Wall_0
-    {
-        type            mappedWall;
-        sampleMode      nearestPatchFace;
-        sampleRegion    region0;
-        samplePatch     baffle1Wall_1;
-        offsetMode      uniform;
-        offset          (0 0 0);
-        faces ();
-    }
-
-    baffle1Wall_1
-    {
-        type            mappedWall;
-        sampleMode      nearestPatchFace;
-        sampleRegion    region0;
-        samplePatch     baffle1Wall_0;
-        offsetMode      uniform;
-        offset          (0 0 0);
-        faces ();
-    }
 );
 
 mergePatchPairs
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/controlDict b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/controlDict
index ace5803d49..acd337d877 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/controlDict
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/controlDict
@@ -44,6 +44,6 @@ timePrecision   6;
 
 runTimeModifiable true;
 
-libs ("libthermoBaffleModels.so");
+libs ("libthermoBaffleModels.so" "libcompressibleRASModels.so");
 
 // ************************************************************************* //
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;
diff --git a/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/U b/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/U
index d1bbe70951..c78eec59c1 100644
--- a/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/U
+++ b/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/U
@@ -45,18 +45,6 @@ boundaryField
         type            movingWallVelocity;
         value           uniform (0 0 0);
     }
-
-    AMI1
-    {
-        type            cyclicAMI;
-        value           uniform (0 0 0);
-    }
-
-    AMI2
-    {
-        type            cyclicAMI;
-        value           uniform (0 0 0);
-    }
 }
 
 
diff --git a/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/epsilon b/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/epsilon
index 4512a15707..c8a9ab3e88 100644
--- a/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/epsilon
+++ b/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/epsilon
@@ -45,18 +45,6 @@ boundaryField
         type            epsilonWallFunction;
         value           $internalField;
     }
-
-    AMI1
-    {
-        type            cyclicAMI;
-        value           $internalField;
-    }
-
-    AMI2
-    {
-        type            cyclicAMI;
-        value           $internalField;
-    }
 }
 
 
diff --git a/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/k b/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/k
index 4539cc1591..fa184b506d 100644
--- a/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/k
+++ b/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/k
@@ -45,18 +45,6 @@ boundaryField
         type            kqRWallFunction;
         value           $internalField;
     }
-
-    AMI1
-    {
-        type            cyclicAMI;
-        value           $internalField;
-    }
-
-    AMI2
-    {
-        type            cyclicAMI;
-        value           $internalField;
-    }
 }
 
 
diff --git a/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/nut b/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/nut
index 3b34b2cbbf..0c69619423 100644
--- a/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/nut
+++ b/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/nut
@@ -44,18 +44,6 @@ boundaryField
         type            nutkWallFunction;
         value           uniform 0;
     }
-
-    AMI1
-    {
-        type            cyclicAMI;
-        value           uniform 0;
-    }
-
-    AMI2
-    {
-        type            cyclicAMI;
-        value           uniform 0;
-    }
 }
 
 
diff --git a/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/p b/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/p
index 3196c16c90..98504cb29b 100644
--- a/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/p
+++ b/tutorials/incompressible/pimpleDyMFoam/propeller/0.org/p
@@ -41,18 +41,6 @@ boundaryField
     {
         type            zeroGradient;
     }
-
-    AMI1
-    {
-        type            cyclicAMI;
-        value           uniform 0;
-    }
-
-    AMI2
-    {
-        type            cyclicAMI;
-        value           uniform 0;
-    }
 }
 
 
diff --git a/tutorials/incompressible/pimpleDyMFoam/propeller/Allrun.pre b/tutorials/incompressible/pimpleDyMFoam/propeller/Allrun.pre
index 472c03a6f7..2abd5b41e1 100755
--- a/tutorials/incompressible/pimpleDyMFoam/propeller/Allrun.pre
+++ b/tutorials/incompressible/pimpleDyMFoam/propeller/Allrun.pre
@@ -17,6 +17,7 @@ runApplication surfaceFeatureExtract
 runApplication snappyHexMesh -overwrite
 # force removal of fields generated by snappy
 \rm -rf 0
+cp -rf 0.org 0
 
 
 # - generate face/cell sets and zones
@@ -37,19 +38,18 @@ runApplication topoSet -dict system/createAMIFaces.topoSetDict
 mv log.topoSet log.createAMIFaces.topoSet
 
 
-# - create the inlet/outlet patches and AMI patches
+# - create the inlet/outlet patches
 
 runApplication createPatch -overwrite
 
 
-# - create the AMI faces by creating baffles, and then splitting the mesh
-
-createBaffles -internalFacesOnly -overwrite innerCylinderSmall '(AMI1 AMI2)'   \
-    > log.createBaffles 2>&1
-
-runApplication mergeOrSplitBaffles -split -overwrite
-
-
 # - apply the initial fields
 
 cp -rf 0.org 0
+
+
+# - create the AMI faces by creating baffles, and then splitting the mesh
+
+runApplication createBaffles -overwrite
+
+runApplication mergeOrSplitBaffles -split -overwrite
diff --git a/tutorials/incompressible/pimpleDyMFoam/propeller/system/controlDict b/tutorials/incompressible/pimpleDyMFoam/propeller/system/controlDict
index babfb5351e..96dc9e9cfd 100644
--- a/tutorials/incompressible/pimpleDyMFoam/propeller/system/controlDict
+++ b/tutorials/incompressible/pimpleDyMFoam/propeller/system/controlDict
@@ -15,6 +15,12 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+libs
+(
+    "libincompressibleTurbulenceModel.so"
+    "libincompressibleRASModels.so"
+);
+
 application     pimpleDyMFoam;
 
 startFrom       startTime;
diff --git a/tutorials/incompressible/pimpleDyMFoam/propeller/system/createBafflesDict b/tutorials/incompressible/pimpleDyMFoam/propeller/system/createBafflesDict
new file mode 100644
index 0000000000..f3a7eadc97
--- /dev/null
+++ b/tutorials/incompressible/pimpleDyMFoam/propeller/system/createBafflesDict
@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      createBafflesDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Whether to convert internal faces only (so leave boundary faces intact).
+// This is only relevant if your face selection type can pick up boundary
+// faces.
+internalFacesOnly true;
+
+
+// Baffles to create.
+baffles
+{
+    baffleFaces
+    {
+        //- Use predefined faceZone to select faces and orientation.
+        type        faceZone;
+        zoneName    innerCylinderSmall;
+
+
+        patches
+        {
+            master
+            {
+                //- Master side patch
+                name            AMI1;
+                type            cyclicAMI;
+                matchTolerance  0.0001;
+                neighbourPatch  AMI2;
+                transform       noOrdering;
+            }
+            slave
+            {
+                //- Slave side patch
+                name            AMI2;
+                type            cyclicAMI;
+                matchTolerance  0.0001;
+                neighbourPatch  AMI1;
+                transform       noOrdering;
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleDyMFoam/propeller/system/createPatchDict b/tutorials/incompressible/pimpleDyMFoam/propeller/system/createPatchDict
index 2400a43509..1559be08dd 100644
--- a/tutorials/incompressible/pimpleDyMFoam/propeller/system/createPatchDict
+++ b/tutorials/incompressible/pimpleDyMFoam/propeller/system/createPatchDict
@@ -46,37 +46,6 @@ patches
 
         set outletFaces;
     }
-
-
-    // Construct empty patches
-    {
-        name AMI1;
-        patchInfo
-        {
-            type            cyclicAMI;
-            neighbourPatch  AMI2;
-            transform       noOrdering;
-            surface
-            {}
-        }
-        constructFrom set;
-        set dummyFaces;
-    }
-    {
-        name AMI2;
-
-        patchInfo
-        {
-            type            cyclicAMI;
-            neighbourPatch  AMI1;
-            transform       noOrdering;
-            surface
-            {}
-        }
-        constructFrom set;
-        set dummyFaces;
-    }
-
 );
 
 // ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/TJunctionFan/0/U b/tutorials/incompressible/pimpleFoam/TJunctionFan/0/U
index 7d734082dd..54e519308c 100644
--- a/tutorials/incompressible/pimpleFoam/TJunctionFan/0/U
+++ b/tutorials/incompressible/pimpleFoam/TJunctionFan/0/U
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [0 1 -1 0 0 0 0];
+dimensions      [m s^-1];
 
 internalField   uniform (0 0 0);
 
@@ -38,19 +38,6 @@ boundaryField
         inletValue      uniform (0 0 0);
         value           uniform (0 0 0);
     }
-    baffles
-    {
-        type            fixedValue;
-        value           uniform (0 0 0);
-    }
-    fan_half0
-    {
-        type            cyclic;
-    }
-    fan_half1
-    {
-        type            cyclic;
-    }
     defaultFaces
     {
         type            fixedValue;
diff --git a/tutorials/incompressible/pimpleFoam/TJunctionFan/0/epsilon b/tutorials/incompressible/pimpleFoam/TJunctionFan/0/epsilon
index b972d9f1d2..3648316bd3 100644
--- a/tutorials/incompressible/pimpleFoam/TJunctionFan/0/epsilon
+++ b/tutorials/incompressible/pimpleFoam/TJunctionFan/0/epsilon
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [0 2 -3 0 0 0 0];
+dimensions      [m^2 s^-3];
 
 internalField   uniform 200;
 
@@ -41,22 +41,6 @@ boundaryField
         inletValue      uniform 200;
         value           uniform 200;
     }
-    baffles
-    {
-        type            epsilonWallFunction;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
-        value           uniform 0;
-    }
-    fan_half0
-    {
-        type            cyclic;
-    }
-    fan_half1
-    {
-        type            cyclic;
-    }
     defaultFaces
     {
         type            epsilonWallFunction;
diff --git a/tutorials/incompressible/pimpleFoam/TJunctionFan/0/k b/tutorials/incompressible/pimpleFoam/TJunctionFan/0/k
index a58acefeab..bd8e325c2f 100644
--- a/tutorials/incompressible/pimpleFoam/TJunctionFan/0/k
+++ b/tutorials/incompressible/pimpleFoam/TJunctionFan/0/k
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [0 2 -2 0 0 0 0];
+dimensions      [m^2 s^-2];
 
 internalField   uniform 0.375;
 
@@ -25,8 +25,6 @@ boundaryField
     {
         type            turbulentIntensityKineticEnergyInlet;
         intensity       0.05;
-        U               U;
-        phi             phi;
         value           uniform 0.375;
     }
     outlet1
@@ -41,19 +39,6 @@ boundaryField
         inletValue      uniform 0.375;
         value           uniform 0.375;
     }
-    baffles
-    {
-        type            kqRWallFunction;
-        value           uniform 0;
-    }
-    fan_half0
-    {
-        type            cyclic;
-    }
-    fan_half1
-    {
-        type            cyclic;
-    }
     defaultFaces
     {
         type            kqRWallFunction;
diff --git a/tutorials/incompressible/pimpleFoam/TJunctionFan/0/nuTilda b/tutorials/incompressible/pimpleFoam/TJunctionFan/0/nuTilda
index dfa76ef6fd..5bd1e94d8b 100644
--- a/tutorials/incompressible/pimpleFoam/TJunctionFan/0/nuTilda
+++ b/tutorials/incompressible/pimpleFoam/TJunctionFan/0/nuTilda
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [0 2 -1 0 0 0 0];
+dimensions      [m^2 s^-1];
 
 internalField   uniform 0;
 
@@ -33,18 +33,6 @@ boundaryField
     {
         type            zeroGradient;
     }
-    baffles
-    {
-        type            zeroGradient;
-    }
-    fan_half0
-    {
-        type            cyclic;
-    }
-    fan_half1
-    {
-        type            cyclic;
-    }
     defaultFaces
     {
         type            zeroGradient;
diff --git a/tutorials/incompressible/pimpleFoam/TJunctionFan/0/nut b/tutorials/incompressible/pimpleFoam/TJunctionFan/0/nut
index 831ea1e70b..894d5b5345 100644
--- a/tutorials/incompressible/pimpleFoam/TJunctionFan/0/nut
+++ b/tutorials/incompressible/pimpleFoam/TJunctionFan/0/nut
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [0 2 -1 0 0 0 0];
+dimensions      [m^2 s^-1];
 
 internalField   uniform 0;
 
@@ -36,22 +36,6 @@ boundaryField
         type            calculated;
         value           uniform 0;
     }
-    baffles
-    {
-        type            nutkWallFunction;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
-        value           uniform 0;
-    }
-    fan_half0
-    {
-        type            cyclic;
-    }
-    fan_half1
-    {
-        type            cyclic;
-    }
     defaultFaces
     {
         type            nutkWallFunction;
diff --git a/tutorials/incompressible/pimpleFoam/TJunctionFan/0/p b/tutorials/incompressible/pimpleFoam/TJunctionFan/0/p
index e69dbaddef..42728c992c 100644
--- a/tutorials/incompressible/pimpleFoam/TJunctionFan/0/p
+++ b/tutorials/incompressible/pimpleFoam/TJunctionFan/0/p
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [0 2 -2 0 0 0 0];
+dimensions      [m^2 s^-2];
 
 internalField   uniform 0;
 
@@ -27,15 +27,13 @@ boundaryField
         rho             none;
         psi             none;
         gamma           1;
-        p0              40;
         pressure        table
-
-2
-(
-(0 10)
-(1 40)
-)
-;
+        2
+        (
+            (0 10)
+            (1 40)
+        )
+        ;
         value           uniform 40;
     }
     outlet1
@@ -48,30 +46,6 @@ boundaryField
         type            fixedValue;
         value           uniform 0;
     }
-    baffles
-    {
-        type            zeroGradient;
-    }
-    fan_half0
-    {
-        type            fan;
-        patchType       cyclic;
-        jumpTable       polynomial
-
-2
-(
-(100 0)
-(-0.1 1)
-)
-;
-        value           uniform 0;
-    }
-    fan_half1
-    {
-        type            fan;
-        patchType       cyclic;
-        value           uniform 0;
-    }
     defaultFaces
     {
         type            zeroGradient;
diff --git a/tutorials/incompressible/pimpleFoam/TJunctionFan/Allrun b/tutorials/incompressible/pimpleFoam/TJunctionFan/Allrun
index 2374d2fb6c..932dded90b 100755
--- a/tutorials/incompressible/pimpleFoam/TJunctionFan/Allrun
+++ b/tutorials/incompressible/pimpleFoam/TJunctionFan/Allrun
@@ -8,21 +8,10 @@ application=`getApplication`
 
 runApplication blockMesh
 
-
-# Creating patches 'out of nothing' so will get undefined value which
-# might cause floating point error
-unset FOAM_SIGFPE
-unset FOAM_SETNAN
-
-# Create faceZones for fan and baffles
+# Create faceZones for baffles and fan
 runApplication topoSet
 
-# Create fan cyclics
-runApplication createBaffles cyclicFaces '(fan_half0 fan_half1)' -overwrite
-mv log.createBaffles log.createBaffles.1
-
-# Create wall baffles
-runApplication createBaffles baffleFaces '(baffles baffles)' -overwrite
-mv log.createBaffles log.createBaffles.2
+# Create wall and cyclic baffles and the fields on them
+runApplication createBaffles -overwrite
 
 runApplication $application
diff --git a/tutorials/incompressible/pimpleFoam/TJunctionFan/README.txt b/tutorials/incompressible/pimpleFoam/TJunctionFan/README.txt
index 0ef40d8110..7c7924cf56 100644
--- a/tutorials/incompressible/pimpleFoam/TJunctionFan/README.txt
+++ b/tutorials/incompressible/pimpleFoam/TJunctionFan/README.txt
@@ -1,3 +1,3 @@
 2010-04-16 T-junction with baffles and fan.
 Inlet on left, one outlet at bottom, one at top.
-To test multiple outlets.
+To test multiple outlets and fan cyclics
diff --git a/tutorials/incompressible/pimpleFoam/TJunctionFan/constant/polyMesh/blockMeshDict b/tutorials/incompressible/pimpleFoam/TJunctionFan/constant/polyMesh/blockMeshDict
index 66f0af7d65..7464688c33 100644
--- a/tutorials/incompressible/pimpleFoam/TJunctionFan/constant/polyMesh/blockMeshDict
+++ b/tutorials/incompressible/pimpleFoam/TJunctionFan/constant/polyMesh/blockMeshDict
@@ -104,25 +104,6 @@ boundary
         faces ((8 18 19 9));
     }
 
-    baffles
-    {
-        type wall;
-        faces ();
-    }
-
-    fan_half0
-    {
-        type cyclic;
-        faces ();
-        neighbourPatch fan_half1;
-    }
-    fan_half1
-    {
-        type cyclic;
-        faces ();
-        neighbourPatch fan_half0;
-    }
-
     defaultFaces
     {
         type wall;
diff --git a/tutorials/incompressible/pimpleFoam/TJunctionFan/system/createBafflesDict b/tutorials/incompressible/pimpleFoam/TJunctionFan/system/createBafflesDict
new file mode 100644
index 0000000000..92c3f3dd82
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/TJunctionFan/system/createBafflesDict
@@ -0,0 +1,186 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      createBafflesDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Whether to convert internal faces only (so leave boundary faces intact).
+// This is only relevant if your face selection type can pick up boundary
+// faces.
+internalFacesOnly true;
+
+
+// Baffles to create.
+baffles
+{
+    baffleFaces
+    {
+        //- Use predefined faceZone to select faces and orientation.
+        type        faceZone;
+        zoneName    baffleFaces;
+
+
+        //- Optional flip
+        //flip        false;
+
+        patches
+        {
+            master
+            {
+                //- Master side patch
+                name            baffles;
+                type            wall;
+
+                patchFields
+                {
+                    epsilon
+                    {
+                        type            epsilonWallFunction;
+                        Cmu             0.09;
+                        kappa           0.41;
+                        E               9.8;
+                        value           uniform 0;
+                    }
+                    k
+                    {
+                        type            kqRWallFunction;
+                        value           uniform 0;
+                    }
+                    nut
+                    {
+                        type            nutkWallFunction;
+                        Cmu             0.09;
+                        kappa           0.41;
+                        E               9.8;
+                        value           uniform 0;
+                    }
+                    nuTilda
+                    {
+                        type            zeroGradient;
+                    }
+                    p
+                    {
+                        type            zeroGradient;
+                    }
+                    U
+                    {
+                        type            fixedValue;
+                        value           uniform (0 0 0);
+                    }
+                }
+            }
+            slave
+            {
+                //- Slave side patch
+                name            baffles;
+                type            wall;
+
+                patchFields
+                {
+                    epsilon
+                    {
+                        type            epsilonWallFunction;
+                        Cmu             0.09;
+                        kappa           0.41;
+                        E               9.8;
+                        value           uniform 0;
+                    }
+                    k
+                    {
+                        type            kqRWallFunction;
+                        value           uniform 0;
+                    }
+                    nut
+                    {
+                        type            nutkWallFunction;
+                        Cmu             0.09;
+                        kappa           0.41;
+                        E               9.8;
+                        value           uniform 0;
+                    }
+                    nuTilda
+                    {
+                        type            zeroGradient;
+                    }
+                    p
+                    {
+                        type            zeroGradient;
+                    }
+                    U
+                    {
+                        type            fixedValue;
+                        value           uniform (0 0 0);
+                    }
+                }
+            }
+        }
+    }
+
+
+
+    cyclicFaces
+    {
+        //- Select faces and orientation through a searchableSurface
+        type        searchableSurface;
+        surface     searchablePlate;
+        origin      (0.099 -0.006 0.004);
+        span        (0 0.012 0.012);
+
+        patches
+        {
+            master
+            {
+                //- Master side patch
+
+                name            fan_half0;
+                type            cyclic;
+                neighbourPatch  fan_half1;
+
+                //- Optional override of added patchfields. If not specified
+                //  any added patchfields are of type calculated.
+                patchFields
+                {
+                    p
+                    {
+                        type            fan;
+                        patchType       cyclic;
+                        jump            uniform 0;
+                        value           uniform 0;
+                        jumpTable       polynomial 1((100 0));
+                    }
+                }
+            }
+            slave
+            {
+                //- Slave side patch
+
+                name            fan_half1;
+                type            cyclic;
+                neighbourPatch  fan_half0;
+
+                patchFields
+                {
+                    p
+                    {
+                        type            fan;
+                        patchType       cyclic;
+                        value           uniform 0;
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict b/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
index 6b36b1ef39..ebed3bc292 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
@@ -193,17 +193,17 @@ snapControls
 
         //- Number of feature edge snapping iterations.
         //  Leave out altogether to disable.
-        nFeatureSnapIter 10;
+        nFeatureSnapIter 0;
 
         //- Detect (geometric only) features by sampling the surface
         //  (default=false).
         implicitFeatureSnap false;
 
         //- Use castellatedMeshControls::features (default = true)
-        explicitFeatureSnap true;
+        explicitFeatureSnap false;
 
         //- Detect points on multiple surfaces (only for explicitFeatureSnap)
-        multiRegionFeatureSnap true;
+        multiRegionFeatureSnap false;
 }
 
 
diff --git a/tutorials/lagrangian/reactingParcelFoam/filter/Allrun b/tutorials/lagrangian/reactingParcelFoam/filter/Allrun
index b82b7aad36..7db81c2676 100755
--- a/tutorials/lagrangian/reactingParcelFoam/filter/Allrun
+++ b/tutorials/lagrangian/reactingParcelFoam/filter/Allrun
@@ -10,23 +10,10 @@ application=`getApplication`
 runApplication blockMesh
 
 # create sets
-#setSet -batch system/sets.setSet > log.setSet1 2>&1
 runApplication topoSet
 
-# create the first cyclic - lhs of porous zone
-# Note that we don't know what value to give these patches-out-of-nothing so
-# - use binary writing to avoid 'nan'
-# - use setFields to set values
-unset FOAM_SIGFPE
-runApplication createBaffles cycLeft '(cycLeft_half0 cycLeft_half1)' -overwrite
-mv log.createBaffles log.createBaffles1
-
-# create the second cyclic - rhs of porous zone
-runApplication createBaffles cycRight '(cycRight_half0 cycRight_half1)' -overwrite
-mv log.createBaffles log.createBaffles2
-
-# Initialise newly created patchFields to 0
-runApplication changeDictionary
+# create baffles and fields
+createBaffles -overwrite
 
 runApplication $application
 
diff --git a/tutorials/lagrangian/reactingParcelFoam/filter/constant/polyMesh/blockMeshDict b/tutorials/lagrangian/reactingParcelFoam/filter/constant/polyMesh/blockMeshDict
index 04b1edd6ab..ae16964f93 100644
--- a/tutorials/lagrangian/reactingParcelFoam/filter/constant/polyMesh/blockMeshDict
+++ b/tutorials/lagrangian/reactingParcelFoam/filter/constant/polyMesh/blockMeshDict
@@ -100,32 +100,6 @@ boundary
         );
     }
 
-    cycLeft_half0
-    {
-        type cyclic;
-        faces    ();
-        neighbourPatch cycLeft_half1;
-    }
-    cycLeft_half1
-    {
-        type cyclic;
-        faces    ();
-        neighbourPatch cycLeft_half0;
-    }
-
-    cycRight_half0
-    {
-        type cyclic;
-        faces    ();
-        neighbourPatch cycRight_half1;
-    }
-    cycRight_half1
-    {
-        type cyclic;
-        faces    ();
-        neighbourPatch cycRight_half0;
-    }
-
     frontAndBack
     {
         type empty;
diff --git a/tutorials/lagrangian/reactingParcelFoam/filter/constant/polyMesh/boundary b/tutorials/lagrangian/reactingParcelFoam/filter/constant/polyMesh/boundary
deleted file mode 100644
index 3a75efa15b..0000000000
--- a/tutorials/lagrangian/reactingParcelFoam/filter/constant/polyMesh/boundary
+++ /dev/null
@@ -1,83 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      binary;
-    class       polyBoundaryMesh;
-    location    "constant/polyMesh";
-    object      boundary;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-8
-(
-    walls
-    {
-        type            wall;
-        nFaces          172;
-        startFace       3294;
-    }
-    inlet
-    {
-        type            patch;
-        nFaces          20;
-        startFace       3466;
-    }
-    outlet
-    {
-        type            patch;
-        nFaces          20;
-        startFace       3486;
-    }
-    cycLeft_half0
-    {
-        type            cyclic;
-        inGroups        1(cyclic);
-        nFaces          20;
-        startFace       3506;
-        matchTolerance  0.0001;
-        neighbourPatch  cycLeft_half1;
-    }
-    cycLeft_half1
-    {
-        type            cyclic;
-        inGroups        1(cyclic);
-        nFaces          20;
-        startFace       3526;
-        matchTolerance  0.0001;
-        neighbourPatch  cycLeft_half0;
-    }
-    cycRight_half0
-    {
-        type            cyclic;
-        inGroups        1(cyclic);
-        nFaces          20;
-        startFace       3546;
-        matchTolerance  0.0001;
-        neighbourPatch  cycRight_half1;
-    }
-    cycRight_half1
-    {
-        type            cyclic;
-        inGroups        1(cyclic);
-        nFaces          20;
-        startFace       3566;
-        matchTolerance  0.0001;
-        neighbourPatch  cycRight_half0;
-    }
-    frontAndBack
-    {
-        type            empty;
-        inGroups        1(empty);
-        nFaces          3440;
-        startFace       3586;
-    }
-)
-
-// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/filter/system/createBafflesDict b/tutorials/lagrangian/reactingParcelFoam/filter/system/createBafflesDict
new file mode 100644
index 0000000000..1caf47dc89
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/filter/system/createBafflesDict
@@ -0,0 +1,81 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      createBafflesDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Whether to convert internal faces only (so leave boundary faces intact).
+// This is only relevant if your face selection type can pick up boundary
+// faces.
+internalFacesOnly true;
+
+// Baffles to create.
+baffles
+{
+    cycLeft
+    {
+        //- Use predefined faceZone to select faces and orientation.
+        type        faceZone;
+        zoneName    cycLeft;
+
+        patches
+        {
+            master
+            {
+                //- Master side patch
+
+                name            cycLeft_half0;
+                type            cyclic;
+                neighbourPatch  cycLeft_half1;
+            }
+            slave
+            {
+                //- Slave side patch
+
+                name            cycLeft_half1;
+                type            cyclic;
+                neighbourPatch  cycLeft_half0;
+            }
+        }
+    }
+
+    cycRight
+    {
+        //- Use predefined faceZone to select faces and orientation.
+        type        faceZone;
+        zoneName    cycRight;
+
+        patches
+        {
+            master
+            {
+                //- Master side patch
+
+                name            cycRight_half0;
+                type            cyclic;
+                neighbourPatch  cycRight_half1;
+            }
+            slave
+            {
+                //- Slave side patch
+
+                name            cycRight_half1;
+                type            cyclic;
+                neighbourPatch  cycRight_half0;
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/0/alpha1 b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/0/alpha1
index 6e0ae20fd8..c1fa13b514 100644
--- a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/0/alpha1
+++ b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/0/alpha1
@@ -39,14 +39,6 @@ boundaryField
         inletValue      uniform 0;
         value           uniform 0;
     }
-    porous_half0
-    {
-        type            cyclic;
-    }
-    porous_half1
-    {
-        type            cyclic;
-    }
     defaultFaces
     {
         type            empty;
diff --git a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/Allrun b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/Allrun
index 32c7e6f948..2b6ce02153 100755
--- a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/Allrun
+++ b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/Allrun
@@ -10,15 +10,7 @@ application=`getApplication`
 runApplication blockMesh
 runApplication setFields
 
-unset FOAM_SIGFPE
-unset FOAM_SETNAN
-
-# Create faceZones for porous baffles
-runApplication topoSet
-
-runApplication createBaffles cyclicZoneFaces '(porous_half0 porous_half1)' -overwrite
-
-runApplication changeDictionary
+runApplication createBaffles -overwrite
 
 runApplication $application
 
diff --git a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/constant/polyMesh/blockMeshDict b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/constant/polyMesh/blockMeshDict
index a60b3257e3..11344be6ac 100644
--- a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/constant/polyMesh/blockMeshDict
+++ b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/constant/polyMesh/blockMeshDict
@@ -99,20 +99,6 @@ boundary
             (10 22 23 11)
         );
     }
-
-    porous_half0
-    {
-        type cyclic;
-        faces ();
-        neighbourPatch porous_half1;
-    }
-
-    porous_half1
-    {
-        type cyclic;
-        faces ();
-        neighbourPatch porous_half0;
-    }
 );
 
 mergePatchPairs
diff --git a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/constant/polyMesh/boundary b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/constant/polyMesh/boundary
deleted file mode 100644
index de4cd85200..0000000000
--- a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/constant/polyMesh/boundary
+++ /dev/null
@@ -1,68 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       polyBoundaryMesh;
-    location    "constant/polyMesh";
-    object      boundary;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-7
-(
-    leftWall
-    {
-        type            wall;
-        nFaces          50;
-        startFace       4419;
-    }
-    rightWall
-    {
-        type            wall;
-        nFaces          50;
-        startFace       4469;
-    }
-    lowerWall
-    {
-        type            wall;
-        nFaces          62;
-        startFace       4519;
-    }
-    atmosphere
-    {
-        type            patch;
-        nFaces          46;
-        startFace       4581;
-    }
-    porous_half0
-    {
-        type            cyclic;
-        nFaces          13;
-        startFace       4627;
-        matchTolerance  0.0001;
-        neighbourPatch  porous_half1;
-    }
-    porous_half1
-    {
-        type            cyclic;
-        nFaces          13;
-        startFace       4640;
-        matchTolerance  0.0001;
-        neighbourPatch  porous_half0;
-    }
-    defaultFaces
-    {
-        type            empty;
-        nFaces          4536;
-        startFace       4653;
-    }
-)
-
-// ************************************************************************* //
diff --git a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/controlDict b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/controlDict
index e433c1ed67..bf6c452b46 100644
--- a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/controlDict
+++ b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/controlDict
@@ -55,6 +55,8 @@ maxDeltaT       1;
 libs
 (
     "libturbulenceDerivedFvPatchFields.so"
+    "libincompressibleTurbulenceModel.so"
+    "libincompressibleRASModels.so"
 );
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/createBafflesDict b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/createBafflesDict
new file mode 100644
index 0000000000..46e0c66bf8
--- /dev/null
+++ b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/createBafflesDict
@@ -0,0 +1,78 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      createBafflesDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Whether to convert internal faces only (so leave boundary faces intact).
+// This is only relevant if your face selection type can pick up boundary
+// faces.
+internalFacesOnly true;
+
+
+// Baffles to create.
+baffles
+{
+    cyclicFaces
+    {
+        //- Select faces and orientation through a searchableSurface
+        type        searchableSurface;
+        surface     searchablePlate;
+        origin      (0.3042 0.0493 -100);
+        span        (0 0.1584  200);
+
+        patches
+        {
+            master
+            {
+                //- Master side patch
+
+                name            porous_half0;
+                type            cyclic;
+                neighbourPatch  porous_half1;
+
+                //- Optional override of added patchfields. If not specified
+                //  any added patchfields are of type calculated.
+                patchFields
+                {
+                    p_rgh
+                    {
+                        type            porousBafflePressure;
+                        patchType       cyclic;
+                        D               700;
+                        I               500;
+                        length          1.05;
+                        jump            uniform 0;
+                        value           uniform 0;
+                    }
+                }
+            }
+            slave
+            {
+                //- Slave side patch
+
+                name            porous_half1;
+                type            cyclic;
+                neighbourPatch  porous_half0;
+
+                patchFields
+                {
+                    ${...master.patchFields}
+                }
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //