Merge branch 'master' into dsmc

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H

@@ -41,24 +41,6 @@
             )
         );
 
-        Info<< "    Adding to KFluid\n" << endl;
-        KFluid.set
-        (
-            i,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "K",
-                    runTime.timeName(),
-                    fluidRegions[i],
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                thermoFluid[i].Cp()*thermoFluid[i].alpha()
-            )
-        );
-
         Info<< "    Adding to UFluid\n" << endl;
         UFluid.set
         (
@@ -129,6 +111,24 @@
             ).ptr()
         );
 
+        Info<< "    Adding to KFluid\n" << endl;
+        KFluid.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "K",
+                    runTime.timeName(),
+                    fluidRegions[i],
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                thermoFluid[i].Cp()*turbulence[i].alphaEff()
+            )
+        );
+
         Info<< "    Adding to DpDtFluid\n" << endl;
         DpDtFluid.set
         (

applications/test/dictionary/calcEntry/calcEntry.H

@@ -64,7 +64,7 @@ class calcEntry
 public:
 
     //- Runtime type information
-    TypeName("calc");
+    ClassName("calc");
 
 
     // Member Functions

applications/utilities/mesh/conversion/writeMeshObj/writeMeshObj.C

@@ -32,7 +32,8 @@ Description
 
     patch_YYY_XXX.obj : all face centres of patch YYY
 
-    Optional: patch faces (as polygons) : patchFaces_YYY_XXX.obj
+    Optional: - patch faces (as polygons) : patchFaces_YYY_XXX.obj
+              - non-manifold edges : patchEdges_YYY_XXX.obj
 
 \*---------------------------------------------------------------------------*/
 
@@ -51,7 +52,7 @@ using namespace Foam;
 
 void writeOBJ(const point& pt, Ostream& os)
 {
-    os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
+    os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
 }
 
 // All edges of mesh
@@ -75,8 +76,7 @@ void writePoints(const polyMesh& mesh, const fileName& timeName)
     {
         const edge& e = mesh.edges()[edgeI];
 
-        pointStream << "l " << e.start() + 1 << ' ' << e.end() + 1
+        pointStream << "l " << e.start() + 1 << ' ' << e.end() + 1 << nl;
-            << endl;
     }
 }
 
@@ -277,7 +277,47 @@ void writePatchFaces
             {
                 patchFaceStream << ' ' << f[fp]+1;
             }
-            patchFaceStream << endl;
+            patchFaceStream << nl;
+        }
+    }
+}
+
+
+void writePatchBoundaryEdges
+(
+    const polyMesh& mesh,
+    const fileName& timeName
+)
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        fileName edgeFile
+        (
+            mesh.time().path()
+          / "patchEdges_" + pp.name() + '_' + timeName + ".obj"
+        );
+
+        Info << "Writing patch edges to " << edgeFile << endl;
+
+        OFstream patchEdgeStream(edgeFile);
+
+        forAll(pp.localPoints(), pointI)
+        {
+            writeOBJ(pp.localPoints()[pointI], patchEdgeStream);
+        }
+
+        for (label edgeI = pp.nInternalEdges(); edgeI < pp.nEdges(); edgeI++)
+        {
+            if (pp.edgeFaces()[edgeI].size() == 1)
+            {
+                const edge& e = pp.edges()[edgeI];
+
+                patchEdgeStream<< "l " << e[0]+1 << ' ' << e[1]+1 << nl;
+            }
         }
     }
 }
@@ -339,6 +379,7 @@ int main(int argc, char *argv[])
 {
     timeSelector::addOptions();
     argList::validOptions.insert("patchFaces", "");
+    argList::validOptions.insert("patchEdges", "");
     argList::validOptions.insert("cell", "cellI");
     argList::validOptions.insert("face", "faceI");
     argList::validOptions.insert("point", "pointI");
@@ -351,6 +392,7 @@ int main(int argc, char *argv[])
     runTime.functionObjects().off();
 
     bool patchFaces = args.optionFound("patchFaces");
+    bool patchEdges = args.optionFound("patchEdges");
     bool doCell     = args.optionFound("cell");
     bool doPoint    = args.optionFound("point");
     bool doFace     = args.optionFound("face");
@@ -381,19 +423,23 @@ int main(int argc, char *argv[])
             {
                 writePatchFaces(mesh, runTime.timeName());
             }
-            else if (doCell)
+            if (patchEdges)
+            {
+                writePatchBoundaryEdges(mesh, runTime.timeName());
+            }
+            if (doCell)
             {
                 label cellI = args.optionRead<label>("cell");
 
                 writePoints(mesh, cellI, runTime.timeName());
             }
-            else if (doPoint)
+            if (doPoint)
             {
                 label pointI = args.optionRead<label>("point");
 
                 writePointCells(mesh, pointI, runTime.timeName());
             }
-            else if (doFace)
+            if (doFace)
             {
                 label faceI = args.optionRead<label>("face");
 
@@ -415,7 +461,7 @@ int main(int argc, char *argv[])
 
                 meshTools::writeOBJ(str, faceList(1, f), mesh.points());
             }
-            else if (doCellSet)
+            if (doCellSet)
             {
                 word setName(args.option("cellSet"));
 
@@ -427,7 +473,7 @@ int main(int argc, char *argv[])
                 writePoints(mesh, cells.toc(), runTime.timeName());
 
             }
-            else if (doFaceSet)
+            if (doFaceSet)
             {
                 word setName(args.option("faceSet"));
 
@@ -458,7 +504,16 @@ int main(int argc, char *argv[])
                     faces.toc()
                 );
             }
-            else
+            else if
+            (
+                !patchFaces 
+             && !patchEdges
+             && !doCell
+             && !doPoint
+             && !doFace
+             && !doCellSet
+             && !doFaceSet
+            )
             {
                 // points & edges
                 writePoints(mesh, runTime.timeName());

applications/utilities/mesh/generation/extrudeMesh/extrudeProperties

@@ -71,7 +71,8 @@ linearNormalCoeffs
 
 linearDirectionCoeffs
 {
-    direction       (0 0 1);
+    direction       (0 1 0);
+    thickness       0.05;
 }
 
 linearRadialCoeffs
@@ -89,7 +90,7 @@ sigmaRadialCoeffs
 
 // Do front and back need to be merged? Usually only makes sense for 360
 // degree wedges.
-mergeFaces true;
+mergeFaces false;   //true;
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/utilities/mesh/manipulation/cellSet/cellSet.C

@@ -37,37 +37,6 @@ using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Copy set
-void backup
-(
-    const polyMesh& mesh,
-    const word& fromName,
-    const topoSet& fromSet,
-    const word& toName
-)
-{
-    Info<< "Backing up " << fromName << " into " << toName << endl;
-
-    topoSet backupSet(mesh, toName, fromSet);
-
-    backupSet.write();
-}
-
-
-// Read and copy set
-void backup
-(
-    const polyMesh& mesh,
-    const word& fromName,
-    const word& toName
-)
-{
-    topoSet fromSet(mesh, fromName, IOobject::READ_IF_PRESENT);
-
-    backup(mesh, fromName, fromSet, toName);
-}
-
-
 // Main program:
 
 int main(int argc, char *argv[])
@@ -114,8 +83,6 @@ int main(int argc, char *argv[])
     {
         r = IOobject::NO_READ;
 
-        backup(mesh, setName, setName + "_old");
-
         currentSetPtr.reset
         (
             new cellSet
@@ -151,7 +118,7 @@ int main(int argc, char *argv[])
     if ((r == IOobject::MUST_READ) && (action != topoSetSource::LIST))
     {
         // currentSet has been read so can make copy.
-        backup(mesh, setName, currentSet, setName + "_old");
+        //backup(mesh, setName, currentSet, setName + "_old");
     }
 
     if (action == topoSetSource::CLEAR)
@@ -173,7 +140,16 @@ int main(int argc, char *argv[])
         forAll(topoSetSources, topoSetSourceI)
         {
             // Backup current set.
-            topoSet oldSet(mesh, currentSet.name() + "_old2", currentSet);
+            autoPtr<topoSet> oldSet
+            (
+                topoSet::New
+                (
+                    currentSet.type(),
+                    mesh,
+                    currentSet.name() + "_old2",
+                    currentSet
+                )
+            );
 
             currentSet.clear();
 

applications/utilities/mesh/manipulation/checkMesh/checkMesh.C

@@ -47,7 +47,7 @@ using namespace Foam;
 
 int main(int argc, char *argv[])
 {
-    timeSelector::addOptions(false);   // no constant
+    timeSelector::addOptions();
 #   include "addRegionOption.H"
     argList::validOptions.insert("noTopology", "");
     argList::validOptions.insert("allGeometry", "");

applications/utilities/mesh/manipulation/faceSet/faceSet.C

@@ -37,37 +37,6 @@ using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Copy set
-void backup
-(
-    const polyMesh& mesh,
-    const word& fromName,
-    const topoSet& fromSet,
-    const word& toName
-)
-{
-    Info<< "Backing up " << fromName << " into " << toName << endl;
-
-    topoSet backupSet(mesh, toName, fromSet);
-
-    backupSet.write();
-}
-
-
-// Read and copy set
-void backup
-(
-    const polyMesh& mesh,
-    const word& fromName,
-    const word& toName
-)
-{
-    topoSet fromSet(mesh, fromName, IOobject::READ_IF_PRESENT);
-
-    backup(mesh, fromName, fromSet, toName);
-}
-
-
 // Main program:
 
 int main(int argc, char *argv[])
@@ -114,8 +83,6 @@ int main(int argc, char *argv[])
     {
         r = IOobject::NO_READ;
 
-        backup(mesh, setName, setName + "_old");
-
         currentSetPtr.reset
         (
             new faceSet
@@ -151,7 +118,7 @@ int main(int argc, char *argv[])
     if ((r == IOobject::MUST_READ) && (action != topoSetSource::LIST))
     {
         // currentSet has been read so can make copy.
-        backup(mesh, setName, currentSet, setName + "_old");
+        //backup(mesh, setName, currentSet, setName + "_old");
     }
 
     if (action == topoSetSource::CLEAR)
@@ -173,7 +140,16 @@ int main(int argc, char *argv[])
         forAll(topoSetSources, topoSetSourceI)
         {
             // Backup current set.
-            topoSet oldSet(mesh, currentSet.name() + "_old2", currentSet);
+            autoPtr<topoSet> oldSet
+            (
+                topoSet::New
+                (
+                    currentSet.type(),
+                    mesh,
+                    currentSet.name() + "_old2",
+                    currentSet
+                )
+            );
 
             currentSet.clear();
 

applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C

@@ -49,7 +49,7 @@ int main(int argc, char *argv[])
     (
         IOobject
         (
-            mergePolyMesh::defaultRegion,
+            masterRegion,
             runTimeMaster.timeName(),
             runTimeMaster
         )
@@ -63,7 +63,7 @@ int main(int argc, char *argv[])
     (
         IOobject
         (
-            mergePolyMesh::defaultRegion,
+            addRegion,
             runTimeToAdd.timeName(),
             runTimeToAdd
         )

applications/utilities/mesh/manipulation/mergeMeshes/setRoots.H

@@ -2,9 +2,11 @@
 
     argList::validArgs.append("master root");
     argList::validArgs.append("master case");
+    argList::validOptions.insert("masterRegion", "name");
 
     argList::validArgs.append("root to add");
     argList::validArgs.append("case to add");
+    argList::validOptions.insert("addRegion", "name");
 
     argList args(argc, argv);
 
@@ -15,9 +17,15 @@
 
     fileName rootDirMaster(args.additionalArgs()[0]);
     fileName caseDirMaster(args.additionalArgs()[1]);
+    word masterRegion = polyMesh::defaultRegion;
+    args.optionReadIfPresent("masterRegion", masterRegion);
 
     fileName rootDirToAdd(args.additionalArgs()[2]);
     fileName caseDirToAdd(args.additionalArgs()[3]);
+    word addRegion = polyMesh::defaultRegion;
+    args.optionReadIfPresent("addRegion", addRegion);
 
-    Info<< "Master:      " << rootDirMaster << " " << caseDirMaster << nl
+    Info<< "Master:      " << rootDirMaster << " " << caseDirMaster
-        << "mesh to add: " << rootDirToAdd << " " << caseDirToAdd << endl;
+        << "  region " << masterRegion << nl
+        << "mesh to add: " << rootDirToAdd << " " << caseDirToAdd
+        << "  region " << addRegion << endl;

applications/utilities/mesh/manipulation/pointSet/pointSet.C

@@ -37,37 +37,6 @@ using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Copy set
-void backup
-(
-    const polyMesh& mesh,
-    const word& fromName,
-    const topoSet& fromSet,
-    const word& toName
-)
-{
-    Info<< "Backing up " << fromName << " into " << toName << endl;
-
-    topoSet backupSet(mesh, toName, fromSet);
-
-    backupSet.write();
-}
-
-
-// Read and copy set
-void backup
-(
-    const polyMesh& mesh,
-    const word& fromName,
-    const word& toName
-)
-{
-    topoSet fromSet(mesh, fromName, IOobject::READ_IF_PRESENT);
-
-    backup(mesh, fromName, fromSet, toName);
-}
-
-
 // Main program:
 
 int main(int argc, char *argv[])
@@ -114,8 +83,6 @@ int main(int argc, char *argv[])
     {
         r = IOobject::NO_READ;
 
-        backup(mesh, setName, setName + "_old");
-
         currentSetPtr.reset
         (
             new pointSet
@@ -151,7 +118,7 @@ int main(int argc, char *argv[])
     if ((r == IOobject::MUST_READ) && (action != topoSetSource::LIST))
     {
         // currentSet has been read so can make copy.
-        backup(mesh, setName, currentSet, setName + "_old");
+        //backup(mesh, setName, currentSet, setName + "_old");
     }
 
     if (action == topoSetSource::CLEAR)
@@ -173,7 +140,16 @@ int main(int argc, char *argv[])
         forAll(topoSetSources, topoSetSourceI)
         {
             // Backup current set.
-            topoSet oldSet(mesh, currentSet.name() + "_old2", currentSet);
+            autoPtr<topoSet> oldSet
+            (
+                topoSet::New
+                (
+                    currentSet.type(),
+                    mesh,
+                    currentSet.name() + "_old2",
+                    currentSet
+                )
+            );
 
             currentSet.clear();
 

applications/utilities/mesh/manipulation/setSet/setSet.C

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Description
-    Manipulate a cell/face/point set interactively.
+    Manipulate a cell/face/point/ set or zone interactively.
 
 \*---------------------------------------------------------------------------*/
 
@@ -42,6 +42,9 @@ Description
 #include "writePatch.H"
 #include "writePointSet.H"
 #include "IOobjectList.H"
+#include "cellZoneSet.H"
+#include "faceZoneSet.H"
+#include "pointZoneSet.H"
 
 #include <stdio.h>
 
@@ -82,6 +85,7 @@ Istream& selectStream(Istream* is0Ptr, Istream* is1Ptr)
 // Copy set
 void backup
 (
+    const word& setType,
     const polyMesh& mesh,
     const word& fromName,
     const topoSet& fromSet,
@@ -92,9 +96,7 @@ void backup
     {
         Pout<< "    Backing up " << fromName << " into " << toName << endl;
 
-        topoSet backupSet(mesh, toName, fromSet);
+        topoSet::New(setType, mesh, toName, fromSet)().write();
-
-        backupSet.write();
     }
 }
 
@@ -102,14 +104,21 @@ void backup
 // Read and copy set
 void backup
 (
+    const word& setType,
     const polyMesh& mesh,
     const word& fromName,
     const word& toName
 )
 {
-    topoSet fromSet(mesh, fromName, IOobject::READ_IF_PRESENT);
+    autoPtr<topoSet> fromSet = topoSet::New
+    (
+        setType,
+        mesh,
+        fromName,
+        IOobject::READ_IF_PRESENT
+    );
 
-    backup(mesh, fromName, fromSet, toName);
+    backup(setType, mesh, fromName, fromSet(), toName);
 }
 
 
@@ -241,11 +250,13 @@ void printHelp(Ostream& os)
         << "A set command should be of the following form" << endl
         << endl
         << "    cellSet|faceSet|pointSet <setName> <action> <source>"
-        << endl << endl
+        << endl
+        << endl
         << "The <action> is one of" << endl
         << "    list            - prints the contents of the set" << endl
         << "    clear           - clears the set" << endl
         << "    invert          - inverts the set" << endl
+        << "    remove          - remove the set" << endl
         << "    new <source>    - sets to set to the source set" << endl
         << "    add <source>    - adds all elements from the source set" << endl
         << "    delete <source> - deletes      ,," << endl
@@ -270,6 +281,20 @@ void printHelp(Ostream& os)
         << "    cellSet c0 add pointToCell p0 any" << endl
         << "List set:" << endl
         << "    cellSet c0 list" << endl
+        << endl
+        << "Zones can be set using zoneSets from corresponding sets:" << endl
+        << "    cellZoneSet c0Zone new setToZone c0" << endl
+        << "    faceZoneSet f0Zone new setToZone f0" << endl 
+        << endl
+        << "or if orientation is important:" << endl
+        << "    faceZoneSet f0Zone new setsToZone f0 c0" << endl 
+        << endl
+        << "ZoneSets can be manipulated using the general actions:" << endl
+        << "    list            - prints the contents of the set" << endl
+        << "    clear           - clears the set" << endl
+        << "    invert          - inverts the set (undefined orientation)"
+        << endl
+        << "    remove          - remove the set" << endl
         << endl;
 }
 
@@ -312,10 +337,126 @@ void printAllSets(const polyMesh& mesh, Ostream& os)
             os  << '\t' << set.name() << "\tsize:" << set.size() << endl;
         }
     }
+
+    const cellZoneMesh& cellZones = mesh.cellZones();
+    if (cellZones.size())
+    {
+        os  << "cellZones:" << endl;
+        forAll(cellZones, i)
+        {
+            const cellZone& zone = cellZones[i];
+            os  << '\t' << zone.name() << "\tsize:" << zone.size() << endl;
+        }
+    }
+    const faceZoneMesh& faceZones = mesh.faceZones();
+    if (faceZones.size())
+    {
+        os  << "faceZones:" << endl;
+        forAll(faceZones, i)
+        {
+            const faceZone& zone = faceZones[i];
+            os  << '\t' << zone.name() << "\tsize:" << zone.size() << endl;
+        }
+    }
+    const pointZoneMesh& pointZones = mesh.pointZones();
+    if (pointZones.size())
+    {
+        os  << "pointZones:" << endl;
+        forAll(pointZones, i)
+        {
+            const pointZone& zone = pointZones[i];
+            os  << '\t' << zone.name() << "\tsize:" << zone.size() << endl;
+        }
+    }
+
     os  << endl;
 }
 
 
+template<class ZoneType>
+void removeZone
+(
+    ZoneMesh<ZoneType, polyMesh>& zones,
+    const word& setName
+)
+{
+    label zoneID = zones.findZoneID(setName);
+
+    if (zoneID != -1)
+    {
+        Info<< "Removing zone " << setName << " at index " << zoneID << endl;
+        // Shuffle to last position
+        labelList oldToNew(zones.size());
+        label newI = 0;
+        forAll(oldToNew, i)
+        {
+            if (i != zoneID)
+            {
+                oldToNew[i] = newI++;
+            }
+        }
+        oldToNew[zoneID] = newI;
+        zones.reorder(oldToNew);
+        // Remove last element
+        zones.setSize(zones.size()-1);
+        zones.clearAddressing();
+        zones.write();
+    }
+}
+
+
+// Physically remove a set
+void removeSet
+(
+    const polyMesh& mesh,
+    const word& setType,
+    const word& setName
+)
+{
+    // Remove the file
+    IOobjectList objects
+    (
+        mesh,
+        mesh.pointsInstance(),
+        polyMesh::meshSubDir/"sets"
+    );
+
+    if (objects.found(setName))
+    {
+        // Remove file
+        fileName object = objects[setName]->objectPath();
+        Info<< "Removing file " << object << endl;
+        rm(object);
+    }
+
+    // See if zone
+    if (setType == cellZoneSet::typeName)
+    {
+        removeZone
+        (
+            const_cast<cellZoneMesh&>(mesh.cellZones()),
+            setName
+        );
+    }
+    else if (setType == faceZoneSet::typeName)
+    {
+        removeZone
+        (
+            const_cast<faceZoneMesh&>(mesh.faceZones()),
+            setName
+        );
+    }
+    else if (setType == pointZoneSet::typeName)
+    {
+        removeZone
+        (
+            const_cast<pointZoneMesh&>(mesh.pointZones()),
+            setName
+        );
+    }
+}
+
+
 // Read command and execute. Return true if ok, false otherwise.
 bool doCommand
 (
@@ -358,38 +499,29 @@ bool doCommand
 
         IOobject::readOption r;
 
-        if
+        if (action == topoSetSource::REMOVE)
+        {
+            removeSet(mesh, setType, setName);
+        }
+        else if
         (
             (action == topoSetSource::NEW)
          || (action == topoSetSource::CLEAR)
         )
         {
             r = IOobject::NO_READ;
-
-            //backup(mesh, setName, setName + "_old");
-
             currentSetPtr = topoSet::New(setType, mesh, setName, typSize);
         }
         else
         {
             r = IOobject::MUST_READ;
-
             currentSetPtr = topoSet::New(setType, mesh, setName, r);
-
             topoSet& currentSet = currentSetPtr();
-
             // Presize it according to current mesh data.
             currentSet.resize(max(currentSet.size(), typSize));
         }
 
-        if (currentSetPtr.empty())
+        if (currentSetPtr.valid())
-        {
-            Pout<< "    Cannot construct/load set "
-                << topoSet::localPath(mesh, setName) << endl;
-
-            ok = false;
-        }
-        else
         {
             topoSet& currentSet = currentSetPtr();
 
@@ -398,12 +530,6 @@ bool doCommand
                 << "  Action:" << actionName
                 << endl;
 
-            //if ((r == IOobject::MUST_READ) && (action != topoSetSource::LIST))
-            //{
-            //    // currentSet has been read so can make copy.
-            //    backup(mesh, setName, currentSet, setName + "_old");
-            //}
-
             switch (action)
             {
                 case topoSetSource::CLEAR:
@@ -437,15 +563,18 @@ bool doCommand
                         );
 
                         // Backup current set.
-                        topoSet oldSet
+                        autoPtr<topoSet> oldSet
                         (
-                            mesh,
+                            topoSet::New
-                            currentSet.name() + "_old2",
+                            (
-                            currentSet
+                                setType,
+                                mesh,
+                                currentSet.name() + "_old2",
+                                currentSet
+                            )
                         );
 
                         currentSet.clear();
-                        currentSet.resize(oldSet.size());
                         setSource().applyToSet(topoSetSource::NEW, currentSet);
 
                         // Combine new value of currentSet with old one.
@@ -547,7 +676,8 @@ enum commandStatus
 {
     QUIT,           // quit program
     INVALID,        // token is not a valid set manipulation command
-    VALID           // ,,    is a valid     ,,
+    VALIDSETCMD,    // ,,    is a valid     ,,
+    VALIDZONECMD    // ,,    is a valid     zone      ,,
 };
 
 
@@ -654,7 +784,16 @@ commandStatus parseType
      || setType == "pointSet"
     )
     {
-        return VALID;
+        return VALIDSETCMD;
+    }
+    else if
+    (
+        setType == "cellZoneSet"
+     || setType == "faceZoneSet"
+     || setType == "pointZoneSet"
+    )
+    {
+        return VALIDZONECMD;
     }
     else
     {
@@ -664,7 +803,7 @@ commandStatus parseType
             ", IStringStream&)"
         )   << "Illegal command " << setType << endl
             << "Should be one of 'help', 'list', 'time' or a set type :"
-            << " 'cellSet', 'faceSet', 'pointSet'"
+            << " 'cellSet', 'faceSet', 'pointSet', 'faceZoneSet'"
             << endl;
 
         return INVALID;
@@ -682,7 +821,7 @@ commandStatus parseAction(const word& actionName)
         {
             (void)topoSetSource::toAction(actionName);
 
-            stat = VALID;
+            stat = VALIDSETCMD;
         }
         catch (Foam::IOerror& fIOErr)
         {
@@ -724,6 +863,10 @@ int main(int argc, char *argv[])
     // Print some mesh info
     printMesh(runTime, mesh);
 
+    // Print current sets
+    printAllSets(mesh, Pout);
+
+
 
     std::ifstream* fileStreamPtr(NULL);
 
@@ -745,11 +888,10 @@ int main(int argc, char *argv[])
 #if READLINE != 0
     else if (!read_history(historyFile))
     {
-        Info<< "Successfully read history from " << historyFile << endl;
+        Pout<< "Successfully read history from " << historyFile << endl;
     }
 #endif
 
-
     Pout<< "Please type 'help', 'quit' or a set command after prompt." << endl;
 
     bool ok = true;
@@ -816,12 +958,12 @@ int main(int argc, char *argv[])
 
         IStringStream is(rawLine + ' ');
 
-        // Type: cellSet, faceSet, pointSet
+        // Type: cellSet, faceSet, pointSet, faceZoneSet
         is >> setType;
 
         stat = parseType(runTime, mesh, setType, is);
 
-        if (stat == VALID)
+        if (stat == VALIDSETCMD || stat == VALIDZONECMD)
         {
             if (is >> setName)
             {
@@ -831,14 +973,13 @@ int main(int argc, char *argv[])
                 }
             }
         }
-
         ok = true;
 
         if (stat == QUIT)
         {
             break;
         }
-        else if (stat == VALID)
+        else if (stat == VALIDSETCMD || stat == VALIDZONECMD)
         {
             ok = doCommand(mesh, setType, setName, actionName, writeVTK, is);
         }

applications/utilities/mesh/manipulation/splitMesh/splitMesh.C

@@ -270,6 +270,10 @@ int main(int argc, char *argv[])
     {
         mesh.setInstance(oldInstance);
     }
+    else
+    {
+        mesh.setInstance(runTime.timeName());
+    }
 
     Info<< "Writing mesh to " << runTime.timeName() << endl;
     if (!mesh.write())

applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C

@@ -43,6 +43,9 @@ Description
     the largest matching region (-sloppyCellZones). This will accept any
     region that covers more than 50% of the zone. It has to be a subset
     so cannot have any cells in any other zone.
+    - useCellZonesOnly does not do a walk and uses the cellZones only. Use
+    this if you don't mind having disconnected domains in a single region.
+    This option requires all cells to be in one (and one only) region.
 
 \*---------------------------------------------------------------------------*/
 
@@ -609,7 +612,7 @@ void getInterfaceSizes
 // Create mesh for region.
 autoPtr<mapPolyMesh> createRegionMesh
 (
-    const regionSplit& cellRegion,
+    const labelList& cellRegion,
     const EdgeMap<label>& interfaceToPatch,
     const fvMesh& mesh,
     const label regionI,
@@ -766,7 +769,7 @@ autoPtr<mapPolyMesh> createRegionMesh
 void createAndWriteRegion
 (
     const fvMesh& mesh,
-    const regionSplit& cellRegion,
+    const labelList& cellRegion,
     const wordList& regionNames,
     const EdgeMap<label>& interfaceToPatch,
     const label regionI,
@@ -1005,7 +1008,8 @@ void createAndWriteRegion
 EdgeMap<label> addRegionPatches
 (
     fvMesh& mesh,
-    const regionSplit& cellRegion,
+    const labelList& cellRegion,
+    const label nCellRegions,
     const edgeList& interfaces,
     const EdgeMap<label>& interfaceSizes,
     const wordList& regionNames
@@ -1016,7 +1020,7 @@ EdgeMap<label> addRegionPatches
 
     Info<< nl << "Adding patches" << nl << endl;
 
-    EdgeMap<label> interfaceToPatch(cellRegion.nRegions());
+    EdgeMap<label> interfaceToPatch(nCellRegions);
 
     forAll(interfaces, interI)
     {
@@ -1108,13 +1112,14 @@ EdgeMap<label> addRegionPatches
 label findCorrespondingRegion
 (
     const labelList& existingZoneID,    // per cell the (unique) zoneID
-    const regionSplit& cellRegion,
+    const labelList& cellRegion,
+    const label nCellRegions,
     const label zoneI,
     const label minOverlapSize
 )
 {
     // Per region the number of cells in zoneI
-    labelList cellsInZone(cellRegion.nRegions(), 0);
+    labelList cellsInZone(nCellRegions, 0);
 
     forAll(cellRegion, cellI)
     {
@@ -1162,7 +1167,8 @@ label findCorrespondingRegion
 //(
 //    const cellZoneMesh& cellZones,
 //    const labelList& existingZoneID,    // per cell the (unique) zoneID
-//    const regionSplit& cellRegion,
+//    const labelList& cellRegion,
+//    const label nCellRegions,
 //    const label zoneI
 //)
 //{
@@ -1227,6 +1233,7 @@ label findCorrespondingRegion
 int main(int argc, char *argv[])
 {
     argList::validOptions.insert("cellZones", "");
+    argList::validOptions.insert("cellZonesOnly", "");
     argList::validOptions.insert("blockedFaces", "faceSet");
     argList::validOptions.insert("makeCellZones", "");
     argList::validOptions.insert("largestOnly", "");
@@ -1249,13 +1256,14 @@ int main(int argc, char *argv[])
             << blockedFacesName << nl << endl;
     }
 
-    bool makeCellZones   = args.optionFound("makeCellZones");
+    bool makeCellZones    = args.optionFound("makeCellZones");
-    bool largestOnly     = args.optionFound("largestOnly");
+    bool largestOnly      = args.optionFound("largestOnly");
-    bool insidePoint     = args.optionFound("insidePoint");
+    bool insidePoint      = args.optionFound("insidePoint");
-    bool useCellZones    = args.optionFound("cellZones");
+    bool useCellZones     = args.optionFound("cellZones");
-    bool overwrite       = args.optionFound("overwrite");
+    bool useCellZonesOnly = args.optionFound("cellZonesOnly");
-    bool detectOnly      = args.optionFound("detectOnly");
+    bool overwrite        = args.optionFound("overwrite");
-    bool sloppyCellZones = args.optionFound("sloppyCellZones");
+    bool detectOnly       = args.optionFound("detectOnly");
+    bool sloppyCellZones  = args.optionFound("sloppyCellZones");
 
     if (insidePoint && largestOnly)
     {
@@ -1370,13 +1378,37 @@ int main(int argc, char *argv[])
     }
 
 
-    // Do the topological walk to determine regions
+    // Determine per cell the region it belongs to
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // regionSplit is the labelList with the region per cell.
+    // cellRegion is the labelList with the region per cell.
-    regionSplit cellRegion(mesh, blockedFace);
+    labelList cellRegion;
+    label nCellRegions = 0;
+    if (useCellZonesOnly)
+    {
+        label unzonedCellI = findIndex(zoneID, -1);
+        if (unzonedCellI != -1)
+        {
+            FatalErrorIn(args.executable())
+                << "For the cellZonesOnly option all cells "
+                << "have to be in a cellZone." << endl
+                << "Cell " << unzonedCellI
+                << " at" << mesh.cellCentres()[unzonedCellI]
+                << " is not in a cellZone. There might be more unzoned cells."
+                << exit(FatalError);
+        }
+        cellRegion = zoneID;
+        nCellRegions = gMax(cellRegion)+1;
+    }
+    else
+    {
+        // Do a topological walk to determine regions
+        regionSplit regions(mesh, blockedFace);
+        nCellRegions = regions.nRegions();
+        cellRegion.transfer(regions);
+    }
 
-    Info<< endl << "Number of regions:" << cellRegion.nRegions() << nl << endl;
+    Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
 
 
     // Write to manual decomposition option
@@ -1429,7 +1461,7 @@ int main(int argc, char *argv[])
     // Sizes per region
     // ~~~~~~~~~~~~~~~~
 
-    labelList regionSizes(cellRegion.nRegions(), 0);
+    labelList regionSizes(nCellRegions, 0);
 
     forAll(cellRegion, cellI)
     {
@@ -1489,9 +1521,9 @@ int main(int argc, char *argv[])
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Region per zone
-    labelList regionToZone(cellRegion.nRegions(), -1);
+    labelList regionToZone(nCellRegions, -1);
     // Name of region
-    wordList regionNames(cellRegion.nRegions());
+    wordList regionNames(nCellRegions);
     // Zone to region
     labelList zoneToRegion(cellZones.size(), -1);
 
@@ -1506,6 +1538,7 @@ int main(int argc, char *argv[])
             (
                 zoneID,
                 cellRegion,
+                nCellRegions,
                 zoneI,
                 label(0.5*zoneSizes[zoneI]) // minimum overlap
             );
@@ -1532,6 +1565,7 @@ int main(int argc, char *argv[])
             (
                 zoneID,
                 cellRegion,
+                nCellRegions,
                 zoneI,
                 1               // minimum overlap
             );
@@ -1660,7 +1694,7 @@ int main(int argc, char *argv[])
     mesh.clearOut();
 
 
-    if (cellRegion.nRegions() == 1)
+    if (nCellRegions == 1)
     {
         Info<< "Only one region. Doing nothing." << endl;
     }
@@ -1671,7 +1705,7 @@ int main(int argc, char *argv[])
 
         // Check if region overlaps with existing zone. If so keep.
 
-        for (label regionI = 0; regionI < cellRegion.nRegions(); regionI++)
+        for (label regionI = 0; regionI < nCellRegions; regionI++)
         {
             label zoneI = regionToZone[regionI];
 
@@ -1759,6 +1793,7 @@ int main(int argc, char *argv[])
             (
                 mesh,
                 cellRegion,
+                nCellRegions,
                 interfaces,
                 interfaceSizes,
                 regionNames
@@ -1837,7 +1872,7 @@ int main(int argc, char *argv[])
         else
         {
             // Split all
-            for (label regionI = 0; regionI < cellRegion.nRegions(); regionI++)
+            for (label regionI = 0; regionI < nCellRegions; regionI++)
             {
                 Info<< nl
                     << "Region " << regionI << nl

applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C

@@ -334,7 +334,7 @@ int main(int argc, char *argv[])
     }
     else
     {
-        mesh.setInstance(oldInstance);
+        subsetter.subMesh().setInstance(oldInstance);
     }
 
     Info<< "Writing subsetted mesh and fields to time " << runTime.timeName()

applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C

@@ -34,6 +34,7 @@ License
 #include "Map.H"
 #include "globalMeshData.H"
 #include "DynamicList.H"
+#include "fvFieldDecomposer.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -306,12 +307,29 @@ bool domainDecomposition::writeDecomposition()
 
         forAll (curPatchSizes, patchi)
         {
+            // Get the face labels consistent with the field mapping
+            // (reuse the patch field mappers)
+            const polyPatch& meshPatch =
+                meshPatches[curBoundaryAddressing[patchi]];
+
+            fvFieldDecomposer::patchFieldDecomposer patchMapper
+            (
+                SubList<label>
+                (
+                    curFaceLabels,
+                    curPatchSizes[patchi],
+                    curPatchStarts[patchi]
+                ),
+                meshPatch.start()
+            );
+
+            // Map existing patches
             procPatches[nPatches] =
                 meshPatches[curBoundaryAddressing[patchi]].clone
                 (
                     procMesh.boundaryMesh(),
                     nPatches,
-                    curPatchSizes[patchi],
+                    patchMapper.directAddressing(),
                     curPatchStarts[patchi]
                 ).ptr();
 

applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C

@@ -38,13 +38,14 @@ Description
 
 int main(int argc, char *argv[])
 {
+#   include "addRegionOption.H"
     timeSelector::addOptions();
     argList::validArgs.append("fieldName");
     argList::validArgs.append("patchName");
 #   include "setRootCase.H"
 #   include "createTime.H"
     instantList timeDirs = timeSelector::select0(runTime, args);
-#   include "createMesh.H"
+#   include "createNamedMesh.H"
 
     word fieldName(args.additionalArgs()[0]);
     word patchName(args.additionalArgs()[1]);

etc/settings.csh

@@ -35,6 +35,8 @@
 alias _foamAddPath 'set path=(\!* $path)'
 # prefix to LD_LIBRARY_PATH
 alias _foamAddLib 'setenv LD_LIBRARY_PATH \!*\:${LD_LIBRARY_PATH}'
+# prefix to MANPATH
+alias _foamAddManPath 'setenv MANPATH \!*\:${MANPATH}'
 
 # location of the jobControl directory
 setenv FOAM_JOB_DIR $WM_PROJECT_INST_DIR/jobControl
@@ -114,9 +116,10 @@ case OpenFOAM:
         echo
     endif
 
-    _foamAddPath ${WM_COMPILER_DIR}/bin
+    _foamAddPath    ${WM_COMPILER_DIR}/bin
-    _foamAddLib  ${WM_COMPILER_DIR}/lib${WM_COMPILER_LIB_ARCH}
+    _foamAddLib     ${WM_COMPILER_DIR}/lib${WM_COMPILER_LIB_ARCH}
-    _foamAddLib  ${WM_COMPILER_DIR}/lib
+    _foamAddLib     ${WM_COMPILER_DIR}/lib
+    _foamAddManPath ${WM_COMPILER_DIR}/man
 
     breaksw
 endsw
@@ -136,21 +139,22 @@ case OPENMPI:
     # Tell OpenMPI where to find its install directory
     setenv OPAL_PREFIX $MPI_ARCH_PATH
 
-    _foamAddPath $MPI_ARCH_PATH/bin
+    _foamAddPath    $MPI_ARCH_PATH/bin
-    _foamAddLib  $MPI_ARCH_PATH/lib
+    _foamAddLib     $MPI_ARCH_PATH/lib
+    _foamAddManPath $MPI_ARCH_PATH/man
 
     setenv FOAM_MPI_LIBBIN $FOAM_LIBBIN/$mpi_version
     unset mpi_version
     breaksw
 
 case MPICH:
-    set mpi_version=mpich-1.2.4
+    set mpi_version=mpich2-1.1.1p1
     setenv MPI_HOME $WM_THIRD_PARTY_DIR/$mpi_version
     setenv MPI_ARCH_PATH $MPI_HOME/platforms/$WM_OPTIONS
-    setenv MPICH_ROOT $MPI_ARCH_PATH
 
-    _foamAddPath $MPI_ARCH_PATH/bin
+    _foamAddPath    $MPI_ARCH_PATH/bin
-    _foamAddLib  $MPI_ARCH_PATH/lib
+    _foamAddLib     $MPI_ARCH_PATH/lib
+    _foamAddManPath $MPI_ARCH_PATH/share/man
 
     setenv FOAM_MPI_LIBBIN $FOAM_LIBBIN/$mpi_version
     unset mpi_version
@@ -159,7 +163,6 @@ case MPICH:
 case MPICH-GM:
     setenv MPI_ARCH_PATH /opt/mpi
     setenv MPICH_PATH $MPI_ARCH_PATH
-    setenv MPICH_ROOT $MPI_ARCH_PATH
     setenv GM_LIB_PATH /opt/gm/lib64
 
     _foamAddPath $MPI_ARCH_PATH/bin
@@ -172,7 +175,6 @@ case MPICH-GM:
 case HPMPI:
     setenv MPI_HOME /opt/hpmpi
     setenv MPI_ARCH_PATH $MPI_HOME
-    setenv MPICH_ROOT=$MPI_ARCH_PATH
 
     _foamAddPath $MPI_ARCH_PATH/bin
 

etc/settings.sh

@@ -51,6 +51,16 @@ _foamAddLib()
     done
 }
 
+# prefix to MANPATH
+_foamAddManPath()
+{
+    while [ $# -ge 1 ]
+    do
+        export MANPATH=$1:$MANPATH
+        shift
+    done
+}
+
 
 # location of the jobControl directory
 export FOAM_JOB_DIR=$WM_PROJECT_INST_DIR/jobControl
@@ -90,7 +100,7 @@ _foamAddLib  $FOAM_LIBBIN $FOAM_SITE_LIBBIN $FOAM_USER_LIBBIN
 
 # Compiler settings
 # ~~~~~~~~~~~~~~~~~
-unset compilerBin compilerLib
+unset compilerBin compilerLib compilerMan
 
 # Select compiler installation
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -133,16 +143,18 @@ OpenFOAM)
 
     compilerBin=$WM_COMPILER_DIR/bin
     compilerLib=$WM_COMPILER_DIR/lib$WM_COMPILER_LIB_ARCH:$WM_COMPILER_DIR/lib
+    compilerMan=$WM_COMPILER_DIR/man
     ;;
 esac
 
 if [ -d "$compilerBin" ]
 then
-    _foamAddPath $compilerBin
+    _foamAddPath    $compilerBin
-    _foamAddLib  $compilerLib
+    _foamAddLib     $compilerLib
+    _foamAddManPath $compilerMan
 fi
 
-unset compilerBin compilerLib compilerInstall
+unset compilerBin compilerLib compilerMan compilerInstall
 
 # Communications library
 # ~~~~~~~~~~~~~~~~~~~~~~
@@ -158,21 +170,22 @@ OPENMPI)
     # Tell OpenMPI where to find its install directory
     export OPAL_PREFIX=$MPI_ARCH_PATH
 
-    _foamAddPath $MPI_ARCH_PATH/bin
+    _foamAddPath    $MPI_ARCH_PATH/bin
-    _foamAddLib  $MPI_ARCH_PATH/lib
+    _foamAddLib     $MPI_ARCH_PATH/lib
+    _foamAddManPath $MPI_ARCH_PATH/man
 
     export FOAM_MPI_LIBBIN=$FOAM_LIBBIN/$mpi_version
     unset mpi_version
     ;;
 
 MPICH)
-    mpi_version=mpich-1.2.4
+    mpi_version=mpich2-1.1.1p1
     export MPI_HOME=$WM_THIRD_PARTY_DIR/$mpi_version
     export MPI_ARCH_PATH=$MPI_HOME/platforms/$WM_OPTIONS
-    export MPICH_ROOT=$MPI_ARCH_PATH
 
-    _foamAddPath $MPI_ARCH_PATH/bin
+    _foamAddPath    $MPI_ARCH_PATH/bin
-    _foamAddLib  $MPI_ARCH_PATH/lib
+    _foamAddLib     $MPI_ARCH_PATH/lib
+    _foamAddManPath $MPI_ARCH_PATH/share/man
 
     export FOAM_MPI_LIBBIN=$FOAM_LIBBIN/$mpi_version
     unset mpi_version
@@ -181,7 +194,6 @@ MPICH)
 MPICH-GM)
     export MPI_ARCH_PATH=/opt/mpi
     export MPICH_PATH=$MPI_ARCH_PATH
-    export MPICH_ROOT=$MPI_ARCH_PATH
     export GM_LIB_PATH=/opt/gm/lib64
 
     _foamAddPath $MPI_ARCH_PATH/bin
@@ -194,7 +206,6 @@ MPICH-GM)
 HPMPI)
     export MPI_HOME=/opt/hpmpi
     export MPI_ARCH_PATH=$MPI_HOME
-    export MPICH_ROOT=$MPI_ARCH_PATH
 
     _foamAddPath $MPI_ARCH_PATH/bin
 

src/OpenFOAM/db/IOstreams/Pstreams/OPstream.C

@@ -85,7 +85,6 @@ inline void Foam::OPstream::writeToBuffer
 }
 
 
-
 // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
 
 Foam::OPstream::OPstream

src/OpenFOAM/db/IOstreams/Pstreams/Pstream.C

@@ -217,7 +217,7 @@ bool Foam::Pstream::parRun_(false);
 Foam::List<int> Foam::Pstream::procIDs_(1, 0);
 
 // Standard transfer message type
-const int Foam::Pstream::msgType_(1);
+int Foam::Pstream::msgType_(1);
 
 // Linear communication schedule
 Foam::List<Foam::Pstream::commsStruct> Foam::Pstream::linearCommunication_(0);

src/OpenFOAM/db/IOstreams/Pstreams/Pstream.H

@@ -163,7 +163,7 @@ private:
         static bool parRun_;
 
         static List<int> procIDs_;
-        static const int msgType_;
+        static int msgType_;
 
         static List<commsStruct> linearCommunication_;
         static List<commsStruct> treeCommunication_;
@@ -337,11 +337,12 @@ public:
         }
 
         //- Message tag of standard messages
-        static int msgType()
+        static int& msgType()
         {
             return msgType_;
         }
 
+
             //- Get the communications type of the stream
             commsTypes commsType() const
             {

src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryIO.C

@@ -81,7 +81,12 @@ bool Foam::primitiveEntry::expandVariable
     word varName = w(1, w.size()-1);
 
     // lookup the variable name in the given dictionary....
-    const entry* ePtr = dict.lookupEntryPtr(varName, true, true);
+    // Note: allow wildcards to match? For now disabled since following
+    // would expand internalField to wildcard match and not expected
+    // internalField:
+    //      internalField XXX;
+    //      boundaryField { ".*" {YYY;} movingWall {value $internalField;}
+    const entry* ePtr = dict.lookupEntryPtr(varName, true, false);
 
     // ...if defined insert its tokens into this
     if (ePtr != NULL)

src/OpenFOAM/fields/Fields/symmTransformField/symmTransformField.H

@@ -38,6 +38,7 @@ SourceFiles
 
 #include "symmTransform.H"
 #include "symmTensorField.H"
+#include "sphericalTensor.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/fields/pointPatchFields/basic/basicSymmetry/basicSymmetryPointPatchField.C

@@ -26,6 +26,7 @@ License
 
 #include "basicSymmetryPointPatchField.H"
 #include "transformField.H"
+#include "symmTransformField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -89,7 +90,12 @@ void basicSymmetryPointPatchField<Type>::evaluate(const Pstream::commsTypes)
     const vectorField& nHat = this->patch().pointNormals();
 
     tmp<Field<Type> > tvalues =
-        transform(I - nHat*nHat, this->patchInternalField());
+    (
+        (
+            this->patchInternalField()
+          + transform(I - 2.0*sqr(nHat), this->patchInternalField())
+        )/2.0
+    );
 
     // Get internal field to insert values into
     Field<Type>& iF = const_cast<Field<Type>&>(this->internalField());

src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C

@@ -445,6 +445,19 @@ Foam::coupledPolyPatch::coupledPolyPatch
 {}
 
 
+Foam::coupledPolyPatch::coupledPolyPatch
+(
+    const coupledPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    polyPatch(pp, bm, index, mapAddressing, newStart)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::coupledPolyPatch::~coupledPolyPatch()

src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H

@@ -221,6 +221,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        coupledPolyPatch
+        (
+            const coupledPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
 
     // Destructor
 

src/OpenFOAM/meshes/polyMesh/polyPatches/basic/generic/genericPolyPatch.C

@@ -94,6 +94,21 @@ Foam::genericPolyPatch::genericPolyPatch
 {}
 
 
+Foam::genericPolyPatch::genericPolyPatch
+(
+    const genericPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    polyPatch(pp, bm, index, mapAddressing, newStart),
+    actualTypeName_(pp.actualTypeName_),
+    dict_(pp.dict_)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::genericPolyPatch::~genericPolyPatch()

src/OpenFOAM/meshes/polyMesh/polyPatches/basic/generic/genericPolyPatch.H

@@ -106,6 +106,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        genericPolyPatch
+        (
+            const genericPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -128,6 +138,22 @@ public:
             );
         }
 
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new genericPolyPatch(*this, bm, index, mapAddressing, newStart)
+            );
+        }
+
 
     // Destructor
 

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C

@@ -826,6 +826,26 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
 {}
 
 
+Foam::cyclicPolyPatch::cyclicPolyPatch
+(
+    const cyclicPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    coupledPolyPatch(pp, bm, index, mapAddressing, newStart),
+    coupledPointsPtr_(NULL),
+    coupledEdgesPtr_(NULL),
+    featureCos_(pp.featureCos_),
+    transform_(pp.transform_),
+    rotationAxis_(pp.rotationAxis_),
+    rotationCentre_(pp.rotationCentre_),
+    separationVector_(pp.separationVector_)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::cyclicPolyPatch::~cyclicPolyPatch()

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H

@@ -234,6 +234,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        cyclicPolyPatch
+        (
+            const cyclicPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -256,6 +266,22 @@ public:
             );
         }
 
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new cyclicPolyPatch(*this, bm, index, mapAddressing, newStart)
+            );
+        }
+
 
     // Destructor
 

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/empty/emptyPolyPatch.C

@@ -87,4 +87,17 @@ Foam::emptyPolyPatch::emptyPolyPatch
 {}
 
 
+Foam::emptyPolyPatch::emptyPolyPatch
+(
+    const emptyPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    polyPatch(pp, bm, index, mapAddressing, newStart)
+{}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/empty/emptyPolyPatch.H

@@ -92,6 +92,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        emptyPolyPatch
+        (
+            const emptyPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -113,6 +123,22 @@ public:
                 new emptyPolyPatch(*this, bm, index, newSize, newStart)
             );
         }
+
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new emptyPolyPatch(*this, bm, index, mapAddressing, newStart)
+            );
+        }
 };
 
 

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C

@@ -124,6 +124,26 @@ Foam::processorPolyPatch::processorPolyPatch
 {}
 
 
+Foam::processorPolyPatch::processorPolyPatch
+(
+    const processorPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+     const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    coupledPolyPatch(pp, bm, index, mapAddressing, newStart),
+    myProcNo_(pp.myProcNo_),
+    neighbProcNo_(pp.neighbProcNo_),
+    neighbFaceCentres_(),
+    neighbFaceAreas_(),
+    neighbFaceCellCentres_(),
+    neighbPointsPtr_(NULL),
+    neighbEdgesPtr_(NULL)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::processorPolyPatch::~processorPolyPatch()

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H

@@ -158,6 +158,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        processorPolyPatch
+        (
+            const processorPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -178,7 +188,7 @@ public:
             (
                 new processorPolyPatch
                 (
-                    refCast<const processorPolyPatch>(*this),
+                    *this,
                     bm,
                     index,
                     newSize,
@@ -187,6 +197,29 @@ public:
             );
         }
 
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new processorPolyPatch
+                (
+                    *this,
+                    bm,
+                    index,
+                    mapAddressing,
+                    newStart
+                )
+            );
+        }
+
 
     // Destructor
 

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/symmetry/symmetryPolyPatch.C

@@ -87,4 +87,17 @@ Foam::symmetryPolyPatch::symmetryPolyPatch
 {}
 
 
+Foam::symmetryPolyPatch::symmetryPolyPatch
+(
+    const symmetryPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    polyPatch(pp, bm, index, mapAddressing, newStart)
+{}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/symmetry/symmetryPolyPatch.H

@@ -92,6 +92,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        symmetryPolyPatch
+        (
+            const symmetryPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -113,6 +123,22 @@ public:
                 new symmetryPolyPatch(*this, bm, index, newSize, newStart)
             );
         }
+
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new symmetryPolyPatch(*this, bm, index, mapAddressing, newStart)
+            );
+        }
 };
 
 

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C

@@ -157,4 +157,19 @@ Foam::wedgePolyPatch::wedgePolyPatch
 }
 
 
+Foam::wedgePolyPatch::wedgePolyPatch
+(
+    const wedgePolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    polyPatch(pp, bm, index, mapAddressing, newStart)
+{
+    initTransforms();
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.H

@@ -113,6 +113,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        wedgePolyPatch
+        (
+            const wedgePolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -135,6 +145,22 @@ public:
             );
         }
 
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new wedgePolyPatch(*this, bm, index, mapAddressing, newStart)
+            );
+        }
+
 
     // Member functions
 

src/OpenFOAM/meshes/polyMesh/polyPatches/derived/wall/wallPolyPatch.C

@@ -87,4 +87,17 @@ Foam::wallPolyPatch::wallPolyPatch
 {}
 
 
+Foam::wallPolyPatch::wallPolyPatch
+(
+    const wallPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    polyPatch(pp, bm, index, mapAddressing, newStart)
+{}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/polyPatches/derived/wall/wallPolyPatch.H

@@ -92,6 +92,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        wallPolyPatch
+        (
+            const wallPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -113,6 +123,22 @@ public:
                 new wallPolyPatch(*this, bm, index, newSize, newStart)
             );
         }
+
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new wallPolyPatch(*this, bm, index, mapAddressing, newStart)
+            );
+        }
 };
 
 

src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C

@@ -167,6 +167,33 @@ Foam::polyPatch::polyPatch
 {}
 
 
+Foam::polyPatch::polyPatch
+(
+    const polyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    patchIdentifier(pp, index),
+    primitivePatch
+    (
+        faceSubList
+        (
+            bm.mesh().faces(),
+            mapAddressing.size(),
+            newStart
+        ),
+        bm.mesh().points()
+    ),
+    start_(newStart),
+    boundaryMesh_(bm),
+    faceCellsPtr_(NULL),
+    mePtr_(NULL)
+{}
+
+
 Foam::polyPatch::polyPatch(const polyPatch& p)
 :
     patchIdentifier(p),

src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H

@@ -199,6 +199,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        polyPatch
+        (
+            const polyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct as copy
         polyPatch(const polyPatch&);
 
@@ -224,6 +234,22 @@ public:
             );
         }
 
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new polyPatch(*this, bm, index, mapAddressing, newStart)
+            );
+        }
+
 
     // Selectors
 

src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C

@@ -134,16 +134,24 @@ void Foam::faceZone::calcCellLayers() const
 
         forAll (mf, faceI)
         {
+            label ownCellI = own[mf[faceI]];
+            label neiCellI =
+            (
+                zoneMesh().mesh().isInternalFace(mf[faceI])
+              ? nei[mf[faceI]]
+              : -1
+            );
+
             if (!faceFlip[faceI])
             {
                 // Face is oriented correctly, no flip needed
-                mc[faceI] = nei[mf[faceI]];
+                mc[faceI] = neiCellI;
-                sc[faceI] = own[mf[faceI]];
+                sc[faceI] = ownCellI;
             }
             else
             {
-                mc[faceI] = own[mf[faceI]];
+                mc[faceI] = ownCellI;
-                sc[faceI] = nei[mf[faceI]];
+                sc[faceI] = neiCellI;
             }
         }
         //Info << "masterCells: " << mc << endl;

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -2509,138 +2509,156 @@ void Foam::autoLayerDriver::addLayers
 (
     const layerParameters& layerParams,
     const dictionary& motionDict,
+    const labelList& patchIDs,
     const label nAllowableErrors,
-    motionSmoother& meshMover,
     decompositionMethod& decomposer,
     fvMeshDistribute& distributor
 )
 {
     fvMesh& mesh = meshRefiner_.mesh();
-    const indirectPrimitivePatch& pp = meshMover.patch();
-    const labelList& meshPoints = pp.meshPoints();
 
-    // Precalculate mesh edge labels for patch edges
+    autoPtr<indirectPrimitivePatch> pp
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    (
-
+        meshRefinement::makePatch
-    labelList meshEdges(pp.nEdges());
-    forAll(pp.edges(), edgeI)
-    {
-        const edge& ppEdge = pp.edges()[edgeI];
-        label v0 = meshPoints[ppEdge[0]];
-        label v1 = meshPoints[ppEdge[1]];
-        meshEdges[edgeI] = meshTools::findEdge
         (
-            mesh.edges(),
+            mesh,
-            mesh.pointEdges()[v0],
+            patchIDs
-            v0,
+        )
-            v1
-        );
-    }
-
-    // Displacement for all pp.localPoints.
-    vectorField patchDisp(pp.nPoints(), vector::one);
-
-    // Number of layers for all pp.localPoints. Note: only valid if
-    // extrudeStatus = EXTRUDE.
-    labelList patchNLayers(pp.nPoints(), 0);
-
-    // Whether to add edge for all pp.localPoints.
-    List<extrudeMode> extrudeStatus(pp.nPoints(), EXTRUDE);
-
-
-    // Get number of layer per point from number of layers per patch
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    setNumLayers
-    (
-        layerParams.numLayers(),    // per patch the num layers
-        meshMover.adaptPatchIDs(),  // patches that are being moved
-        pp,                         // indirectpatch for all faces moving
-
-        patchDisp,
-        patchNLayers,
-        extrudeStatus
     );
 
-    // Disable extrusion on non-manifold points
+    // Construct iterative mesh mover.
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    Info<< "Constructing mesh displacer ..." << endl;
 
-    handleNonManifolds
+    autoPtr<motionSmoother> meshMover
     (
-        pp,
+        new motionSmoother
-        meshEdges,
+        (
-
+            mesh,
-        patchDisp,
+            pp(),
-        patchNLayers,
+            patchIDs,
-        extrudeStatus
+            meshRefinement::makeDisplacementField
+            (
+                pointMesh::New(mesh),
+                patchIDs
+            ),
+            motionDict
+        )
     );
 
-    // Disable extrusion on feature angles
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    handleFeatureAngle
-    (
-        pp,
-        meshEdges,
-        layerParams.featureAngle()*constant::math::pi/180.0,
-
-        patchDisp,
-        patchNLayers,
-        extrudeStatus
-    );
-
-    // Disable extrusion on warped faces
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Undistorted edge length
     const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
     const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
 
-    handleWarpedFaces
-    (
-        pp,
-        layerParams.maxFaceThicknessRatio(),
-        edge0Len,
-        cellLevel,
 
-        patchDisp,
+    // Point-wise extrusion data
-        patchNLayers,
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-        extrudeStatus
-    );
 
-    //// Disable extrusion on cells with multiple patch faces
+    // Displacement for all pp.localPoints.
-    //// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    vectorField patchDisp(pp().nPoints(), vector::one);
-    //
-    //handleMultiplePatchFaces
-    //(
-    //    pp,
-    //
-    //    patchDisp,
-    //    patchNLayers,
-    //    extrudeStatus
-    //);
 
+    // Number of layers for all pp.localPoints. Note: only valid if
+    // extrudeStatus = EXTRUDE.
+    labelList patchNLayers(pp().nPoints(), 0);
+
+    // Whether to add edge for all pp.localPoints.
+    List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
 
-    // Grow out region of non-extrusion
-    for (label i = 0; i < layerParams.nGrow(); i++)
     {
-        growNoExtrusion
+        // Get number of layer per point from number of layers per patch
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        setNumLayers
         (
-            pp,
+            layerParams.numLayers(),    // per patch the num layers
+            meshMover().adaptPatchIDs(),// patches that are being moved
+            pp,                         // indirectpatch for all faces moving
+
             patchDisp,
             patchNLayers,
             extrudeStatus
         );
+
+        // Precalculate mesh edge labels for patch edges
+        labelList meshEdges(pp().meshEdges(mesh.edges(), mesh.pointEdges()));
+
+        // Disable extrusion on non-manifold points
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        handleNonManifolds
+        (
+            pp,
+            meshEdges,
+
+            patchDisp,
+            patchNLayers,
+            extrudeStatus
+        );
+
+        // Disable extrusion on feature angles
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        handleFeatureAngle
+        (
+            pp,
+            meshEdges,
+            layerParams.featureAngle()*constant::math::pi/180.0,
+
+            patchDisp,
+            patchNLayers,
+            extrudeStatus
+        );
+
+        // Disable extrusion on warped faces
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        handleWarpedFaces
+        (
+            pp,
+            layerParams.maxFaceThicknessRatio(),
+            edge0Len,
+            cellLevel,
+
+            patchDisp,
+            patchNLayers,
+            extrudeStatus
+        );
+
+        //// Disable extrusion on cells with multiple patch faces
+        //// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        //
+        //handleMultiplePatchFaces
+        //(
+        //    pp,
+        //
+        //    patchDisp,
+        //    patchNLayers,
+        //    extrudeStatus
+        //);
+
+
+        // Grow out region of non-extrusion
+        for (label i = 0; i < layerParams.nGrow(); i++)
+        {
+            growNoExtrusion
+            (
+                pp,
+                patchDisp,
+                patchNLayers,
+                extrudeStatus
+            );
+        }
     }
 
+
     // Determine (wanted) point-wise layer thickness and expansion ratio
-    scalarField thickness(pp.nPoints());
+    scalarField thickness(pp().nPoints());
-    scalarField minThickness(pp.nPoints());
+    scalarField minThickness(pp().nPoints());
-    scalarField expansionRatio(pp.nPoints());
+    scalarField expansionRatio(pp().nPoints());
     calculateLayerThickness
     (
         pp,
-        meshMover.adaptPatchIDs(),
+        meshMover().adaptPatchIDs(),
         layerParams.expansionRatio(),
 
         layerParams.relativeSizes(),        // thickness relative to cellsize?
@@ -2663,11 +2681,11 @@ void Foam::autoLayerDriver::addLayers
 
         // Find maximum length of a patch name, for a nicer output
         label maxPatchNameLen = 0;
-        forAll(meshMover.adaptPatchIDs(), i)
+        forAll(meshMover().adaptPatchIDs(), i)
         {
-            label patchI = meshMover.adaptPatchIDs()[i];
+            label patchI = meshMover().adaptPatchIDs()[i];
             word patchName = patches[patchI].name();
-            maxPatchNameLen = max(maxPatchNameLen,label(patchName.size()));
+            maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
         }
 
         Info<< nl 
@@ -2678,9 +2696,9 @@ void Foam::autoLayerDriver::addLayers
             << setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
             << setw(0) << " -----    ------ --------- -------" << endl;
 
-        forAll(meshMover.adaptPatchIDs(), i)
+        forAll(meshMover().adaptPatchIDs(), i)
         {
-            label patchI = meshMover.adaptPatchIDs()[i];
+            label patchI = meshMover().adaptPatchIDs()[i];
 
             const labelList& meshPoints = patches[patchI].meshPoints();
 
@@ -2691,7 +2709,7 @@ void Foam::autoLayerDriver::addLayers
 
             forAll(meshPoints, patchPointI)
             {
-                label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
+                label ppPointI = pp().meshPointMap()[meshPoints[patchPointI]];
 
                 //maxThickness = max(maxThickness, thickness[ppPointI]);
                 //minThickness = min(minThickness, thickness[ppPointI]);
@@ -2761,7 +2779,7 @@ void Foam::autoLayerDriver::addLayers
             IOobject::NO_WRITE,
             false
         ),
-        meshMover.pMesh(),
+        meshMover().pMesh(),
         dimensionedScalar("pointMedialDist", dimless, 0.0)
     );
 
@@ -2776,7 +2794,7 @@ void Foam::autoLayerDriver::addLayers
             IOobject::NO_WRITE,
             false
         ),
-        meshMover.pMesh(),
+        meshMover().pMesh(),
         dimensionedVector("dispVec", dimless, vector::zero)
     );
 
@@ -2791,7 +2809,7 @@ void Foam::autoLayerDriver::addLayers
             IOobject::NO_WRITE,
             false
         ),
-        meshMover.pMesh(),
+        meshMover().pMesh(),
         dimensionedScalar("medialRatio", dimless, 0.0)
     );
 
@@ -2871,7 +2889,7 @@ void Foam::autoLayerDriver::addLayers
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
         {
-            pointField oldPatchPos(pp.localPoints());
+            pointField oldPatchPos(pp().localPoints());
 
             //// Laplacian displacement shrinking.
             //shrinkMeshDistance
@@ -2886,7 +2904,7 @@ void Foam::autoLayerDriver::addLayers
             // Medial axis based shrinking
             shrinkMeshMedialDistance
             (
-                meshMover,
+                meshMover(),
                 meshQualityDict,
 
                 layerParams.nSmoothThickness(),
@@ -2908,7 +2926,7 @@ void Foam::autoLayerDriver::addLayers
             );
 
             // Update patchDisp (since not all might have been honoured)
-            patchDisp = oldPatchPos - pp.localPoints();
+            patchDisp = oldPatchPos - pp().localPoints();
         }
 
         // Truncate displacements that are too small (this will do internal
@@ -2916,7 +2934,7 @@ void Foam::autoLayerDriver::addLayers
         faceSet dummySet(mesh, "wrongPatchFaces", 0);
         truncateDisplacement
         (
-            meshMover,
+            meshMover(),
             minThickness,
             dummySet,
             patchDisp,
@@ -2931,7 +2949,7 @@ void Foam::autoLayerDriver::addLayers
             dumpDisplacement
             (
                 mesh.time().path()/"layer",
-                pp,
+                pp(),
                 patchDisp,
                 extrudeStatus
             );
@@ -2955,11 +2973,11 @@ void Foam::autoLayerDriver::addLayers
         // Determine per point/per face number of layers to extrude. Also
         // handles the slow termination of layers when going switching layers
 
-        labelList nPatchPointLayers(pp.nPoints(),-1);
+        labelList nPatchPointLayers(pp().nPoints(),-1);
-        labelList nPatchFaceLayers(pp.localFaces().size(),-1);
+        labelList nPatchFaceLayers(pp().localFaces().size(),-1);
         setupLayerInfoTruncation
         (
-            meshMover,
+            meshMover(),
             patchNLayers,
             extrudeStatus,
             layerParams.nBufferCellsNoExtrude(),
@@ -2998,7 +3016,7 @@ void Foam::autoLayerDriver::addLayers
         addLayer.setRefinement
         (
             invExpansionRatio,
-            pp,
+            pp(),
             labelList(0),       // exposed patchIDs, not used for adding layers
             nPatchFaceLayers,   // layers per face
             nPatchPointLayers,  // layers per point
@@ -3050,8 +3068,8 @@ void Foam::autoLayerDriver::addLayers
         addLayer.updateMesh
         (
             map,
-            identity(pp.size()),
+            identity(pp().size()),
-            identity(pp.nPoints())
+            identity(pp().nPoints())
         );
 
         // Collect layer faces and cells for outside loop.
@@ -3098,7 +3116,7 @@ void Foam::autoLayerDriver::addLayers
         (
             addLayer,
             meshQualityDict,
-            pp,
+            pp(),
             newMesh,
 
             patchDisp,
@@ -3107,7 +3125,7 @@ void Foam::autoLayerDriver::addLayers
         );
 
         Info<< "Extruding " << countExtrusion(pp, extrudeStatus)
-            << " out of " << returnReduce(pp.size(), sumOp<label>())
+            << " out of " << returnReduce(pp().size(), sumOp<label>())
             << " faces. Removed extrusion at " << nTotChanged << " faces."
             << endl;
 
@@ -3117,8 +3135,8 @@ void Foam::autoLayerDriver::addLayers
         }
 
         // Reset mesh points and start again
-        meshMover.movePoints(oldPoints);
+        meshMover().movePoints(oldPoints);
-        meshMover.correct();
+        meshMover().correct();
 
         Info<< endl;
     }
@@ -3173,6 +3191,7 @@ void Foam::autoLayerDriver::addLayers
         (
             false,
             false,
+            scalarField(mesh.nCells(), 1.0),
             decomposer,
             distributor
         );
@@ -3211,7 +3230,7 @@ void Foam::autoLayerDriver::doLayers
     fvMeshDistribute& distributor
 )
 {
-    fvMesh& mesh = meshRefiner_.mesh();
+    const fvMesh& mesh = meshRefiner_.mesh();
 
     Info<< nl
         << "Shrinking and layer addition phase" << nl
@@ -3245,59 +3264,80 @@ void Foam::autoLayerDriver::doLayers
     }
     else
     {
-        autoPtr<indirectPrimitivePatch> ppPtr
+        // Check that outside of mesh is not multiply connected.
+        checkMeshManifold();
+
+        // Check initial mesh
+        Info<< "Checking initial mesh ..." << endl;
+        labelHashSet wrongFaces(mesh.nFaces()/100);
+        motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
+        const label nInitErrors = returnReduce
         (
-            meshRefinement::makePatch
+            wrongFaces.size(),
-            (
+            sumOp<label>()
-                mesh,
-                patchIDs
-            )
         );
-        indirectPrimitivePatch& pp = ppPtr();
 
-        // Construct iterative mesh mover.
+        Info<< "Detected " << nInitErrors << " illegal faces"
-        Info<< "Constructing mesh displacer ..." << endl;
+            << " (concave, zero area or negative cell pyramid volume)"
+            << endl;
 
+
+        // Balance
+        if (Pstream::parRun())
         {
-            const pointMesh& pMesh = pointMesh::New(mesh);
+            Info<< nl
-
+                << "Doing initial balancing" << nl
-            motionSmoother meshMover
+                << "-----------------------" << nl
-            (
-                mesh,
-                pp,
-                patchIDs,
-                meshRefinement::makeDisplacementField(pMesh, patchIDs),
-                motionDict
-            );
-
-            // Check that outside of mesh is not multiply connected.
-            checkMeshManifold();
-
-            // Check initial mesh
-            Info<< "Checking initial mesh ..." << endl;
-            labelHashSet wrongFaces(mesh.nFaces()/100);
-            motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
-            const label nInitErrors = returnReduce
-            (
-                wrongFaces.size(),
-                sumOp<label>()
-            );
-
-            Info<< "Detected " << nInitErrors << " illegal faces"
-                << " (concave, zero area or negative cell pyramid volume)"
                 << endl;
-
+        
-            // Do all topo changes
+            scalarField cellWeights(mesh.nCells(), 1);
-            addLayers
+            forAll(numLayers, patchI)
+            {
+                if (numLayers[patchI] > 0)
+                {
+                    const polyPatch& pp = mesh.boundaryMesh()[patchI];
+                    forAll(pp.faceCells(), i)
+                    {
+                        cellWeights[pp.faceCells()[i]] += numLayers[patchI];
+                    }
+                }
+            }
+        
+            // Balance mesh (and meshRefinement). No restriction on face zones
+            // and baffles.
+            autoPtr<mapDistributePolyMesh> map = meshRefiner_.balance
             (
-                layerParams,
+                false,
-                motionDict,
+                false,
-                nInitErrors,
+                cellWeights,
-                meshMover,
                 decomposer,
                 distributor
             );
+        
+            //{
+            //    globalIndex globalCells(mesh.nCells());
+            //
+            //    Info<< "** Distribution after balancing:" << endl;
+            //    for (label procI = 0; procI < Pstream::nProcs(); procI++)
+            //    {
+            //        Info<< "    " << procI << '\t'
+            //            << globalCells.localSize(procI) << endl;
+            //    }
+            //    Info<< endl;
+            //}
         }
+
+
+        // Do all topo changes
+        addLayers
+        (
+            layerParams,
+            motionDict,
+            patchIDs,
+            nInitErrors,
+            decomposer,
+            distributor
+        );
     }
 }
 

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H

@@ -523,8 +523,8 @@ public:
             (
                 const layerParameters& layerParams,
                 const dictionary& motionDict,
+                const labelList& patchIDs,
                 const label nAllowableErrors,
-                motionSmoother& meshMover,
                 decompositionMethod& decomposer,
                 fvMeshDistribute& distributor
             );

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C

@@ -35,7 +35,6 @@ License
 #include "refinementSurfaces.H"
 #include "shellSurfaces.H"
 #include "mapDistributePolyMesh.H"
-#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -210,7 +209,8 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
                 "feature refinement iteration " + name(iter),
                 decomposer_,
                 distributor_,
-                cellsToRefine
+                cellsToRefine,
+                refineParams.maxLoadUnbalance()
             );
         }
     }
@@ -310,7 +310,8 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
             "surface refinement iteration " + name(iter),
             decomposer_,
             distributor_,
-            cellsToRefine
+            cellsToRefine,
+            refineParams.maxLoadUnbalance()
         );
     }
     return iter;
@@ -492,7 +493,8 @@ Foam::label Foam::autoRefineDriver::shellRefine
             "shell refinement iteration " + name(iter),
             decomposer_,
             distributor_,
-            cellsToRefine
+            cellsToRefine,
+            refineParams.maxLoadUnbalance()
         );
     }
     meshRefiner_.userFaceData().clear();
@@ -776,11 +778,14 @@ void Foam::autoRefineDriver::doRefine
             const_cast<Time&>(mesh.time())++;
         }
 
-        // Do final balancing. Keep zoned faces on one processor.
+        // Do final balancing. Keep zoned faces on one processor since the
+        // snap phase will convert them to baffles and this only works for
+        // internal faces.
         meshRefiner_.balance
         (
             true,
             false,
+            scalarField(mesh.nCells(), 1), // dummy weights
             decomposer_,
             distributor_
         );

src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C

@@ -43,7 +43,8 @@ Foam::refinementParameters::refinementParameters
     minRefineCells_(readLabel(dict.lookup("minimumRefine"))),
     curvature_(readScalar(dict.lookup("curvature"))),
     nBufferLayers_(readLabel(dict.lookup("nBufferLayers"))),
-    keepPoints_(dict.lookup("keepPoints"))
+    keepPoints_(dict.lookup("keepPoints")),
+    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0))
 {}
 
 
@@ -53,7 +54,8 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
     maxLocalCells_(readLabel(dict.lookup("maxLocalCells"))),
     minRefineCells_(readLabel(dict.lookup("minRefinementCells"))),
     nBufferLayers_(readLabel(dict.lookup("nCellsBetweenLevels"))),
-    keepPoints_(pointField(1, dict.lookup("locationInMesh")))
+    keepPoints_(pointField(1, dict.lookup("locationInMesh"))),
+    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0))
 {
     scalar featAngle(readScalar(dict.lookup("resolveFeatureAngle")));
 

src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H

@@ -73,6 +73,9 @@ class refinementParameters
         //- Areas to keep
         const pointField keepPoints_;
 
+        //- Allowed load unbalance
+        scalar maxLoadUnbalance_;
+
 
     // Private Member Functions
 
@@ -134,6 +137,12 @@ public:
                 return keepPoints_;
             }
 
+            //- Allowed load unbalance
+            scalar maxLoadUnbalance() const
+            {
+                return maxLoadUnbalance_;
+            }
+
 
         // Other
 

src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C

@@ -552,13 +552,16 @@ void Foam::meshRefinement::calcLocalRegions
     const globalIndex& globalCells,
     const labelList& globalRegion,
     const Map<label>& coupledRegionToMaster,
+    const scalarField& cellWeights,
 
     Map<label>& globalToLocalRegion,
-    pointField& localPoints
+    pointField& localPoints,
+    scalarField& localWeights
 ) const
 {
     globalToLocalRegion.resize(globalRegion.size());
     DynamicList<point> localCc(globalRegion.size()/2);
+    DynamicList<scalar> localWts(globalRegion.size()/2);
 
     forAll(globalRegion, cellI)
     {
@@ -573,6 +576,7 @@ void Foam::meshRefinement::calcLocalRegions
                 // I am master. Allocate region for me.
                 globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
                 localCc.append(mesh_.cellCentres()[cellI]);
+                localWts.append(cellWeights[cellI]);
             }
         }
         else
@@ -581,11 +585,13 @@ void Foam::meshRefinement::calcLocalRegions
             if (globalToLocalRegion.insert(globalRegion[cellI], localCc.size()))
             {
                 localCc.append(mesh_.cellCentres()[cellI]);
+                localWts.append(cellWeights[cellI]);
             }
         }
     }
 
     localPoints.transfer(localCc);
+    localWeights.transfer(localWts);
 
     if (localPoints.size() != globalToLocalRegion.size())
     {
@@ -924,6 +930,7 @@ Foam::label Foam::meshRefinement::countHits() const
 // Determine distribution to move connected regions onto one processor.
 Foam::labelList Foam::meshRefinement::decomposeCombineRegions
 (
+    const scalarField& cellWeights,
     const boolList& blockedFace,
     const List<labelPair>& explicitConnections,
     decompositionMethod& decomposer
@@ -965,14 +972,17 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
 
     Map<label> globalToLocalRegion;
     pointField localPoints;
+    scalarField localWeights;
     calcLocalRegions
     (
         globalCells,
         globalRegion,
         coupledRegionToMaster,
+        cellWeights,
 
         globalToLocalRegion,
-        localPoints
+        localPoints,
+        localWeights
     );
 
 
@@ -984,7 +994,7 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
 
     if (isA<geomDecomp>(decomposer))
     {
-        regionDistribution = decomposer.decompose(localPoints);
+        regionDistribution = decomposer.decompose(localPoints, localWeights);
     }
     else
     {
@@ -998,7 +1008,12 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
             regionRegions
         );
 
-        regionDistribution = decomposer.decompose(regionRegions, localPoints);
+        regionDistribution = decomposer.decompose
+        (
+            regionRegions,
+            localPoints,
+            localWeights
+        );
     }
 
 
@@ -1058,6 +1073,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
 (
     const bool keepZoneFaces,
     const bool keepBaffles,
+    const scalarField& cellWeights,
     decompositionMethod& decomposer,
     fvMeshDistribute& distributor
 )
@@ -1118,6 +1134,18 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
                         }
                     }
                 }
+
+
+                // If the faceZones are not synchronised the blockedFace
+                // might not be synchronised. If you are sure the faceZones
+                // are synchronised remove below check.
+                syncTools::syncFaceList
+                (
+                    mesh_,
+                    blockedFace,
+                    andEqOp<bool>(),    // combine operator
+                    false               // separation
+                );
             }
             reduce(nUnblocked, sumOp<label>());
 
@@ -1151,6 +1179,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
 
                 distribution = decomposeCombineRegions
                 (
+                    cellWeights,
                     blockedFace,
                     couples,
                     decomposer
@@ -1170,13 +1199,21 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
             else
             {
                 // Normal decomposition
-                distribution = decomposer.decompose(mesh_.cellCentres());
+                distribution = decomposer.decompose
+                (
+                    mesh_.cellCentres(),
+                    cellWeights
+                );
             }
         }
         else
         {
             // Normal decomposition
-            distribution = decomposer.decompose(mesh_.cellCentres());
+            distribution = decomposer.decompose
+            (
+                mesh_.cellCentres(),
+                cellWeights
+            );
         }
 
         if (debug)

src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H

@@ -214,9 +214,11 @@ private:
                 const globalIndex& globalCells,
                 const labelList& globalRegion,
                 const Map<label>& coupledRegionToMaster,
+                const scalarField& cellWeights,
 
                 Map<label>& globalToLocalRegion,
-                pointField& localPoints
+                pointField& localPoints,
+                scalarField& localWeights
             ) const;
 
             //- Convert region into global index.
@@ -574,6 +576,7 @@ public:
             //  (e.g. to keep baffles together)
             labelList decomposeCombineRegions
             (
+                const scalarField& cellWeights,
                 const boolList& blockedFace,
                 const List<labelPair>& explicitConnections,
                 decompositionMethod&
@@ -587,6 +590,7 @@ public:
             (
                 const bool keepZoneFaces,
                 const bool keepBaffles,
+                const scalarField& cellWeights,
                 decompositionMethod& decomposer,
                 fvMeshDistribute& distributor
             );
@@ -648,7 +652,8 @@ public:
                 const string& msg,
                 decompositionMethod& decomposer,
                 fvMeshDistribute& distributor,
-                const labelList& cellsToRefine
+                const labelList& cellsToRefine,
+                const scalar maxLoadUnbalance
             );
 
 

src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C

@@ -37,6 +37,7 @@ License
 #include "mapDistributePolyMesh.H"
 #include "featureEdgeMesh.H"
 #include "Cloud.H"
+//#include "globalIndex.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -1244,40 +1245,122 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::refine
 }
 
 
-// Do refinement of consistent set of cells followed by truncation and
+//// Do refinement of consistent set of cells followed by truncation and
-// load balancing.
+//// load balancing.
+//Foam::autoPtr<Foam::mapDistributePolyMesh>
+//Foam::meshRefinement::refineAndBalance
+//(
+//    const string& msg,
+//    decompositionMethod& decomposer,
+//    fvMeshDistribute& distributor,
+//    const labelList& cellsToRefine
+//)
+//{
+//    // Do all refinement
+//    refine(cellsToRefine);
+//
+//    if (debug)
+//    {
+//        Pout<< "Writing refined but unbalanced " << msg
+//            << " mesh to time " << timeName() << endl;
+//        write
+//        (
+//            debug,
+//            mesh_.time().path()
+//           /timeName()
+//        );
+//        Pout<< "Dumped debug data in = "
+//            << mesh_.time().cpuTimeIncrement() << " s" << endl;
+//
+//        // test all is still synced across proc patches
+//        checkData();
+//    }
+//
+//    Info<< "Refined mesh in = "
+//        << mesh_.time().cpuTimeIncrement() << " s" << endl;
+//    printMeshInfo(debug, "After refinement " + msg);
+//
+//
+//    // Load balancing
+//    // ~~~~~~~~~~~~~~
+//
+//    autoPtr<mapDistributePolyMesh> distMap;
+//
+//    if (Pstream::nProcs() > 1)
+//    {
+//        scalarField cellWeights(mesh_.nCells(), 1);
+//
+//        distMap = balance
+//        (
+//            false,  //keepZoneFaces
+//            false,  //keepBaffles
+//            cellWeights,
+//            decomposer,
+//            distributor
+//        );
+//
+//        Info<< "Balanced mesh in = "
+//            << mesh_.time().cpuTimeIncrement() << " s" << endl;
+//
+//        printMeshInfo(debug, "After balancing " + msg);
+//
+//
+//        if (debug)
+//        {
+//            Pout<< "Writing balanced " << msg
+//                << " mesh to time " << timeName() << endl;
+//            write
+//            (
+//                debug,
+//                mesh_.time().path()/timeName()
+//            );
+//            Pout<< "Dumped debug data in = "
+//                << mesh_.time().cpuTimeIncrement() << " s" << endl;
+//
+//            // test all is still synced across proc patches
+//            checkData();
+//        }
+//    }
+//
+//    return distMap;
+//}
+
+
+// Do load balancing followed by refinement of consistent set of cells.
 Foam::autoPtr<Foam::mapDistributePolyMesh>
 Foam::meshRefinement::refineAndBalance
 (
     const string& msg,
     decompositionMethod& decomposer,
     fvMeshDistribute& distributor,
-    const labelList& cellsToRefine
+    const labelList& initCellsToRefine,
+    const scalar maxLoadUnbalance
 )
 {
-    // Do all refinement
+    labelList cellsToRefine(initCellsToRefine);
-    refine(cellsToRefine);
 
-    if (debug)
+    //{
-    {
+    //    globalIndex globalCells(mesh_.nCells());
-        Pout<< "Writing refined but unbalanced " << msg
+    //
-            << " mesh to time " << timeName() << endl;
+    //    Info<< "** Distribution before balancing/refining:" << endl;
-        write
+    //    for (label procI = 0; procI < Pstream::nProcs(); procI++)
-        (
+    //    {
-            debug,
+    //        Info<< "    " << procI << '\t'
-            mesh_.time().path()
+    //            << globalCells.localSize(procI) << endl;
-           /timeName()
+    //    }
-        );
+    //    Info<< endl;
-        Pout<< "Dumped debug data in = "
+    //}
-            << mesh_.time().cpuTimeIncrement() << " s" << endl;
+    //{
-
+    //    globalIndex globalCells(cellsToRefine.size());
-        // test all is still synced across proc patches
+    //
-        checkData();
+    //    Info<< "** Cells to be refined:" << endl;
-    }
+    //    for (label procI = 0; procI < Pstream::nProcs(); procI++)
-
+    //    {
-    Info<< "Refined mesh in = "
+    //        Info<< "    " << procI << '\t'
-        << mesh_.time().cpuTimeIncrement() << " s" << endl;
+    //            << globalCells.localSize(procI) << endl;
-    printMeshInfo(debug, "After refinement " + msg);
+    //    }
+    //    Info<< endl;
+    //}
 
 
     // Load balancing
@@ -1287,20 +1370,62 @@ Foam::meshRefinement::refineAndBalance
 
     if (Pstream::nProcs() > 1)
     {
-        distMap = balance
+        // First check if we need to balance at all. Precalculate number of
+        // cells after refinement and see what maximum difference is.
+        scalar nNewCells = scalar(mesh_.nCells() + 7*cellsToRefine.size());
+        scalar nIdealNewCells =
+            returnReduce(nNewCells, sumOp<scalar>())
+          / Pstream::nProcs();
+        scalar unbalance = returnReduce
         (
-            false,  //keepZoneFaces
+            mag(1.0-nNewCells/nIdealNewCells),
-            false,  //keepBaffles
+            maxOp<scalar>()
-            decomposer,
-            distributor
         );
 
-        Info<< "Balanced mesh in = "
+        if (unbalance <= maxLoadUnbalance)
-            << mesh_.time().cpuTimeIncrement() << " s" << endl;
+        {
+            Info<< "Skipping balancing since max unbalance " << unbalance
+                << " in = "
+                << mesh_.time().cpuTimeIncrement() << " s" << endl;
+        }
+        else
+        {
+            scalarField cellWeights(mesh_.nCells(), 1);
+            forAll(cellsToRefine, i)
+            {
+                cellWeights[cellsToRefine[i]] += 7;
+            }
+
+            distMap = balance
+            (
+                false,  //keepZoneFaces
+                false,  //keepBaffles
+                cellWeights,
+                decomposer,
+                distributor
+            );
+
+            // Update cells to refine
+            distMap().distributeCellIndices(cellsToRefine);
+
+            Info<< "Balanced mesh in = "
+                << mesh_.time().cpuTimeIncrement() << " s" << endl;
+        }
+
+        //{
+        //    globalIndex globalCells(mesh_.nCells());
+        //
+        //    Info<< "** Distribution after balancing:" << endl;
+        //    for (label procI = 0; procI < Pstream::nProcs(); procI++)
+        //    {
+        //        Info<< "    " << procI << '\t'
+        //            << globalCells.localSize(procI) << endl;
+        //    }
+        //    Info<< endl;
+        //}
 
         printMeshInfo(debug, "After balancing " + msg);
 
-
         if (debug)
         {
             Pout<< "Writing balanced " << msg
@@ -1318,6 +1443,46 @@ Foam::meshRefinement::refineAndBalance
         }
     }
 
+
+    // Refinement
+    // ~~~~~~~~~~
+
+    refine(cellsToRefine);
+
+    if (debug)
+    {
+        Pout<< "Writing refined " << msg
+            << " mesh to time " << timeName() << endl;
+        write
+        (
+            debug,
+            mesh_.time().path()
+           /timeName()
+        );
+        Pout<< "Dumped debug data in = "
+            << mesh_.time().cpuTimeIncrement() << " s" << endl;
+
+        // test all is still synced across proc patches
+        checkData();
+    }
+
+    Info<< "Refined mesh in = "
+        << mesh_.time().cpuTimeIncrement() << " s" << endl;
+
+    //{
+    //    globalIndex globalCells(mesh_.nCells());
+    //
+    //    Info<< "** After refinement distribution:" << endl;
+    //    for (label procI = 0; procI < Pstream::nProcs(); procI++)
+    //    {
+    //        Info<< "    " << procI << '\t'
+    //            << globalCells.localSize(procI) << endl;
+    //    }
+    //    Info<< endl;
+    //}
+
+    printMeshInfo(debug, "After refinement " + msg);
+
     return distMap;
 }
 

src/finiteVolume/fields/fvPatchFields/basic/basicSymmetry/basicSymmetryFvPatchField.C

@@ -95,11 +95,11 @@ basicSymmetryFvPatchField<Type>::basicSymmetryFvPatchField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-// return gradient at boundary
 template<class Type>
 tmp<Field<Type> > basicSymmetryFvPatchField<Type>::snGrad() const
 {
     vectorField nHat = this->patch().nf();
+
     return
     (
         transform(I - 2.0*sqr(nHat), this->patchInternalField())
@@ -108,7 +108,6 @@ tmp<Field<Type> > basicSymmetryFvPatchField<Type>::snGrad() const
 }
 
 
-// Evaluate the field on the patch
 template<class Type>
 void basicSymmetryFvPatchField<Type>::evaluate(const Pstream::commsTypes)
 {
@@ -118,6 +117,7 @@ void basicSymmetryFvPatchField<Type>::evaluate(const Pstream::commsTypes)
     }
 
     vectorField nHat = this->patch().nf();
+
     Field<Type>::operator=
     (
         (
@@ -130,11 +130,11 @@ void basicSymmetryFvPatchField<Type>::evaluate(const Pstream::commsTypes)
 }
 
 
-// Return defining fields
 template<class Type>
 tmp<Field<Type> > basicSymmetryFvPatchField<Type>::snGradTransformDiag() const
 {
     vectorField nHat = this->patch().nf();
+
     vectorField diag(nHat.size());
 
     diag.replace(vector::X, mag(nHat.component(vector::X)));

src/finiteVolume/fields/fvPatchFields/derived/totalTemperature/totalTemperatureFvPatchScalarField.C

@@ -188,7 +188,7 @@ void Foam::totalTemperatureFvPatchScalarField::write(Ostream& os) const
     {
         os.writeKeyword("phi") << phiName_ << token::END_STATEMENT << nl;
     }
-    if (phiName_ != "psi")
+    if (psiName_ != "psi")
     {
         os.writeKeyword("psi") << psiName_ << token::END_STATEMENT << nl;
     }

src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C

@@ -329,7 +329,7 @@ surfaceDisplacementPointPatchVectorField
     surfacesDict_(dict.subDict("geometry")),
     projectMode_(projectModeNames_.read(dict.lookup("projectMode"))),
     projectDir_(dict.lookup("projectDirection")),
-    wedgePlane_(readLabel(dict.lookup("wedgePlane"))),
+    wedgePlane_(dict.lookupOrDefault(dict.lookup("wedgePlane"), -1)),
     frozenPointsZone_(dict.lookupOrDefault("frozenPointsZone", word::null))
 {
     if (velocity_.x() < 0 || velocity_.y() < 0 || velocity_.z() < 0)

src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C

@@ -326,7 +326,7 @@ surfaceSlipDisplacementPointPatchVectorField
     surfacesDict_(dict.subDict("geometry")),
     projectMode_(projectModeNames_.read(dict.lookup("projectMode"))),
     projectDir_(dict.lookup("projectDirection")),
-    wedgePlane_(readLabel(dict.lookup("wedgePlane"))),
+    wedgePlane_(dict.lookupOrDefault("wedgePlane", -1)),
     frozenPointsZone_(dict.lookupOrDefault("frozenPointsZone", word::null))
 {}
 

src/meshTools/Make/files

@@ -71,6 +71,9 @@ $(topoSets)/cellSet.C
 $(topoSets)/topoSet.C
 $(topoSets)/faceSet.C
 $(topoSets)/pointSet.C
+$(topoSets)/cellZoneSet.C
+$(topoSets)/faceZoneSet.C
+$(topoSets)/pointZoneSet.C
 
 sets/topoSetSource/topoSetSource.C
 
@@ -113,6 +116,18 @@ $(pointSources)/surfaceToPoint/surfaceToPoint.C
 $(pointSources)/zoneToPoint/zoneToPoint.C
 $(pointSources)/nearestToPoint/nearestToPoint.C
 
+faceZoneSources = sets/faceZoneSources
+$(faceZoneSources)/faceZoneToFaceZone/faceZoneToFaceZone.C
+$(faceZoneSources)/setsToFaceZone/setsToFaceZone.C
+$(faceZoneSources)/setToFaceZone/setToFaceZone.C
+
+cellZoneSources = sets/cellZoneSources
+$(cellZoneSources)/setToCellZone/setToCellZone.C
+
+pointZoneSources = sets/pointZoneSources
+$(pointZoneSources)/setToPointZone/setToPointZone.C
+
+
 surfaceSets/surfaceSets.C
 
 triSurface/orientedSurface/orientedSurface.C

src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C

@@ -50,6 +50,34 @@ namespace Foam
 
     const NamedEnum<directMappedPatchBase::sampleMode, 3>
         directMappedPatchBase::sampleModeNames_;
+
+
+    //- Private class for finding nearest
+    //  - point+local index
+    //  - sqr(distance)
+    //  - processor
+    typedef Tuple2<pointIndexHit, Tuple2<scalar, label> > nearInfo;
+
+    class nearestEqOp
+    {
+
+    public:
+
+        void operator()(nearInfo& x, const nearInfo& y) const
+        {
+            if (y.first().hit())
+            {
+                if (!x.first().hit())
+                {
+                    x = y;
+                }
+                else if (y.second().first() < x.second().first())
+                {
+                    x = y;
+                }
+            }
+        }
+    };
 }
 
 
@@ -70,7 +98,7 @@ void Foam::directMappedPatchBase::collectSamples
     labelListList globalFaces(Pstream::nProcs());
 
     globalFc[Pstream::myProcNo()] = patch_.faceCentres();
-    globalSamples[Pstream::myProcNo()] = globalFc[Pstream::myProcNo()]+offset_;
+    globalSamples[Pstream::myProcNo()] = globalFc[Pstream::myProcNo()]+offsets_;
     globalFaces[Pstream::myProcNo()] = identity(patch_.size());
 
     // Distribute to all processors
@@ -365,17 +393,17 @@ void Foam::directMappedPatchBase::calcMapping() const
 
     if
     (
-        offset_ == vector::zero
+        gAverage(mag(offsets_)) <= ROOTVSMALL
      && mode_ == NEARESTPATCHFACE
      && sampleRegion_ == patch_.boundaryMesh().mesh().name()
      && samplePatch_ == patch_.name()
     )
     {
         WarningIn("directMappedPatchBase::calcMapping() const")
-            << "Invalid offset " << offset_ << endl
+            << "Invalid offset " << offsets_ << endl
             << "Offset is the vector added to the patch face centres to"
             << " find the patch face supplying the data." << endl
-            << "Setting it to " << offset_
+            << "Setting it to " << offsets_
             << " on the same patch, on the same region"
             << " will find the faces themselves which does not make sense"
             << " for anything but testing." << endl
@@ -383,7 +411,7 @@ void Foam::directMappedPatchBase::calcMapping() const
             << "sampleRegion_:" << sampleRegion_ << endl
             << "mode_:" << sampleModeNames_[mode_] << endl
             << "samplePatch_:" << samplePatch_ << endl
-            << "offset_:" << offset_ << endl;
+            << "offsets_:" << offsets_ << endl;
     }
 
 
@@ -447,34 +475,38 @@ void Foam::directMappedPatchBase::calcMapping() const
     }
 
 
-    // Check that actual offset vector (sampleLocations - patchFc) is more or
+    //// Check that actual offset vector (sampleLocations - patchFc) is more or
-    // less constant.
+    //// less constant.
-    if (Pstream::master())
+    //if (Pstream::master())
-    {
+    //{
-        const scalarField magOffset(mag(sampleLocations - patchFc));
+    //    const scalarField magOffset(mag(sampleLocations - patchFc));
-        const scalar avgOffset(average(magOffset));
+    //    const scalar avgOffset(average(magOffset));
-
+    //
-        forAll(magOffset, sampleI)
+    //    forAll(magOffset, sampleI)
-        {
+    //    {
-            if (mag(magOffset[sampleI]-avgOffset) > max(SMALL, 0.001*avgOffset))
+    //        if
-            {
+    //        (
-                WarningIn("directMappedPatchBase::calcMapping() const")
+    //            mag(magOffset[sampleI]-avgOffset)
-                    << "The actual cell/face centres picked up using offset "
+    //          > max(SMALL, 0.001*avgOffset)
-                    << offset_ << " are not" << endl
+    //        )
-                    << "    on a single plane."
+    //        {
-                    << " This might give numerical problems." << endl
+    //            WarningIn("directMappedPatchBase::calcMapping() const")
-                    << "    At patchface " << patchFc[sampleI]
+    //                << "The actual cell/face centres picked up using offset "
-                    << " the sampled cell/face " << sampleLocations[sampleI]
+    //                << offsets_ << " are not" << endl
-                    << endl
+    //                << "    on a single plane."
-                    << "    is not on a plane " << avgOffset
+    //                << " This might give numerical problems." << endl
-                    << " offset from the patch." << endl
+    //                << "    At patchface " << patchFc[sampleI]
-                    << "    You might want to shift your plane offset."
+    //                << " the sampled cell/face " << sampleLocations[sampleI]
-                    << " Set the debug flag to get a dump of sampled cells."
+    //                << endl
-                    << endl;
+    //                << "    is not on a plane " << avgOffset
-                break;
+    //                << " offset from the patch." << endl
-            }
+    //                << "    You might want to shift your plane offset."
-        }
+    //                << " Set the debug flag to get a dump of sampled cells."
-    }
+    //                << endl;
+    //            break;
+    //        }
+    //    }
+    //}
 
 
     // Determine schedule.
@@ -556,14 +588,36 @@ void Foam::directMappedPatchBase::calcMapping() const
 
 Foam::directMappedPatchBase::directMappedPatchBase
 (
-     const polyPatch& pp
+    const polyPatch& pp
 )
 :
     patch_(pp),
     sampleRegion_(patch_.boundaryMesh().mesh().name()),
     mode_(NEARESTPATCHFACE),
     samplePatch_("none"),
+    uniformOffset_(true),
     offset_(vector::zero),
+    offsets_(pp.size(), offset_),
+    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
+    mapPtr_(NULL)
+{}
+
+
+Foam::directMappedPatchBase::directMappedPatchBase
+(
+    const polyPatch& pp,
+    const word& sampleRegion,
+    const sampleMode mode,
+    const word& samplePatch,
+    const vectorField& offsets
+)
+:
+    patch_(pp),
+    sampleRegion_(sampleRegion),
+    mode_(mode),
+    samplePatch_(samplePatch),
+    uniformOffset_(false),
+    offsets_(offsets),
     sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
     mapPtr_(NULL)
 {}
@@ -586,7 +640,19 @@ Foam::directMappedPatchBase::directMappedPatchBase
     ),
     mode_(sampleModeNames_.read(dict.lookup("sampleMode"))),
     samplePatch_(dict.lookup("samplePatch")),
-    offset_(dict.lookup("offset")),
+    uniformOffset_(dict.found("offset")),
+    offset_
+    (
+        uniformOffset_
+      ? point(dict.lookup("offset"))
+      : vector::zero
+    ),
+    offsets_
+    (
+        uniformOffset_
+      ? pointField(patch_.size(), offset_)
+      : dict.lookup("offsets")
+    ),
     sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
     mapPtr_(NULL)
 {}
@@ -602,7 +668,28 @@ Foam::directMappedPatchBase::directMappedPatchBase
     sampleRegion_(dmp.sampleRegion_),
     mode_(dmp.mode_),
     samplePatch_(dmp.samplePatch_),
+    uniformOffset_(dmp.uniformOffset_),
     offset_(dmp.offset_),
+    offsets_(dmp.offsets_),
+    sameRegion_(dmp.sameRegion_),
+    mapPtr_(NULL)
+{}
+
+
+Foam::directMappedPatchBase::directMappedPatchBase
+(
+    const polyPatch& pp,
+    const directMappedPatchBase& dmp,
+    const unallocLabelList& mapAddressing
+)
+:
+    patch_(pp),
+    sampleRegion_(dmp.sampleRegion_),
+    mode_(dmp.mode_),
+    samplePatch_(dmp.samplePatch_),
+    uniformOffset_(dmp.uniformOffset_),
+    offset_(dmp.offset_),
+    offsets_(dmp.offsets_, mapAddressing),
     sameRegion_(dmp.sameRegion_),
     mapPtr_(NULL)
 {}
@@ -660,7 +747,14 @@ void Foam::directMappedPatchBase::write(Ostream& os) const
         << token::END_STATEMENT << nl;
     os.writeKeyword("samplePatch") << samplePatch_
         << token::END_STATEMENT << nl;
-    os.writeKeyword("offset") << offset_ << token::END_STATEMENT << nl;
+    if (uniformOffset_)
+    {
+        os.writeKeyword("offset") << offset_ << token::END_STATEMENT << nl;
+    }
+    else
+    {
+        os.writeKeyword("offsets") << offsets_ << token::END_STATEMENT << nl;
+    }
 }
 
 

src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.H

@@ -87,12 +87,18 @@ private:
         //- What to sample
         const sampleMode mode_;
 
-        //- Patch (only if NEARESTBOUNDARY)
+        //- Patch (only if NEARESTPATCHFACE)
         const word samplePatch_;
 
-        //- Offset vector
+        //- For backwards compatibility : reading/writing of uniform offset.
+        const bool uniformOffset_;
+
+        //- Offset vector (uniform)
         const vector offset_;
 
+        //- Offset vector
+        const vectorField offsets_;
+
         //- Same region
         const bool sameRegion_;
 
@@ -130,36 +136,6 @@ private:
         void calcMapping() const;
 
 
-    // Private class
-
-        //- Private class for finding nearest
-        //  - point+local index
-        //  - sqr(distance)
-        //  - processor
-        typedef Tuple2<pointIndexHit, Tuple2<scalar, label> > nearInfo;
-
-        class nearestEqOp
-        {
-
-        public:
-
-            void operator()(nearInfo& x, const nearInfo& y) const
-            {
-                if (y.first().hit())
-                {
-                    if (!x.first().hit())
-                    {
-                        x = y;
-                    }
-                    else if (y.second().first() < x.second().first())
-                    {
-                        x = y;
-                    }
-                }
-            }
-        };
-
-
 public:
 
     //- Runtime type information
@@ -168,15 +144,33 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from patch
         directMappedPatchBase(const polyPatch&);
 
+        //- Construct from components
+        directMappedPatchBase
+        (
+            const polyPatch& pp,
+            const word& sampleRegion,
+            const sampleMode sampleMode,
+            const word& samplePatch,
+            const vectorField& offset
+        );
+
         //- Construct from dictionary
         directMappedPatchBase(const polyPatch&, const dictionary&);
 
         //- Construct as copy, resetting patch
         directMappedPatchBase(const polyPatch&, const directMappedPatchBase&);
 
+        //- Construct as copy, resetting patch, map original data
+        directMappedPatchBase
+        (
+            const polyPatch&,
+            const directMappedPatchBase&,
+            const unallocLabelList& mapAddressing
+        );
+
 
     //- Destructor
     virtual ~directMappedPatchBase();
@@ -210,6 +204,12 @@ public:
             return offset_;
         }
 
+        //- Offset vector (from patch faces to destination mesh objects)
+        const vectorField& offsets() const
+        {
+            return offsets_;
+        }
+
         //- Return reference to the parallel distribution map
         const mapDistribute& map() const
         {

src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C

@@ -57,6 +57,31 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 {}
 
 
+Foam::directMappedPolyPatch::directMappedPolyPatch
+(
+    const word& name,
+    const label size,
+    const label start,
+    const label index,
+    const word& sampleRegion,
+    const directMappedPatchBase::sampleMode mode,
+    const word& samplePatch,
+    const vectorField& offset,
+    const polyBoundaryMesh& bm
+)
+:
+    polyPatch(name, size, start, index, bm),
+    directMappedPatchBase
+    (
+        static_cast<const polyPatch&>(*this),
+        sampleRegion,
+        mode,
+        samplePatch,
+        offset
+    )
+{}
+
+
 Foam::directMappedPolyPatch::directMappedPolyPatch
 (
     const word& name,
@@ -95,6 +120,20 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 {}
 
 
+Foam::directMappedPolyPatch::directMappedPolyPatch
+(
+    const directMappedPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    polyPatch(pp, bm, index, mapAddressing, newStart),
+    directMappedPatchBase(*this, pp, mapAddressing)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::directMappedPolyPatch::~directMappedPolyPatch()

src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H

@@ -101,6 +101,20 @@ public:
             const polyBoundaryMesh& bm
         );
 
+        //- Construct from components
+        directMappedPolyPatch
+        (
+            const word& name,
+            const label size,
+            const label start,
+            const label index,
+            const word& sampleRegion,
+            const directMappedPatchBase::sampleMode mode,
+            const word& samplePatch,
+            const vectorField& offset,
+            const polyBoundaryMesh& bm
+        );
+
         //- Construct from dictionary
         directMappedPolyPatch
         (
@@ -128,6 +142,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        directMappedPolyPatch
+        (
+            const directMappedPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -150,6 +174,29 @@ public:
             );
         }
 
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new directMappedPolyPatch
+                (
+                    *this,
+                    bm,
+                    index,
+                    mapAddressing,
+                    newStart
+                )
+            );
+        }
+
 
     //- Destructor
     virtual ~directMappedPolyPatch();

src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.C

@@ -62,6 +62,31 @@ Foam::directMappedWallPolyPatch::directMappedWallPolyPatch
 {}
 
 
+Foam::directMappedWallPolyPatch::directMappedWallPolyPatch
+(
+    const word& name,
+    const label size,
+    const label start,
+    const label index,
+    const word& sampleRegion,
+    const directMappedPatchBase::sampleMode mode,
+    const word& samplePatch,
+    const vectorField& offset,
+    const polyBoundaryMesh& bm
+)
+:
+    wallPolyPatch(name, size, start, index, bm),
+    directMappedPatchBase
+    (
+        static_cast<const polyPatch&>(*this),
+        sampleRegion,
+        mode,
+        samplePatch,
+        offset
+    )
+{}
+
+
 Foam::directMappedWallPolyPatch::directMappedWallPolyPatch
 (
     const word& name,
@@ -100,6 +125,20 @@ Foam::directMappedWallPolyPatch::directMappedWallPolyPatch
 {}
 
 
+Foam::directMappedWallPolyPatch::directMappedWallPolyPatch
+(
+    const directMappedWallPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const unallocLabelList& mapAddressing,
+    const label newStart
+)
+:
+    wallPolyPatch(pp, bm, index, mapAddressing, newStart),
+    directMappedPatchBase(*this, pp, mapAddressing)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::directMappedWallPolyPatch::~directMappedWallPolyPatch()

src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.H

@@ -101,6 +101,20 @@ public:
             const polyBoundaryMesh& bm
         );
 
+        //- Construct from components
+        directMappedWallPolyPatch
+        (
+            const word& name,
+            const label size,
+            const label start,
+            const label index,
+            const word& sampleRegion,
+            const directMappedPatchBase::sampleMode mode,
+            const word& samplePatch,
+            const vectorField& offset,
+            const polyBoundaryMesh& bm
+        );
+
         //- Construct from dictionary
         directMappedWallPolyPatch
         (
@@ -128,6 +142,16 @@ public:
             const label newStart
         );
 
+        //- Construct given the original patch and a map
+        directMappedWallPolyPatch
+        (
+            const directMappedWallPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        );
+
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
         {
@@ -157,6 +181,29 @@ public:
             );
         }
 
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const unallocLabelList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new directMappedWallPolyPatch
+                (
+                    *this,
+                    bm,
+                    index,
+                    mapAddressing,
+                    newStart
+                )
+            );
+        }
+
 
     //- Destructor
     virtual ~directMappedWallPolyPatch();

src/meshTools/searchableSurface/searchableSurfaceCollection.H

@@ -125,6 +125,31 @@ public:
 
     // Member Functions
 
+        //- scaling vector per subsurface
+        const vectorField& scale() const
+        {
+            return scale_;
+        }
+
+        //- scaling vector per subsurface
+        vectorField& scale()
+        {
+            return scale_;
+        }
+
+        //- coordinate system per subsurface
+        const PtrList<coordinateSystem>& transform() const
+        {
+            return transform_;
+        }
+
+        //- coordinate system per subsurface
+        PtrList<coordinateSystem>& transform()
+        {
+            return transform_;
+        }
+
+
         virtual const wordList& regions() const;
 
         //- Whether supports volume type below

src/meshTools/sets/cellSources/cellToCell/cellToCell.C

@@ -106,7 +106,7 @@ void Foam::cellToCell::applyToSet
 {
     if ((action == topoSetSource::ADD) || (action == topoSetSource::NEW))
     {
-        Pout<< "    Adding all elements of cellSet " << setName_ << " ..."
+        Info<< "    Adding all elements of cellSet " << setName_ << " ..."
             << endl;
 
         // Load the set
@@ -116,7 +116,7 @@ void Foam::cellToCell::applyToSet
     }
     else if (action == topoSetSource::DELETE)
     {
-        Pout<< "    Removing all elements of cellSet " << setName_ << " ..."
+        Info<< "    Removing all elements of cellSet " << setName_ << " ..."
             << endl;
 
         // Load the set

src/meshTools/sets/cellZoneSources/setToCellZone/setToCellZone.C

@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "setToCellZone.H"
+#include "polyMesh.H"
+#include "cellZoneSet.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(setToCellZone, 0);
+
+addToRunTimeSelectionTable(topoSetSource, setToCellZone, word);
+
+addToRunTimeSelectionTable(topoSetSource, setToCellZone, istream);
+
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::setToCellZone::usage_
+(
+    setToCellZone::typeName,
+    "\n    Usage: setToCellZone <cellSet>\n\n"
+    "    Select all cells in the cellSet.\n\n"
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::setToCellZone::setToCellZone
+(
+    const polyMesh& mesh,
+    const word& setName
+)
+:
+    topoSetSource(mesh),
+    setName_(setName)
+{}
+
+
+// Construct from dictionary
+Foam::setToCellZone::setToCellZone
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    topoSetSource(mesh),
+    setName_(dict.lookup("set"))
+{}
+
+
+// Construct from Istream
+Foam::setToCellZone::setToCellZone
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetSource(mesh),
+    setName_(checkIs(is))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::setToCellZone::~setToCellZone()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::setToCellZone::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if (!isA<cellZoneSet>(set))
+    {
+        WarningIn
+        (
+            "setToCellZone::applyToSet(const topoSetSource::setAction"
+            ", topoSet"
+        )   << "Operation only allowed on a cellZoneSet." << endl;
+    }
+    else
+    {
+        cellZoneSet& fzSet = refCast<cellZoneSet>(set);
+
+        if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+        {
+            Info<< "    Adding all cells from cellSet " << setName_
+                << " ..." << endl;
+
+            // Load the sets
+            cellSet fSet(mesh_, setName_);
+
+            // Start off from copy
+            DynamicList<label> newAddressing(fzSet.addressing());
+
+            forAllConstIter(cellSet, fSet, iter)
+            {
+                label cellI = iter.key();
+
+                if (!fzSet.found(cellI))
+                {
+                    newAddressing.append(cellI);
+                }
+            }
+
+            fzSet.addressing().transfer(newAddressing);
+            fzSet.updateSet();
+        }
+        else if (action == topoSetSource::DELETE)
+        {
+            Info<< "    Removing all cells from cellSet " << setName_
+                << " ..." << endl;
+
+            // Load the set
+            cellSet loadedSet(mesh_, setName_);
+
+            // Start off empty
+            DynamicList<label> newAddressing(fzSet.addressing().size());
+
+            forAll(fzSet.addressing(), i)
+            {
+                if (!loadedSet.found(fzSet.addressing()[i]))
+                {
+                    newAddressing.append(fzSet.addressing()[i]);
+                }
+            }
+            fzSet.addressing().transfer(newAddressing);
+            fzSet.updateSet();
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/meshTools/sets/cellZoneSources/setToCellZone/setToCellZone.H

@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::setToCellZone
+
+Description
+    A topoSetSource to select cells based on usage in a cellSet.
+
+SourceFiles
+    setToCellZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef setToCellZone_H
+#define setToCellZone_H
+
+#include "topoSetSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class setToCellZone Declaration
+\*---------------------------------------------------------------------------*/
+
+class setToCellZone
+:
+    public topoSetSource
+{
+    // Private data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- Name of set to use
+        word setName_;
+
+public:
+
+    //- Runtime type information
+    TypeName("setToCellZone");
+
+    // Constructors
+
+        //- Construct from components
+        setToCellZone
+        (
+            const polyMesh& mesh,
+            const word& setName
+        );
+
+        //- Construct from dictionary
+        setToCellZone
+        (
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct from Istream
+        setToCellZone
+        (
+            const polyMesh& mesh,
+            Istream&
+        );
+
+
+    // Destructor
+
+        virtual ~setToCellZone();
+
+
+    // Member Functions
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet&
+        ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/sets/faceSources/cellToFace/cellToFace.C

@@ -217,14 +217,14 @@ void Foam::cellToFace::applyToSet
 {
     if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
     {
-        Pout<< "    Adding faces according to cellSet " << setName_
+        Info<< "    Adding faces according to cellSet " << setName_
             << " ..." << endl;
 
         combine(set, true);
     }
     else if (action == topoSetSource::DELETE)
     {
-        Pout<< "    Removing faces according to cellSet " << setName_
+        Info<< "    Removing faces according to cellSet " << setName_
             << " ..." << endl;
 
         combine(set, false);

src/meshTools/sets/faceZoneSources/faceZoneToFaceZone/faceZoneToFaceZone.C

@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "faceZoneToFaceZone.H"
+#include "polyMesh.H"
+#include "faceZoneSet.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(faceZoneToFaceZone, 0);
+
+addToRunTimeSelectionTable(topoSetSource, faceZoneToFaceZone, word);
+
+addToRunTimeSelectionTable(topoSetSource, faceZoneToFaceZone, istream);
+
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::faceZoneToFaceZone::usage_
+(
+    faceZoneToFaceZone::typeName,
+    "\n    Usage: faceZoneToFaceZone <faceZone>\n\n"
+    "    Select all faces in the faceZone\n\n"
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::faceZoneToFaceZone::faceZoneToFaceZone
+(
+    const polyMesh& mesh,
+    const word& setName
+)
+:
+    topoSetSource(mesh),
+    setName_(setName)
+{}
+
+
+// Construct from dictionary
+Foam::faceZoneToFaceZone::faceZoneToFaceZone
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    topoSetSource(mesh),
+    setName_(dict.lookup("zone"))
+{}
+
+
+// Construct from Istream
+Foam::faceZoneToFaceZone::faceZoneToFaceZone
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetSource(mesh),
+    setName_(checkIs(is))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::faceZoneToFaceZone::~faceZoneToFaceZone()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::faceZoneToFaceZone::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if (!isA<faceZoneSet>(set))
+    {
+        WarningIn
+        (
+            "faceZoneToFaceZone::applyToSet(const topoSetSource::setAction"
+            ", topoSet"
+        )   << "Operation only allowed on a faceZoneSet." << endl;
+    }
+    else
+    {
+        faceZoneSet& fSet = refCast<faceZoneSet>(set);
+
+        if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+        {
+            Info<< "    Adding all faces from faceZone " << setName_ << " ..."
+                << endl;
+
+            // Load the set
+            faceZoneSet loadedSet(mesh_, setName_);
+
+            DynamicList<label> newAddressing(fSet.addressing());
+            DynamicList<bool> newFlipMap(fSet.flipMap());
+
+            forAll(loadedSet.addressing(), i)
+            {
+                if (!fSet.found(loadedSet.addressing()[i]))
+                {
+                    newAddressing.append(loadedSet.addressing()[i]);
+                    newFlipMap.append(loadedSet.flipMap()[i]);
+                }
+            }
+            fSet.addressing().transfer(newAddressing);
+            fSet.flipMap().transfer(newFlipMap);
+            fSet.updateSet();
+        }
+        else if (action == topoSetSource::DELETE)
+        {
+            Info<< "    Removing all faces from faceZone " << setName_ << " ..."
+                << endl;
+
+            // Load the set
+            faceZoneSet loadedSet(mesh_, setName_);
+
+            DynamicList<label> newAddressing(fSet.addressing().size());
+            DynamicList<bool> newFlipMap(fSet.flipMap().size());
+
+            forAll(fSet.addressing(), i)
+            {
+                if (!loadedSet.found(fSet.addressing()[i]))
+                {
+                    newAddressing.append(fSet.addressing()[i]);
+                    newFlipMap.append(fSet.flipMap()[i]);
+                }
+            }
+            fSet.addressing().transfer(newAddressing);
+            fSet.flipMap().transfer(newFlipMap);
+            fSet.updateSet();
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/meshTools/sets/faceZoneSources/faceZoneToFaceZone/faceZoneToFaceZone.H

@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::faceZoneToFaceZone
+
+Description
+    A topoSetSource to select faces based on usage in another faceSet.
+
+SourceFiles
+    faceZoneToFaceZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faceZoneToFaceZone_H
+#define faceZoneToFaceZone_H
+
+#include "topoSetSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class faceZoneToFaceZone Declaration
+\*---------------------------------------------------------------------------*/
+
+class faceZoneToFaceZone
+:
+    public topoSetSource
+{
+    // Private data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- Name of set to use
+        word setName_;
+
+public:
+
+    //- Runtime type information
+    TypeName("faceZoneToFaceZone");
+
+    // Constructors
+
+        //- Construct from components
+        faceZoneToFaceZone
+        (
+            const polyMesh& mesh,
+            const word& setName
+        );
+
+        //- Construct from dictionary
+        faceZoneToFaceZone
+        (
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct from Istream
+        faceZoneToFaceZone
+        (
+            const polyMesh& mesh,
+            Istream&
+        );
+
+
+    // Destructor
+
+        virtual ~faceZoneToFaceZone();
+
+
+    // Member Functions
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet&
+        ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/sets/faceZoneSources/setToFaceZone/setToFaceZone.C

@@ -0,0 +1,175 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "setToFaceZone.H"
+#include "polyMesh.H"
+#include "faceZoneSet.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(setToFaceZone, 0);
+
+addToRunTimeSelectionTable(topoSetSource, setToFaceZone, word);
+
+addToRunTimeSelectionTable(topoSetSource, setToFaceZone, istream);
+
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::setToFaceZone::usage_
+(
+    setToFaceZone::typeName,
+    "\n    Usage: setToFaceZone <faceSet>\n\n"
+    "    Select all faces in the faceSet."
+    " Sets flipMap.\n\n"
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::setToFaceZone::setToFaceZone
+(
+    const polyMesh& mesh,
+    const word& setName
+)
+:
+    topoSetSource(mesh),
+    setName_(setName)
+{}
+
+
+// Construct from dictionary
+Foam::setToFaceZone::setToFaceZone
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    topoSetSource(mesh),
+    setName_(dict.lookup("faceSet"))
+{}
+
+
+// Construct from Istream
+Foam::setToFaceZone::setToFaceZone
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetSource(mesh),
+    setName_(checkIs(is))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::setToFaceZone::~setToFaceZone()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::setToFaceZone::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if (!isA<faceZoneSet>(set))
+    {
+        WarningIn
+        (
+            "setToFaceZone::applyToSet(const topoSetSource::setAction"
+            ", topoSet"
+        )   << "Operation only allowed on a faceZoneSet." << endl;
+    }
+    else
+    {
+        faceZoneSet& fzSet = refCast<faceZoneSet>(set);
+
+        if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+        {
+            Info<< "    Adding all faces from faceSet " << setName_
+                << " ..." << endl;
+
+            // Load the sets
+            faceSet fSet(mesh_, setName_);
+
+            // Start off from copy
+            DynamicList<label> newAddressing(fzSet.addressing());
+            DynamicList<bool> newFlipMap(fzSet.flipMap());
+
+            forAllConstIter(faceSet, fSet, iter)
+            {
+                label faceI = iter.key();
+
+                if (!fzSet.found(faceI))
+                {
+                    newAddressing.append(faceI);
+                    newFlipMap.append(false);
+                }
+            }
+
+            fzSet.addressing().transfer(newAddressing);
+            fzSet.flipMap().transfer(newFlipMap);
+            fzSet.updateSet();
+        }
+        else if (action == topoSetSource::DELETE)
+        {
+            Info<< "    Removing all faces from faceSet " << setName_
+                << " ..." << endl;
+
+            // Load the set
+            faceSet loadedSet(mesh_, setName_);
+
+            // Start off empty
+            DynamicList<label> newAddressing(fzSet.addressing().size());
+            DynamicList<bool> newFlipMap(fzSet.flipMap().size());
+
+            forAll(fzSet.addressing(), i)
+            {
+                if (!loadedSet.found(fzSet.addressing()[i]))
+                {
+                    newAddressing.append(fzSet.addressing()[i]);
+                    newFlipMap.append(fzSet.flipMap()[i]);
+                }
+            }
+            fzSet.addressing().transfer(newAddressing);
+            fzSet.flipMap().transfer(newFlipMap);
+            fzSet.updateSet();
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/meshTools/sets/faceZoneSources/setToFaceZone/setToFaceZone.H

@@ -0,0 +1,116 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::setToFaceZone
+
+Description
+    A topoSetSource to select faces based on usage in a faceSet. Sets flipMap
+    to true.
+
+SourceFiles
+    setToFaceZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef setToFaceZone_H
+#define setToFaceZone_H
+
+#include "topoSetSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class setToFaceZone Declaration
+\*---------------------------------------------------------------------------*/
+
+class setToFaceZone
+:
+    public topoSetSource
+{
+    // Private data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- Name of set to use
+        word setName_;
+
+public:
+
+    //- Runtime type information
+    TypeName("setToFaceZone");
+
+    // Constructors
+
+        //- Construct from components
+        setToFaceZone
+        (
+            const polyMesh& mesh,
+            const word& setName
+        );
+
+        //- Construct from dictionary
+        setToFaceZone
+        (
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct from Istream
+        setToFaceZone
+        (
+            const polyMesh& mesh,
+            Istream&
+        );
+
+
+    // Destructor
+
+        virtual ~setToFaceZone();
+
+
+    // Member Functions
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet&
+        ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.C

@@ -0,0 +1,222 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "setsToFaceZone.H"
+#include "polyMesh.H"
+#include "faceZoneSet.H"
+#include "cellSet.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(setsToFaceZone, 0);
+
+addToRunTimeSelectionTable(topoSetSource, setsToFaceZone, word);
+
+addToRunTimeSelectionTable(topoSetSource, setsToFaceZone, istream);
+
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::setsToFaceZone::usage_
+(
+    setsToFaceZone::typeName,
+    "\n    Usage: setsToFaceZone <faceSet> <slaveCellSet>\n\n"
+    "    Select all faces in the faceSet."
+    " Orientated so slave side is in cellSet.\n\n"
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::setsToFaceZone::setsToFaceZone
+(
+    const polyMesh& mesh,
+    const word& faceSetName,
+    const word& cellSetName
+)
+:
+    topoSetSource(mesh),
+    faceSetName_(faceSetName),
+    cellSetName_(cellSetName)
+{}
+
+
+// Construct from dictionary
+Foam::setsToFaceZone::setsToFaceZone
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    topoSetSource(mesh),
+    faceSetName_(dict.lookup("faceSet")),
+    cellSetName_(dict.lookup("cellSet"))
+{}
+
+
+// Construct from Istream
+Foam::setsToFaceZone::setsToFaceZone
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetSource(mesh),
+    faceSetName_(checkIs(is)),
+    cellSetName_(checkIs(is))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::setsToFaceZone::~setsToFaceZone()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::setsToFaceZone::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if (!isA<faceZoneSet>(set))
+    {
+        WarningIn
+        (
+            "setsToFaceZone::applyToSet(const topoSetSource::setAction"
+            ", topoSet"
+        )   << "Operation only allowed on a faceZoneSet." << endl;
+    }
+    else
+    {
+        faceZoneSet& fzSet = refCast<faceZoneSet>(set);
+
+        if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+        {
+            Info<< "    Adding all faces from faceSet " << faceSetName_
+                << " ..." << endl;
+
+            // Load the sets
+            faceSet fSet(mesh_, faceSetName_);
+            cellSet cSet(mesh_, cellSetName_);
+
+            // Start off from copy
+            DynamicList<label> newAddressing(fzSet.addressing());
+            DynamicList<bool> newFlipMap(fzSet.flipMap());
+
+            forAllConstIter(faceSet, fSet, iter)
+            {
+                label faceI = iter.key();
+
+                if (!fzSet.found(faceI))
+                {
+                    bool flip = false;
+
+                    label own = mesh_.faceOwner()[faceI];
+                    bool ownFound = cSet.found(own);
+
+                    if (mesh_.isInternalFace(faceI))
+                    {
+                        label nei = mesh_.faceNeighbour()[faceI];
+                        bool neiFound = cSet.found(nei);
+
+                        if (ownFound && !neiFound)
+                        {
+                            flip = false;
+                        }
+                        else if(!ownFound && neiFound)
+                        {
+                            flip = true;
+                        }
+                        else
+                        {
+                            WarningIn
+                            (
+                                "setsToFaceZone::applyToSet"
+                                "(const topoSetSource::setAction, topoSet)"
+                            )   << "One of owner or neighbour of internal face "
+                                << faceI << " should be in cellSet "
+                                << cSet.name()
+                                << " to be able to determine orientation."
+                                << endl
+                                << "Face:" << faceI << " own:" << own
+                                << " OwnInCellSet:" << ownFound
+                                << " nei:" << nei
+                                << " NeiInCellSet:" << neiFound
+                                << endl;
+                        }
+                    }
+                    else
+                    {
+                        flip = !ownFound;
+                    }
+
+                    newAddressing.append(faceI);
+                    newFlipMap.append(flip);
+                }
+            }
+
+            fzSet.addressing().transfer(newAddressing);
+            fzSet.flipMap().transfer(newFlipMap);
+            fzSet.updateSet();
+        }
+        else if (action == topoSetSource::DELETE)
+        {
+            Info<< "    Removing all faces from faceSet " << faceSetName_
+                << " ..." << endl;
+
+            // Load the set
+            faceZoneSet loadedSet(mesh_, faceSetName_);
+
+            // Start off empty
+            DynamicList<label> newAddressing(fzSet.addressing().size());
+            DynamicList<bool> newFlipMap(fzSet.flipMap().size());
+
+            forAll(fzSet.addressing(), i)
+            {
+                if (!loadedSet.found(fzSet.addressing()[i]))
+                {
+                    newAddressing.append(fzSet.addressing()[i]);
+                    newFlipMap.append(fzSet.flipMap()[i]);
+                }
+            }
+            fzSet.addressing().transfer(newAddressing);
+            fzSet.flipMap().transfer(newFlipMap);
+            fzSet.updateSet();
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.H

@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::setsToFaceZone
+
+Description
+    A topoSetSource to select faces based on usage in a faceSet and cellSet
+
+SourceFiles
+    setsToFaceZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef setsToFaceZone_H
+#define setsToFaceZone_H
+
+#include "topoSetSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class setsToFaceZone Declaration
+\*---------------------------------------------------------------------------*/
+
+class setsToFaceZone
+:
+    public topoSetSource
+{
+    // Private data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- Name of set to use
+        word faceSetName_;
+
+        //- Name of set to use
+        word cellSetName_;
+
+public:
+
+    //- Runtime type information
+    TypeName("setsToFaceZone");
+
+    // Constructors
+
+        //- Construct from components
+        setsToFaceZone
+        (
+            const polyMesh& mesh,
+            const word& faceSetName,
+            const word& cellSetName
+        );
+
+        //- Construct from dictionary
+        setsToFaceZone
+        (
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct from Istream
+        setsToFaceZone
+        (
+            const polyMesh& mesh,
+            Istream&
+        );
+
+
+    // Destructor
+
+        virtual ~setsToFaceZone();
+
+
+    // Member Functions
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet&
+        ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/sets/pointSources/cellToPoint/cellToPoint.C

@@ -151,13 +151,13 @@ void Foam::cellToPoint::applyToSet
 {
     if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
     {
-        Pout<< "    Adding from " << setName_ << " ..." << endl;
+        Info<< "    Adding from " << setName_ << " ..." << endl;
 
         combine(set, true);
     }
     else if (action == topoSetSource::DELETE)
     {
-        Pout<< "    Removing from " << setName_ << " ..." << endl;
+        Info<< "    Removing from " << setName_ << " ..." << endl;
 
         combine(set, false);
     }

src/meshTools/sets/pointZoneSources/setToPointZone/setToPointZone.C

@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "setToPointZone.H"
+#include "polyMesh.H"
+#include "pointZoneSet.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(setToPointZone, 0);
+
+addToRunTimeSelectionTable(topoSetSource, setToPointZone, word);
+
+addToRunTimeSelectionTable(topoSetSource, setToPointZone, istream);
+
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::setToPointZone::usage_
+(
+    setToPointZone::typeName,
+    "\n    Usage: setToPointZone <pointSet>\n\n"
+    "    Select all points in the pointSet.\n\n"
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::setToPointZone::setToPointZone
+(
+    const polyMesh& mesh,
+    const word& setName
+)
+:
+    topoSetSource(mesh),
+    setName_(setName)
+{}
+
+
+// Construct from dictionary
+Foam::setToPointZone::setToPointZone
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    topoSetSource(mesh),
+    setName_(dict.lookup("set"))
+{}
+
+
+// Construct from Istream
+Foam::setToPointZone::setToPointZone
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetSource(mesh),
+    setName_(checkIs(is))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::setToPointZone::~setToPointZone()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::setToPointZone::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if (!isA<pointZoneSet>(set))
+    {
+        WarningIn
+        (
+            "setToPointZone::applyToSet(const topoSetSource::setAction"
+            ", topoSet"
+        )   << "Operation only allowed on a pointZoneSet." << endl;
+    }
+    else
+    {
+        pointZoneSet& fzSet = refCast<pointZoneSet>(set);
+
+        if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+        {
+            Info<< "    Adding all points from pointSet " << setName_
+                << " ..." << endl;
+
+            // Load the sets
+            pointSet fSet(mesh_, setName_);
+
+            // Start off from copy
+            DynamicList<label> newAddressing(fzSet.addressing());
+
+            forAllConstIter(pointSet, fSet, iter)
+            {
+                label pointI = iter.key();
+
+                if (!fzSet.found(pointI))
+                {
+                    newAddressing.append(pointI);
+                }
+            }
+
+            fzSet.addressing().transfer(newAddressing);
+            fzSet.updateSet();
+        }
+        else if (action == topoSetSource::DELETE)
+        {
+            Info<< "    Removing all points from pointSet " << setName_
+                << " ..." << endl;
+
+            // Load the set
+            pointSet loadedSet(mesh_, setName_);
+
+            // Start off empty
+            DynamicList<label> newAddressing(fzSet.addressing().size());
+
+            forAll(fzSet.addressing(), i)
+            {
+                if (!loadedSet.found(fzSet.addressing()[i]))
+                {
+                    newAddressing.append(fzSet.addressing()[i]);
+                }
+            }
+            fzSet.addressing().transfer(newAddressing);
+            fzSet.updateSet();
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/meshTools/sets/pointZoneSources/setToPointZone/setToPointZone.H

@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::setToPointZone
+
+Description
+    A topoSetSource to select points based on usage in a pointSet.
+
+SourceFiles
+    setToPointZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef setToPointZone_H
+#define setToPointZone_H
+
+#include "topoSetSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class setToPointZone Declaration
+\*---------------------------------------------------------------------------*/
+
+class setToPointZone
+:
+    public topoSetSource
+{
+    // Private data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- Name of set to use
+        word setName_;
+
+public:
+
+    //- Runtime type information
+    TypeName("setToPointZone");
+
+    // Constructors
+
+        //- Construct from components
+        setToPointZone
+        (
+            const polyMesh& mesh,
+            const word& setName
+        );
+
+        //- Construct from dictionary
+        setToPointZone
+        (
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct from Istream
+        setToPointZone
+        (
+            const polyMesh& mesh,
+            Istream&
+        );
+
+
+    // Destructor
+
+        virtual ~setToPointZone();
+
+
+    // Member Functions
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet&
+        ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/sets/topoSetSource/topoSetSource.C

@@ -101,7 +101,7 @@ autoPtr<topoSetSource> topoSetSource::New
 Foam::HashTable<Foam::string>* Foam::topoSetSource::usageTablePtr_ = NULL;
 
 template<>
-const char* Foam::NamedEnum<Foam::topoSetSource::setAction, 7>::names[] =
+const char* Foam::NamedEnum<Foam::topoSetSource::setAction, 8>::names[] =
 {
     "clear",
     "new",
@@ -109,11 +109,12 @@ const char* Foam::NamedEnum<Foam::topoSetSource::setAction, 7>::names[] =
     "add",
     "delete",
     "subset",
-    "list"
+    "list",
+    "remove"
 };
 
 
-const Foam::NamedEnum<Foam::topoSetSource::setAction, 7>
+const Foam::NamedEnum<Foam::topoSetSource::setAction, 8>
     Foam::topoSetSource::actionNames_;
 
 

src/meshTools/sets/topoSetSource/topoSetSource.H

@@ -77,7 +77,8 @@ public:
             ADD,
             DELETE,
             SUBSET,
-            LIST
+            LIST,
+            REMOVE
         };
 
 protected:
@@ -120,7 +121,7 @@ protected:
 
 private:
 
-        static const NamedEnum<setAction, 7> actionNames_;
+        static const NamedEnum<setAction, 8> actionNames_;
 
         static const string illegalSource_;
 

src/meshTools/sets/topoSets/cellSet.C

@@ -41,6 +41,7 @@ defineTypeNameAndDebug(cellSet, 0);
 
 addToRunTimeSelectionTable(topoSet, cellSet, word);
 addToRunTimeSelectionTable(topoSet, cellSet, size);
+addToRunTimeSelectionTable(topoSet, cellSet, set);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -78,6 +79,18 @@ cellSet::cellSet
 {}
 
 
+cellSet::cellSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const topoSet& set,
+    writeOption w
+)
+:
+    topoSet(mesh, name, set, w)
+{}
+
+
 cellSet::cellSet
 (
     const polyMesh& mesh,

src/meshTools/sets/topoSets/cellSet.H

@@ -88,6 +88,15 @@ public:
             writeOption w=NO_WRITE
         );
 
+        //- Construct from existing set
+        cellSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const topoSet&,
+            writeOption w=NO_WRITE
+        );
+
         //- Construct from labelHashSet
         cellSet
         (

src/meshTools/sets/topoSets/cellZoneSet.C

@@ -0,0 +1,314 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "cellZoneSet.H"
+#include "mapPolyMesh.H"
+#include "polyMesh.H"
+#include "processorPolyPatch.H"
+#include "cyclicPolyPatch.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(cellZoneSet, 0);
+
+addToRunTimeSelectionTable(topoSet, cellZoneSet, word);
+addToRunTimeSelectionTable(topoSet, cellZoneSet, size);
+addToRunTimeSelectionTable(topoSet, cellZoneSet, set);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void cellZoneSet::updateSet()
+{
+    labelList order;
+    sortedOrder(addressing_, order);
+    inplaceReorder(order, addressing_);
+
+    cellSet::clearStorage();
+    cellSet::resize(2*addressing_.size());
+    forAll(addressing_, i)
+    {
+        cellSet::insert(addressing_[i]);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+cellZoneSet::cellZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    readOption r,
+    writeOption w
+)
+:
+    cellSet(mesh, name, 1000),  // do not read cellSet
+    mesh_(mesh),
+    addressing_(0)
+{
+    const cellZoneMesh& cellZones = mesh.cellZones();
+    label zoneID = cellZones.findZoneID(name);
+
+    if
+    (
+        (r == IOobject::MUST_READ)
+     || (r == IOobject::READ_IF_PRESENT && zoneID != -1)
+    )
+    {
+        const cellZone& fz = cellZones[zoneID];
+        addressing_ = fz;
+    }
+
+    updateSet();
+
+    check(mesh.nCells());
+}
+
+
+cellZoneSet::cellZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const label size,
+    writeOption w
+)
+:
+    cellSet(mesh, name, size, w),
+    mesh_(mesh),
+    addressing_(0)
+{
+    updateSet();
+}
+
+
+cellZoneSet::cellZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const topoSet& set,
+    writeOption w
+)
+:
+    cellSet(mesh, name, set.size(), w),
+    mesh_(mesh),
+    addressing_(refCast<const cellZoneSet>(set).addressing())
+{
+    updateSet();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+cellZoneSet::~cellZoneSet()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void cellZoneSet::invert(const label maxLen)
+{
+    label n = 0;
+
+    for (label cellI = 0; cellI < maxLen; cellI++)
+    {
+        if (!found(cellI))
+        {
+            addressing_[n] = cellI;
+            n++;
+        }
+    }
+    addressing_.setSize(n);
+    updateSet();
+}
+
+
+void cellZoneSet::subset(const topoSet& set)
+{
+    DynamicList<label> newAddressing(addressing_.size());
+
+    const cellZoneSet& fSet = refCast<const cellZoneSet>(set);
+
+    forAll(fSet.addressing(), i)
+    {
+        label cellI = fSet.addressing()[i];
+
+        if (found(cellI))
+        {
+            newAddressing.append(cellI);
+        }
+    }
+
+    addressing_.transfer(newAddressing);
+    updateSet();
+}
+
+
+void cellZoneSet::addSet(const topoSet& set)
+{
+    DynamicList<label> newAddressing(addressing_);
+
+    const cellZoneSet& fSet = refCast<const cellZoneSet>(set);
+
+    forAll(fSet.addressing(), i)
+    {
+        label cellI = fSet.addressing()[i];
+
+        if (!found(cellI))
+        {
+            newAddressing.append(cellI);
+        }
+    }
+
+    addressing_.transfer(newAddressing);
+    updateSet();
+}
+
+
+void cellZoneSet::deleteSet(const topoSet& set)
+{
+    DynamicList<label> newAddressing(addressing_.size());
+
+    const cellZoneSet& fSet = refCast<const cellZoneSet>(set);
+
+    forAll(addressing_, i)
+    {
+        label cellI = addressing_[i];
+
+        if (!fSet.found(cellI))
+        {
+            // Not found in fSet so add
+            newAddressing.append(cellI);
+        }
+    }
+
+    addressing_.transfer(newAddressing);
+    updateSet();
+}
+
+
+void cellZoneSet::sync(const polyMesh& mesh)
+{}
+
+
+label cellZoneSet::maxSize(const polyMesh& mesh) const
+{
+    return mesh.nCells();
+}
+
+
+//- Write using given format, version and compression
+bool cellZoneSet::writeObject
+(
+    IOstream::streamFormat s,
+    IOstream::versionNumber v,
+    IOstream::compressionType c
+) const
+{
+    // Write shadow cellSet
+    word oldTypeName = typeName;
+    const_cast<word&>(type()) = cellSet::typeName;
+    bool ok = cellSet::writeObject(s, v, c);
+    const_cast<word&>(type()) = oldTypeName;
+
+    // Modify cellZone
+    cellZoneMesh& cellZones = const_cast<polyMesh&>(mesh_).cellZones();
+    label zoneID = cellZones.findZoneID(name());
+
+    if (zoneID == -1)
+    {
+        zoneID = cellZones.size();
+
+        cellZones.setSize(zoneID+1);
+        cellZones.set
+        (
+            zoneID,
+            new cellZone
+            (
+                name(),
+                addressing_,
+                zoneID,
+                cellZones
+            )
+        );
+    }
+    else
+    {
+        cellZones[zoneID] = addressing_;
+    }
+    cellZones.clearAddressing();
+
+    return ok && cellZones.write();
+}
+
+
+void cellZoneSet::updateMesh(const mapPolyMesh& morphMap)
+{
+    // cellZone
+    labelList newAddressing(addressing_.size());
+
+    label n = 0;
+    forAll(addressing_, i)
+    {
+        label cellI = addressing_[i];
+        label newCellI = morphMap.reverseCellMap()[cellI];
+        if (newCellI >= 0)
+        {
+            newAddressing[n] = newCellI;
+            n++;
+        }
+    }
+    newAddressing.setSize(n);
+
+    addressing_.transfer(newAddressing);
+
+    updateSet();
+}
+
+
+void cellZoneSet::writeDebug
+(
+    Ostream& os,
+    const primitiveMesh& mesh,
+    const label maxLen
+) const
+{
+    cellSet::writeDebug(os, mesh, maxLen);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/meshTools/sets/topoSets/cellZoneSet.H

@@ -0,0 +1,170 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::cellZoneSet
+
+Description
+    Like cellSet but updates cellZone when writing.
+
+SourceFiles
+    cellZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cellZoneSet_H
+#define cellZoneSet_H
+
+#include "cellSet.H"
+#include "boolList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class cellZoneSet Declaration
+\*---------------------------------------------------------------------------*/
+
+class cellZoneSet
+:
+    public cellSet
+{
+    // Private data
+
+        const polyMesh& mesh_;
+
+        labelList addressing_;
+
+public:
+
+    //- Runtime type information
+    TypeName("cellZoneSet");
+
+
+    // Constructors
+
+        //- Construct from objectRegistry and name
+        cellZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            readOption r=MUST_READ,
+            writeOption w=NO_WRITE
+        );
+
+        //- Construct from additional size of labelHashSet
+        cellZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const label,
+            writeOption w=NO_WRITE
+        );
+
+        //- Construct from existing set
+        cellZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const topoSet&,
+            writeOption w=NO_WRITE
+        );
+
+
+
+    // Destructor
+
+        virtual ~cellZoneSet();
+
+
+    // Member functions
+
+        const labelList& addressing() const
+        {
+            return addressing_;
+        }
+
+        labelList& addressing()
+        {
+            return addressing_;
+        }
+
+        //- Sort addressing and make cellSet part consistent with addressing
+        void updateSet();
+
+        //- Invert contents. (insert all members 0..maxLen-1 which were not in
+        //  set)
+        virtual void invert(const label maxLen);
+
+        //- Subset contents. Only elements present in both sets remain.
+        virtual void subset(const topoSet& set);
+
+        //- Add elements present in set.
+        virtual void addSet(const topoSet& set);
+
+        //- Delete elements present in set.
+        virtual void deleteSet(const topoSet& set);
+
+        //- Sync cellZoneSet across coupled patches.
+        virtual void sync(const polyMesh& mesh);
+
+        //- Write maxLen items with label and coordinates. 
+        virtual void writeDebug
+        (
+            Ostream& os,
+            const primitiveMesh&,
+            const label maxLen
+        ) const;
+
+        //- Write cellZone
+        virtual bool writeObject
+        (
+            IOstream::streamFormat,
+            IOstream::versionNumber,
+            IOstream::compressionType
+        ) const;
+
+        //- Update any stored data for new labels
+        virtual void updateMesh(const mapPolyMesh& morphMap);
+
+        //- Return max index+1.
+        virtual label maxSize(const polyMesh& mesh) const;
+
+
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/sets/topoSets/faceSet.C

@@ -43,6 +43,7 @@ defineTypeNameAndDebug(faceSet, 0);
 
 addToRunTimeSelectionTable(topoSet, faceSet, word);
 addToRunTimeSelectionTable(topoSet, faceSet, size);
+addToRunTimeSelectionTable(topoSet, faceSet, set);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -79,6 +80,18 @@ faceSet::faceSet
 {}
 
 
+faceSet::faceSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const topoSet& set,
+    writeOption w
+)
+:
+    topoSet(mesh, name, set, w)
+{}
+
+
 faceSet::faceSet
 (
     const polyMesh& mesh,

src/meshTools/sets/topoSets/faceSet.H

@@ -83,6 +83,15 @@ public:
             writeOption w=NO_WRITE
         );
 
+        //- Construct from existing set
+        faceSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const topoSet&,
+            writeOption w=NO_WRITE
+        );
+
         //- Construct from additional labelHashSet
         faceSet
         (

src/meshTools/sets/topoSets/faceZoneSet.C

@@ -0,0 +1,413 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "faceZoneSet.H"
+#include "mapPolyMesh.H"
+#include "polyMesh.H"
+#include "processorPolyPatch.H"
+#include "cyclicPolyPatch.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(faceZoneSet, 0);
+
+addToRunTimeSelectionTable(topoSet, faceZoneSet, word);
+addToRunTimeSelectionTable(topoSet, faceZoneSet, size);
+addToRunTimeSelectionTable(topoSet, faceZoneSet, set);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void faceZoneSet::updateSet()
+{
+    labelList order;
+    sortedOrder(addressing_, order);
+    inplaceReorder(order, addressing_);
+    inplaceReorder(order, flipMap_);
+
+    faceSet::clearStorage();
+    faceSet::resize(2*addressing_.size());
+    forAll(addressing_, i)
+    {
+        faceSet::insert(addressing_[i]);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+faceZoneSet::faceZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    readOption r,
+    writeOption w
+)
+:
+    faceSet(mesh, name, 1000),  // do not read faceSet
+    mesh_(mesh),
+    addressing_(0),
+    flipMap_(0)
+{
+    const faceZoneMesh& faceZones = mesh.faceZones();
+    label zoneID = faceZones.findZoneID(name);
+
+    if
+    (
+        (r == IOobject::MUST_READ)
+     || (r == IOobject::READ_IF_PRESENT && zoneID != -1)
+    )
+    {
+        const faceZone& fz = faceZones[zoneID];
+        addressing_ = fz;
+        flipMap_ = fz.flipMap();
+    }
+
+    updateSet();
+
+    check(mesh.nFaces());
+}
+
+
+faceZoneSet::faceZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const label size,
+    writeOption w
+)
+:
+    faceSet(mesh, name, size, w),
+    mesh_(mesh),
+    addressing_(0),
+    flipMap_(0)
+{
+    updateSet();
+}
+
+
+faceZoneSet::faceZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const topoSet& set,
+    writeOption w
+)
+:
+    faceSet(mesh, name, set.size(), w),
+    mesh_(mesh),
+    addressing_(refCast<const faceZoneSet>(set).addressing()),
+    flipMap_(refCast<const faceZoneSet>(set).flipMap())
+{
+    updateSet();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+faceZoneSet::~faceZoneSet()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void faceZoneSet::invert(const label maxLen)
+{
+    label n = 0;
+
+    for (label faceI = 0; faceI < maxLen; faceI++)
+    {
+        if (!found(faceI))
+        {
+            addressing_[n] = faceI;
+            flipMap_[n] = false;         //? or true?
+            n++;
+        }
+    }
+    addressing_.setSize(n);
+    flipMap_.setSize(n);
+    updateSet();
+}
+
+
+void faceZoneSet::subset(const topoSet& set)
+{
+    label nConflict = 0;
+
+    DynamicList<label> newAddressing(addressing_.size());
+    DynamicList<bool> newFlipMap(flipMap_.size());
+
+    Map<label> faceToIndex(addressing_.size());
+    forAll(addressing_, i)
+    {
+        faceToIndex.insert(addressing_[i], i);
+    }
+
+    const faceZoneSet& fSet = refCast<const faceZoneSet>(set);
+
+    forAll(fSet.addressing(), i)
+    {
+        label faceI = fSet.addressing()[i];
+
+        Map<label>::const_iterator iter = faceToIndex.find(faceI);
+
+        if (iter != faceToIndex.end())
+        {
+            label index = iter();
+
+            if (fSet.flipMap()[i] != flipMap_[index])
+            {
+                nConflict++;
+            }
+            newAddressing.append(faceI);
+            newFlipMap.append(flipMap_[index]);
+        }
+    }
+
+    if (nConflict > 0)
+    {
+        WarningIn(" faceZoneSet::subset(const topoSet&)")
+            << "subset : there are " << nConflict
+            << " faces with different orientation in faceZonesSets "
+            << name() << " and " << set.name() << endl;
+    }
+
+    addressing_.transfer(newAddressing);
+    flipMap_.transfer(newFlipMap);
+    updateSet();
+}
+
+
+void faceZoneSet::addSet(const topoSet& set)
+{
+    label nConflict = 0;
+
+    DynamicList<label> newAddressing(addressing_);
+    DynamicList<bool> newFlipMap(flipMap_);
+
+    Map<label> faceToIndex(addressing_.size());
+    forAll(addressing_, i)
+    {
+        faceToIndex.insert(addressing_[i], i);
+    }
+
+    const faceZoneSet& fSet = refCast<const faceZoneSet>(set);
+
+    forAll(fSet.addressing(), i)
+    {
+        label faceI = fSet.addressing()[i];
+
+        Map<label>::const_iterator iter = faceToIndex.find(faceI);
+
+        if (iter != faceToIndex.end())
+        {
+            label index = iter();
+
+            if (fSet.flipMap()[i] != flipMap_[index])
+            {
+                nConflict++;
+            }
+        }
+        else
+        {
+            newAddressing.append(faceI);
+            newFlipMap.append(fSet.flipMap()[i]);
+        }
+    }
+
+    if (nConflict > 0)
+    {
+        WarningIn("faceZoneSet::addSet(const topoSet&)")
+            << "addSet : there are " << nConflict
+            << " faces with different orientation in faceZonesSets "
+            << name() << " and " << set.name() << endl;
+    }
+
+    addressing_.transfer(newAddressing);
+    flipMap_.transfer(newFlipMap);
+    updateSet();
+}
+
+
+void faceZoneSet::deleteSet(const topoSet& set)
+{
+    label nConflict = 0;
+
+    DynamicList<label> newAddressing(addressing_.size());
+    DynamicList<bool> newFlipMap(flipMap_.size());
+
+    const faceZoneSet& fSet = refCast<const faceZoneSet>(set);
+
+    Map<label> faceToIndex(fSet.addressing().size());
+    forAll(fSet.addressing(), i)
+    {
+        faceToIndex.insert(fSet.addressing()[i], i);
+    }
+
+    forAll(addressing_, i)
+    {
+        label faceI = addressing_[i];
+
+        Map<label>::const_iterator iter = faceToIndex.find(faceI);
+
+        if (iter != faceToIndex.end())
+        {
+            label index = iter();
+
+            if (fSet.flipMap()[index] != flipMap_[i])
+            {
+                nConflict++;
+            }
+        }
+        else
+        {
+            // Not found in fSet so add
+            newAddressing.append(faceI);
+            newFlipMap.append(fSet.flipMap()[i]);
+        }
+    }
+
+    if (nConflict > 0)
+    {
+        WarningIn("faceZoneSet::deleteSet(const topoSet&)")
+            << "deleteSet : there are " << nConflict
+            << " faces with different orientation in faceZonesSets "
+            << name() << " and " << set.name() << endl;
+    }
+
+    addressing_.transfer(newAddressing);
+    flipMap_.transfer(newFlipMap);
+    updateSet();
+}
+
+
+void faceZoneSet::sync(const polyMesh& mesh)
+{}
+
+
+label faceZoneSet::maxSize(const polyMesh& mesh) const
+{
+    return mesh.nFaces();
+}
+
+
+//- Write using given format, version and compression
+bool faceZoneSet::writeObject
+(
+    IOstream::streamFormat s,
+    IOstream::versionNumber v,
+    IOstream::compressionType c
+) const
+{
+    // Write shadow faceSet
+    word oldTypeName = typeName;
+    const_cast<word&>(type()) = faceSet::typeName;
+    bool ok = faceSet::writeObject(s, v, c);
+    const_cast<word&>(type()) = oldTypeName;
+
+    // Modify faceZone
+    faceZoneMesh& faceZones = const_cast<polyMesh&>(mesh_).faceZones();
+    label zoneID = faceZones.findZoneID(name());
+
+    if (zoneID == -1)
+    {
+        zoneID = faceZones.size();
+
+        faceZones.setSize(zoneID+1);
+        faceZones.set
+        (
+            zoneID,
+            new faceZone
+            (
+                name(),
+                addressing_,
+                flipMap_,
+                zoneID,
+                faceZones
+            )
+        );
+    }
+    else
+    {
+        faceZones[zoneID].resetAddressing(addressing_, flipMap_);
+    }
+    faceZones.clearAddressing();
+
+    return ok && faceZones.write();
+}
+
+
+void faceZoneSet::updateMesh(const mapPolyMesh& morphMap)
+{
+    // faceZone
+    labelList newAddressing(addressing_.size());
+    boolList newFlipMap(flipMap_.size());
+
+    label n = 0;
+    forAll(addressing_, i)
+    {
+        label faceI = addressing_[i];
+        label newFaceI = morphMap.reverseFaceMap()[faceI];
+        if (newFaceI >= 0)
+        {
+            newAddressing[n] = newFaceI;
+            newFlipMap[n] = flipMap_[i];
+            n++;
+        }
+    }
+    newAddressing.setSize(n);
+    newFlipMap.setSize(n);
+
+    addressing_.transfer(newAddressing);
+    flipMap_.transfer(newFlipMap);
+
+    updateSet();
+}
+
+
+void faceZoneSet::writeDebug
+(
+    Ostream& os,
+    const primitiveMesh& mesh,
+    const label maxLen
+) const
+{
+    faceSet::writeDebug(os, mesh, maxLen);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/meshTools/sets/topoSets/faceZoneSet.H

@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::faceZoneSet
+
+Description
+    Like faceSet but updates faceZone when writing.
+
+SourceFiles
+    faceZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faceZoneSet_H
+#define faceZoneSet_H
+
+#include "faceSet.H"
+#include "boolList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class faceZoneSet Declaration
+\*---------------------------------------------------------------------------*/
+
+class faceZoneSet
+:
+    public faceSet
+{
+    // Private data
+
+        const polyMesh& mesh_;
+
+        labelList addressing_;
+
+        boolList flipMap_;
+
+   // Private Member Functions
+
+
+public:
+
+    //- Runtime type information
+    TypeName("faceZoneSet");
+
+
+    // Constructors
+
+        //- Construct from objectRegistry and name
+        faceZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            readOption r=MUST_READ,
+            writeOption w=NO_WRITE
+        );
+
+        //- Construct from additional size of labelHashSet
+        faceZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const label,
+            writeOption w=NO_WRITE
+        );
+
+        //- Construct from existing set
+        faceZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const topoSet&,
+            writeOption w=NO_WRITE
+        );
+
+
+
+    // Destructor
+
+        virtual ~faceZoneSet();
+
+
+    // Member functions
+
+        const labelList& addressing() const
+        {
+            return addressing_;
+        }
+
+        labelList& addressing()
+        {
+            return addressing_;
+        }
+
+
+        const boolList& flipMap() const
+        {
+            return flipMap_;
+        }
+
+        boolList& flipMap()
+        {
+            return flipMap_;
+        }
+
+
+        //- Sort addressing and make faceSet part consistent with addressing
+        void updateSet();
+
+        //- Invert contents. (insert all members 0..maxLen-1 which were not in
+        //  set)
+        virtual void invert(const label maxLen);
+
+        //- Subset contents. Only elements present in both sets remain.
+        virtual void subset(const topoSet& set);
+
+        //- Add elements present in set.
+        virtual void addSet(const topoSet& set);
+
+        //- Delete elements present in set.
+        virtual void deleteSet(const topoSet& set);
+
+        //- Sync faceZoneSet across coupled patches.
+        virtual void sync(const polyMesh& mesh);
+
+        //- Write maxLen items with label and coordinates. 
+        virtual void writeDebug
+        (
+            Ostream& os,
+            const primitiveMesh&,
+            const label maxLen
+        ) const;
+
+        //- Write faceZone
+        virtual bool writeObject
+        (
+            IOstream::streamFormat,
+            IOstream::versionNumber,
+            IOstream::compressionType
+        ) const;
+
+        //- Update any stored data for new labels
+        virtual void updateMesh(const mapPolyMesh& morphMap);
+
+        //- Return max index+1.
+        virtual label maxSize(const polyMesh& mesh) const;
+
+
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/sets/topoSets/pointSet.C

@@ -42,6 +42,7 @@ defineTypeNameAndDebug(pointSet, 0);
 
 addToRunTimeSelectionTable(topoSet, pointSet, word);
 addToRunTimeSelectionTable(topoSet, pointSet, size);
+addToRunTimeSelectionTable(topoSet, pointSet, set);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -78,6 +79,18 @@ pointSet::pointSet
 {}
 
 
+pointSet::pointSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const topoSet& set,
+    writeOption w
+)
+:
+    topoSet(mesh, name, set, w)
+{}
+
+
 pointSet::pointSet
 (
     const polyMesh& mesh,

src/meshTools/sets/topoSets/pointSet.H

@@ -83,6 +83,15 @@ public:
             writeOption w=NO_WRITE
         );
 
+        //- Construct from additional labelHashSet
+        pointSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const topoSet&,
+            writeOption w=NO_WRITE
+        );
+
         //- Construct from additional labelHashSet
         pointSet
         (

src/meshTools/sets/topoSets/pointZoneSet.C

@@ -0,0 +1,314 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "pointZoneSet.H"
+#include "mapPolyMesh.H"
+#include "polyMesh.H"
+#include "processorPolyPatch.H"
+#include "cyclicPolyPatch.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(pointZoneSet, 0);
+
+addToRunTimeSelectionTable(topoSet, pointZoneSet, word);
+addToRunTimeSelectionTable(topoSet, pointZoneSet, size);
+addToRunTimeSelectionTable(topoSet, pointZoneSet, set);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void pointZoneSet::updateSet()
+{
+    labelList order;
+    sortedOrder(addressing_, order);
+    inplaceReorder(order, addressing_);
+
+    pointSet::clearStorage();
+    pointSet::resize(2*addressing_.size());
+    forAll(addressing_, i)
+    {
+        pointSet::insert(addressing_[i]);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+pointZoneSet::pointZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    readOption r,
+    writeOption w
+)
+:
+    pointSet(mesh, name, 1000),  // do not read pointSet
+    mesh_(mesh),
+    addressing_(0)
+{
+    const pointZoneMesh& pointZones = mesh.pointZones();
+    label zoneID = pointZones.findZoneID(name);
+
+    if
+    (
+        (r == IOobject::MUST_READ)
+     || (r == IOobject::READ_IF_PRESENT && zoneID != -1)
+    )
+    {
+        const pointZone& fz = pointZones[zoneID];
+        addressing_ = fz;
+    }
+
+    updateSet();
+
+    check(mesh.nPoints());
+}
+
+
+pointZoneSet::pointZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const label size,
+    writeOption w
+)
+:
+    pointSet(mesh, name, size, w),
+    mesh_(mesh),
+    addressing_(0)
+{
+    updateSet();
+}
+
+
+pointZoneSet::pointZoneSet
+(
+    const polyMesh& mesh,
+    const word& name,
+    const topoSet& set,
+    writeOption w
+)
+:
+    pointSet(mesh, name, set.size(), w),
+    mesh_(mesh),
+    addressing_(refCast<const pointZoneSet>(set).addressing())
+{
+    updateSet();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+pointZoneSet::~pointZoneSet()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void pointZoneSet::invert(const label maxLen)
+{
+    label n = 0;
+
+    for (label pointI = 0; pointI < maxLen; pointI++)
+    {
+        if (!found(pointI))
+        {
+            addressing_[n] = pointI;
+            n++;
+        }
+    }
+    addressing_.setSize(n);
+    updateSet();
+}
+
+
+void pointZoneSet::subset(const topoSet& set)
+{
+    DynamicList<label> newAddressing(addressing_.size());
+
+    const pointZoneSet& fSet = refCast<const pointZoneSet>(set);
+
+    forAll(fSet.addressing(), i)
+    {
+        label pointI = fSet.addressing()[i];
+
+        if (found(pointI))
+        {
+            newAddressing.append(pointI);
+        }
+    }
+
+    addressing_.transfer(newAddressing);
+    updateSet();
+}
+
+
+void pointZoneSet::addSet(const topoSet& set)
+{
+    DynamicList<label> newAddressing(addressing_);
+
+    const pointZoneSet& fSet = refCast<const pointZoneSet>(set);
+
+    forAll(fSet.addressing(), i)
+    {
+        label pointI = fSet.addressing()[i];
+
+        if (!found(pointI))
+        {
+            newAddressing.append(pointI);
+        }
+    }
+
+    addressing_.transfer(newAddressing);
+    updateSet();
+}
+
+
+void pointZoneSet::deleteSet(const topoSet& set)
+{
+    DynamicList<label> newAddressing(addressing_.size());
+
+    const pointZoneSet& fSet = refCast<const pointZoneSet>(set);
+
+    forAll(addressing_, i)
+    {
+        label pointI = addressing_[i];
+
+        if (!fSet.found(pointI))
+        {
+            // Not found in fSet so add
+            newAddressing.append(pointI);
+        }
+    }
+
+    addressing_.transfer(newAddressing);
+    updateSet();
+}
+
+
+void pointZoneSet::sync(const polyMesh& mesh)
+{}
+
+
+label pointZoneSet::maxSize(const polyMesh& mesh) const
+{
+    return mesh.nPoints();
+}
+
+
+//- Write using given format, version and compression
+bool pointZoneSet::writeObject
+(
+    IOstream::streamFormat s,
+    IOstream::versionNumber v,
+    IOstream::compressionType c
+) const
+{
+    // Write shadow pointSet
+    word oldTypeName = typeName;
+    const_cast<word&>(type()) = pointSet::typeName;
+    bool ok = pointSet::writeObject(s, v, c);
+    const_cast<word&>(type()) = oldTypeName;
+
+    // Modify pointZone
+    pointZoneMesh& pointZones = const_cast<polyMesh&>(mesh_).pointZones();
+    label zoneID = pointZones.findZoneID(name());
+
+    if (zoneID == -1)
+    {
+        zoneID = pointZones.size();
+
+        pointZones.setSize(zoneID+1);
+        pointZones.set
+        (
+            zoneID,
+            new pointZone
+            (
+                name(),
+                addressing_,
+                zoneID,
+                pointZones
+            )
+        );
+    }
+    else
+    {
+        pointZones[zoneID] = addressing_;
+    }
+    pointZones.clearAddressing();
+
+    return ok && pointZones.write();
+}
+
+
+void pointZoneSet::updateMesh(const mapPolyMesh& morphMap)
+{
+    // pointZone
+    labelList newAddressing(addressing_.size());
+
+    label n = 0;
+    forAll(addressing_, i)
+    {
+        label pointI = addressing_[i];
+        label newPointI = morphMap.reversePointMap()[pointI];
+        if (newPointI >= 0)
+        {
+            newAddressing[n] = newPointI;
+            n++;
+        }
+    }
+    newAddressing.setSize(n);
+
+    addressing_.transfer(newAddressing);
+
+    updateSet();
+}
+
+
+void pointZoneSet::writeDebug
+(
+    Ostream& os,
+    const primitiveMesh& mesh,
+    const label maxLen
+) const
+{
+    pointSet::writeDebug(os, mesh, maxLen);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/meshTools/sets/topoSets/pointZoneSet.H

@@ -0,0 +1,173 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::pointZoneSet
+
+Description
+    Like pointSet but updates pointZone when writing.
+
+SourceFiles
+    pointZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pointZoneSet_H
+#define pointZoneSet_H
+
+#include "pointSet.H"
+#include "boolList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class pointZoneSet Declaration
+\*---------------------------------------------------------------------------*/
+
+class pointZoneSet
+:
+    public pointSet
+{
+    // Private data
+
+        const polyMesh& mesh_;
+
+        labelList addressing_;
+
+   // Private Member Functions
+
+
+public:
+
+    //- Runtime type information
+    TypeName("pointZoneSet");
+
+
+    // Constructors
+
+        //- Construct from objectRegistry and name
+        pointZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            readOption r=MUST_READ,
+            writeOption w=NO_WRITE
+        );
+
+        //- Construct from additional size of labelHashSet
+        pointZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const label,
+            writeOption w=NO_WRITE
+        );
+
+        //- Construct from existing set
+        pointZoneSet
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const topoSet&,
+            writeOption w=NO_WRITE
+        );
+
+
+
+    // Destructor
+
+        virtual ~pointZoneSet();
+
+
+    // Member functions
+
+        const labelList& addressing() const
+        {
+            return addressing_;
+        }
+
+        labelList& addressing()
+        {
+            return addressing_;
+        }
+
+        //- Sort addressing and make pointSet part consistent with addressing
+        void updateSet();
+
+        //- Invert contents. (insert all members 0..maxLen-1 which were not in
+        //  set)
+        virtual void invert(const label maxLen);
+
+        //- Subset contents. Only elements present in both sets remain.
+        virtual void subset(const topoSet& set);
+
+        //- Add elements present in set.
+        virtual void addSet(const topoSet& set);
+
+        //- Delete elements present in set.
+        virtual void deleteSet(const topoSet& set);
+
+        //- Sync pointZoneSet across coupled patches.
+        virtual void sync(const polyMesh& mesh);
+
+        //- Write maxLen items with label and coordinates. 
+        virtual void writeDebug
+        (
+            Ostream& os,
+            const primitiveMesh&,
+            const label maxLen
+        ) const;
+
+        //- Write pointZone
+        virtual bool writeObject
+        (
+            IOstream::streamFormat,
+            IOstream::versionNumber,
+            IOstream::compressionType
+        ) const;
+
+        //- Update any stored data for new labels
+        virtual void updateMesh(const mapPolyMesh& morphMap);
+
+        //- Return max index+1.
+        virtual label maxSize(const polyMesh& mesh) const;
+
+
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/sets/topoSets/topoSet.C

@@ -39,6 +39,7 @@ namespace Foam
 defineTypeNameAndDebug(topoSet, 0);
 defineRunTimeSelectionTable(topoSet, word);
 defineRunTimeSelectionTable(topoSet, size);
+defineRunTimeSelectionTable(topoSet, set);
 
 
 // Construct named object from existing set.
@@ -103,6 +104,37 @@ autoPtr<topoSet> topoSet::New
 }
 
 
+// Construct named object from existing set.
+autoPtr<topoSet> topoSet::New
+(
+    const word& setType,
+    const polyMesh& mesh,
+    const word& name,
+    const topoSet& set,
+    writeOption w
+)
+{
+    setConstructorTable::iterator cstrIter =
+        setConstructorTablePtr_
+            ->find(setType);
+
+    if (cstrIter == setConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "topoSet::New(const word&, "
+            "const polyMesh&, const word&, const topoSet&, writeOption)"
+        )   << "Unknown set type " << setType
+            << endl << endl
+            << "Valid set types : " << endl
+            << setConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<topoSet>(cstrIter()(mesh, name, set, w));
+}
+
+
 Foam::fileName topoSet::topoSet::localPath
 (
     const polyMesh& mesh,
@@ -255,9 +287,9 @@ void topoSet::writeDebug
 ) const
 {
     // Bounding box of contents.
-    boundBox bb(pointField(coords, toc()));
+    boundBox bb(pointField(coords, toc()), true);
 
-    Pout<< "Set bounding box: min = "
+    os  << "Set bounding box: min = "
         << bb.min() << "    max = " << bb.max() << " meters. " << endl << endl;
 
     label n = 0;
@@ -538,18 +570,18 @@ void topoSet::writeDebug(Ostream& os, const label maxLen) const
 }
 
 
-void topoSet::writeDebug
+//void topoSet::writeDebug
-(
+//(
-    Ostream&,
+//    Ostream&,
-    const primitiveMesh&,
+//    const primitiveMesh&,
-    const label
+//    const label
-) const
+//) const
-{
+//{
-    notImplemented
+//    notImplemented
-    (
+//    (
-        "topoSet::writeDebug(Ostream&, const primitiveMesh&, const label)"
+//        "topoSet::writeDebug(Ostream&, const primitiveMesh&, const label)"
-    );
+//    );
-}
+//}
 
 
 bool topoSet::writeData(Ostream& os) const
@@ -564,13 +596,13 @@ void topoSet::updateMesh(const mapPolyMesh&)
 }
 
 
-//- Return max index+1.
+////- Return max index+1.
-label topoSet::maxSize(const polyMesh&) const
+//label topoSet::maxSize(const polyMesh&) const
-{
+//{
-    notImplemented("topoSet::maxSize(const polyMesh&)");
+//    notImplemented("topoSet::maxSize(const polyMesh&)");
-
+//
-    return -1;
+//    return -1;
-}
+//}
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 

src/meshTools/sets/topoSets/topoSet.H

@@ -55,7 +55,6 @@ namespace Foam
 class mapPolyMesh;
 class polyMesh;
 class primitiveMesh;
-//class directPolyTopoChange;
 
 /*---------------------------------------------------------------------------*\
                            Class topoSet Declaration
@@ -154,6 +153,21 @@ public:
             (mesh, name, size, w)
         );
 
+        // For the constructor as copy
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            topoSet,
+            set,
+            (
+                const polyMesh& mesh,
+                const word& name,
+                const topoSet& set,
+                writeOption w
+            ),
+            (mesh, name, set, w)
+        );
+
 
     // Constructors
 
@@ -208,7 +222,7 @@ public:
 
     // Selectors
 
-        //- Return a reference to the selected source
+        //- Return a pointer to a toposet read from file
         static autoPtr<topoSet> New
         (
             const word& setType,
@@ -218,7 +232,7 @@ public:
             writeOption w=NO_WRITE
         );
 
-        //- Return a reference to the selected source
+        //- Return a pointer to a new toposet of given size
         static autoPtr<topoSet> New
         (
             const word& setType,
@@ -228,6 +242,16 @@ public:
             writeOption w=NO_WRITE
         );
 
+        //- Return a pointer to a new toposet as copy of another toposet
+        static autoPtr<topoSet> New
+        (
+            const word& setType,
+            const polyMesh& mesh,
+            const word& name,
+            const topoSet& set,
+            writeOption w=NO_WRITE
+        );
+
 
     // Destructor
 
@@ -256,13 +280,13 @@ public:
         virtual void writeDebug(Ostream& os, const label maxLen) const;
 
         //- Like above but also writes mesh related quantity
-        //  (usually coordinate). Not implemented.
+        //  (usually coordinate).
         virtual void writeDebug
         (
             Ostream& os,
             const primitiveMesh&,
             const label maxLen
-        ) const;
+        ) const = 0;
 
         //- Write contents.
         virtual bool writeData(Ostream&) const;
@@ -271,7 +295,7 @@ public:
         virtual void updateMesh(const mapPolyMesh& morphMap);
 
         //- Return max allowable index (+1). Not implemented.
-        virtual label maxSize(const polyMesh& mesh) const;
+        virtual label maxSize(const polyMesh& mesh) const = 0;
 
 
 

src/thermophysicalModels/radiation/radiationModel/fvDOM/absorptionCoeffs/absorptionCoeffs.C

@@ -99,7 +99,7 @@ void Foam::radiation::absorptionCoeffs::initialise(const dictionary& dict)
 {
     dict.lookup("Tcommon") >> Tcommon_;
     dict.lookup("Tlow") >> Tlow_;
-    dict.lookup("Tlow") >> Thigh_;
+    dict.lookup("Thigh") >> Thigh_;
     dict.lookup("invTemp") >> invTemp_;
 
     dict.lookup("loTcoeffs") >> lowACoeffs_;

src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.C

@@ -26,7 +26,6 @@ License
 
 #include "fvDOM.H"
 #include "addToRunTimeSelectionTable.H"
-#include "fvm.H"
 
 #include "absorptionEmissionModel.H"
 #include "scatterModel.H"
@@ -65,7 +64,7 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
         mesh_,
         dimensionedScalar("G", dimMass/pow3(dimTime), 0.0)
@@ -347,7 +346,7 @@ Foam::radiation::fvDOM::Ru() const
     const DimensionedField<scalar, volMesh> a =
         a_.dimensionedInternalField(); //absorptionEmission_->aCont()()
 
-    return  a*G - 4.0*E;
+    return  a*G - E;
 }
 
 

src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C

@@ -27,12 +27,9 @@ License
 #include "radiativeIntensityRay.H"
 #include "fvm.H"
 #include "fvDOM.H"
+#include "constants.H"
 
-#include "absorptionEmissionModel.H"
+using namespace Foam::constant;
-#include "scatterModel.H"
-#include "mathConstants.H"
-#include "radiationModel.H"
-#include "Vector2D.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -202,9 +199,9 @@ Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
             fvm::div(Ji, ILambda_[lambdaI], "div(Ji,Ii_h)")
           + fvm::Sp(k*omega_, ILambda_[lambdaI])
          ==
-            1.0/constant::math::pi
+            1.0/constant::math::pi*omega_
            *(
-                k*omega_*blackBody_.bLambda(lambdaI)
+                k*blackBody_.bLambda(lambdaI)
               + absorptionEmission_.ECont(lambdaI)
             )
         );

src/thermophysicalModels/radiation/radiationModel/radiationModel/newRadiationModel.C

@@ -43,19 +43,19 @@ autoPtr<radiationModel> radiationModel::New
 {
     word radiationModelTypeName;
 
-    // Enclose the creation of the radiationPropertiesDict to ensure it is
+    // Note: no need to register/keep radiationProperties since models read
-    // deleted before the radiationModel is created otherwise the dictionary
+    // it themselves.
-    // is entered in the database twice
     {
         IOdictionary radiationPropertiesDict
         (
             IOobject
             (
                 "radiationProperties",
-                T.mesh().time().constant(),
+                T.time().constant(),
-                T.mesh().objectRegistry::db(),
+                T.mesh(),
                 IOobject::MUST_READ,
-                IOobject::NO_WRITE
+                IOobject::NO_WRITE,
+                false
             )
         );
 

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C

@@ -259,10 +259,11 @@ Foam::radiation::greyMeanAbsorptionEmission::ECont(const label bandI) const
         )
     );
 
-    if (mesh_.foundObject<volScalarField>("hrr"))
+    if (mesh_.foundObject<volScalarField>("dQ"))
     {
-        const volScalarField& hrr = mesh_.lookupObject<volScalarField>("hrr");
+        const volScalarField& dQ =
-        E().internalField() = EhrrCoeff_*hrr.internalField();
+            mesh_.lookupObject<volScalarField>("dQ");
+        E().internalField() = EhrrCoeff_*dQ;
     }
 
     return E;

src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffle/turbulentTemperatureCoupledBaffleFvPatchScalarField.C

@@ -215,6 +215,8 @@ turbulentTemperatureCoupledBaffleFvPatchScalarField
 tmp<scalarField>
 turbulentTemperatureCoupledBaffleFvPatchScalarField::K() const
 {
+    const fvMesh& mesh = patch().boundaryMesh().mesh();
+
     if (KName_ == "none")
     {
         const compressible::RASModel& model =
@@ -227,12 +229,33 @@ turbulentTemperatureCoupledBaffleFvPatchScalarField::K() const
 
         return
             talpha().boundaryField()[patch().index()]
-           *thermo.rho()().boundaryField()[patch().index()]
            *thermo.Cp()().boundaryField()[patch().index()];
     }
+    else if (mesh.objectRegistry::foundObject<volScalarField>(KName_))
+    {
+        return patch().lookupPatchField<volScalarField, scalar>(KName_);
+    }
+    else if (mesh.objectRegistry::foundObject<volSymmTensorField>(KName_))
+    {
+        const symmTensorField& KWall =
+            patch().lookupPatchField<volSymmTensorField, scalar>(KName_);
+
+        vectorField n = patch().nf();
+
+        return n & KWall & n;
+    }
     else
     {
-        return patch().lookupPatchField<volScalarField, scalar>(KName_);
+        FatalErrorIn
+        (
+            "turbulentTemperatureCoupledBaffleFvPatchScalarField::K() const"
+        )   << "Did not find field " << KName_
+            << " on mesh " << mesh.name() << " patch " << patch().name()
+            << endl
+            << "Please set 'K' to 'none', a valid volScalarField"
+            << " or a valid volSymmTensorField." << exit(FatalError);
+
+        return scalarField(0);
     }
 }
 

src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutkRoughWallFunction/mutkRoughWallFunctionFvPatchScalarField.C

@@ -77,7 +77,7 @@ tmp<scalarField> mutkRoughWallFunctionFvPatchScalarField::calcMut() const
 
     const scalar Cmu25 = pow025(Cmu_);
 
-    tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
+    tmp<scalarField> tmutw(new scalarField(*this));
     scalarField& mutw = tmutw();
 
     forAll(mutw, faceI)

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.C

@@ -76,7 +76,7 @@ tmp<scalarField> nutkRoughWallFunctionFvPatchScalarField::calcNut() const
 
     const scalar Cmu25 = pow025(Cmu_);
 
-    tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
+    tmp<scalarField> tnutw(new scalarField(*this));
     scalarField& nutw = tnutw();
 
     forAll(nutw, faceI)