diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
index a7a65b17bd4c9863a23d0749d64e3cbcb5af5348..6e87fb92034ea9ec69234ae97651637117795afd 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
@@ -2,7 +2,9 @@ scalar DiNum = -GREAT;
 
 forAll(solidRegions, i)
 {
-    #include "setRegionSolidFields.H"
+    //- Note: do not use setRegionSolidFields.H to avoid double registering Cp
+    //#include "setRegionSolidFields.H"
+    const solidThermo& thermo = thermos[i];
 
     tmp<volScalarField> magKappa;
     if (thermo.isotropic())
@@ -14,6 +16,12 @@ forAll(solidRegions, i)
         magKappa = mag(thermo.Kappa());
     }
 
+    tmp<volScalarField> tcp = thermo.Cp();
+    const volScalarField& cp = tcp();
+
+    tmp<volScalarField> trho = thermo.rho();
+    const volScalarField& rho = trho();
+
     DiNum = max
     (
         solidRegionDiffNo
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index c349740ea7185ca30e4587003a07c9925909720c..b0cf1879c1e43567ed9a80e8d9c9111bb68ff73a 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -576,25 +576,30 @@ void writeMesh
 (
     const string& msg,
     const meshRefinement& meshRefiner,
-    const bool writeLevel,
-    const label debug
+    const meshRefinement::debugType debugLevel,
+    const meshRefinement::writeType writeLevel
 )
 {
     const fvMesh& mesh = meshRefiner.mesh();
 
-    meshRefiner.printMeshInfo(debug, msg);
+    meshRefiner.printMeshInfo(debugLevel, msg);
     Info<< "Writing mesh to time " << meshRefiner.timeName() << endl;
 
-    label flag = meshRefinement::MESH;
-    if (writeLevel)
-    {
-        flag |= meshRefinement::SCALARLEVELS;
-    }
-    if (debug & meshRefinement::OBJINTERSECTIONS)
-    {
-        flag |= meshRefinement::OBJINTERSECTIONS;
-    }
-    meshRefiner.write(flag, mesh.time().path()/meshRefiner.timeName());
+    //label flag = meshRefinement::MESH;
+    //if (writeLevel)
+    //{
+    //    flag |= meshRefinement::SCALARLEVELS;
+    //}
+    //if (debug & meshRefinement::OBJINTERSECTIONS)
+    //{
+    //    flag |= meshRefinement::OBJINTERSECTIONS;
+    //}
+    meshRefiner.write
+    (
+        debugLevel,
+        meshRefinement::writeType(writeLevel | meshRefinement::WRITEMESH),
+        mesh.time().path()/meshRefiner.timeName()
+    );
     Info<< "Wrote mesh in = "
         << mesh.time().cpuTimeIncrement() << " s." << endl;
 }
@@ -837,16 +842,74 @@ int main(int argc, char *argv[])
     // Debug
     // ~~~~~
 
-    const label debug = meshDict.lookupOrDefault<label>("debug", 0);
-    if (debug > 0)
+    // Set debug level
+    meshRefinement::debugType debugLevel = meshRefinement::debugType
+    (
+        meshDict.lookupOrDefault<label>
+        (
+            "debug",
+            0
+        )
+    );
     {
-        meshRefinement::debug   = debug;
-        autoRefineDriver::debug = debug;
-        autoSnapDriver::debug   = debug;
-        autoLayerDriver::debug  = debug;
+        wordList flags;
+        if (meshDict.readIfPresent("debugFlags", flags))
+        {
+            debugLevel = meshRefinement::debugType
+            (
+                meshRefinement::readFlags
+                (
+                    meshRefinement::IOdebugTypeNames,
+                    flags
+                )
+            );
+        }
+    }
+    if (debugLevel > 0)
+    {
+        meshRefinement::debug   = debugLevel;
+        autoRefineDriver::debug = debugLevel;
+        autoSnapDriver::debug   = debugLevel;
+        autoLayerDriver::debug  = debugLevel;
+    }
+
+    // Set file writing level
+    {
+        wordList flags;
+        if (meshDict.readIfPresent("writeFlags", flags))
+        {
+            meshRefinement::writeLevel
+            (
+                meshRefinement::writeType
+                (
+                    meshRefinement::readFlags
+                    (
+                        meshRefinement::IOwriteTypeNames,
+                        flags
+                    )
+                )
+            );
+        }
     }
 
-    const bool writeLevel = meshDict.lookupOrDefault<bool>("writeLevel", false);
+    // Set output level
+    {
+        wordList flags;
+        if (meshDict.readIfPresent("outputFlags", flags))
+        {
+            meshRefinement::outputLevel
+            (
+                meshRefinement::outputType
+                (
+                    meshRefinement::readFlags
+                    (
+                        meshRefinement::IOoutputTypeNames,
+                        flags
+                    )
+                )
+            );
+        }
+    }
 
 
     // Read geometry
@@ -1047,11 +1110,12 @@ int main(int argc, char *argv[])
         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
 
     // Some stats
-    meshRefiner.printMeshInfo(debug, "Initial mesh");
+    meshRefiner.printMeshInfo(debugLevel, "Initial mesh");
 
     meshRefiner.write
     (
-        debug & meshRefinement::OBJINTERSECTIONS,
+        meshRefinement::debugType(debugLevel&meshRefinement::OBJINTERSECTIONS),
+        meshRefinement::writeType(0),
         mesh.time().path()/meshRefiner.timeName()
     );
 
@@ -1271,7 +1335,7 @@ int main(int argc, char *argv[])
         );
 
 
-        if (!overwrite && !debug)
+        if (!overwrite && !debugLevel)
         {
             const_cast<Time&>(mesh.time())++;
         }
@@ -1289,8 +1353,8 @@ int main(int argc, char *argv[])
         (
             "Refined mesh",
             meshRefiner,
-            writeLevel,
-            debug
+            debugLevel,
+            meshRefinement::writeLevel()
         );
 
         Info<< "Mesh refined in = "
@@ -1308,7 +1372,7 @@ int main(int argc, char *argv[])
             globalToSlavePatch
         );
 
-        if (!overwrite && !debug)
+        if (!overwrite && !debugLevel)
         {
             const_cast<Time&>(mesh.time())++;
         }
@@ -1330,8 +1394,8 @@ int main(int argc, char *argv[])
         (
             "Snapped mesh",
             meshRefiner,
-            writeLevel,
-            debug
+            debugLevel,
+            meshRefinement::writeLevel()
         );
 
         Info<< "Mesh snapped in = "
@@ -1357,7 +1421,7 @@ int main(int argc, char *argv[])
         );
 
 
-        if (!overwrite &&  !debug)
+        if (!overwrite && !debugLevel)
         {
             const_cast<Time&>(mesh.time())++;
         }
@@ -1376,8 +1440,8 @@ int main(int argc, char *argv[])
         (
             "Layer mesh",
             meshRefiner,
-            writeLevel,
-            debug
+            debugLevel,
+            meshRefinement::writeLevel()
         );
 
         Info<< "Layers added in = "
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index b850033a44db51f27f5eb0882b33e453642d5361..507891d356f316fdee9d2a2227c86c19fc49f1aa 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -475,14 +475,22 @@ meshQualityControls
 
 // Advanced
 
-// Flags for optional output
-//  0 : only write final meshes
-//  1 : write intermediate meshes
-//  2 : write volScalarField with cellLevel for postprocessing
-//  4 : write current mesh intersections as .obj files
-//  8 : write information about explicit feature edge refinement
-// 16 : write information about layers
-debug 0;
+//// Debug flags
+//debugFlags
+//(
+//    mesh            // write intermediate meshes
+//    intersections   // write current mesh intersections as .obj files
+//    featureSeeds,   // write information about explicit feature edge refinement
+//    layerInfo       // write information about layers
+//);
+//
+//// Write flags
+//writeFlags
+//(
+//    scalarLevels    // write volScalarField with cellLevel for postprocessing
+//    layerSets       // write cellSets, faceSets of faces in layer
+//    layerFields     // write volScalarField for layer coverage
+//);
 
 // Merge tolerance. Is fraction of overall bounding box of initial mesh.
 // Note: the write tolerance needs to be higher than this.
diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/orientFaceZone.C b/applications/utilities/mesh/manipulation/orientFaceZone/orientFaceZone.C
index 5f1a985116d2d967e3b289f1a495617ec680473a..45458315ab0261b30d4be313267538140820304b 100644
--- a/applications/utilities/mesh/manipulation/orientFaceZone/orientFaceZone.C
+++ b/applications/utilities/mesh/manipulation/orientFaceZone/orientFaceZone.C
@@ -46,12 +46,13 @@ using namespace Foam;
 
 int main(int argc, char *argv[])
 {
+#   include "addRegionOption.H"
     argList::validArgs.append("faceZone");
     argList::validArgs.append("outsidePoint");
 
 #   include "setRootCase.H"
 #   include "createTime.H"
-#   include "createPolyMesh.H"
+#   include "createNamedPolyMesh.H"
 
     const word zoneName  = args[1];
     const point outsidePoint = args.argRead<point>(2);
diff --git a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
index aa7fa5554bb9ef148ab28ef463ef18156697c0dd..e7e1d52255bb1f9422a8aadbfcc45c856d815c75 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
+++ b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
@@ -104,8 +104,8 @@ int main(int argc, char *argv[])
     );
     argList::addBoolOption
     (
-        "sets",
-        "reconstruct cellSets, faceSets, pointSets"
+        "noSets",
+        "skip reconstructing cellSets, faceSets, pointSets"
     );
     argList::addBoolOption
     (
@@ -122,10 +122,23 @@ int main(int argc, char *argv[])
         args.optionLookup("fields")() >> selectedFields;
     }
 
-    const bool reconstructSets = args.optionFound("sets");
+    const bool noLagrangian = args.optionFound("noLagrangian");
 
+    if (noLagrangian)
+    {
+        Info<< "Skipping reconstructing lagrangian positions and fields"
+            << nl << endl;
+    }
+
+
+    const bool noReconstructSets = args.optionFound("noSets");
+
+    if (noReconstructSets)
+    {
+        Info<< "Skipping reconstructing cellSets, faceSets and pointSets"
+            << nl << endl;
+    }
 
-    const bool noLagrangian = args.optionFound("noLagrangian");
 
     HashSet<word> selectedLagrangianFields;
     if (args.optionFound("lagrangianFields"))
@@ -682,7 +695,7 @@ int main(int argc, char *argv[])
             }
 
 
-            if (reconstructSets)
+            if (!noReconstructSets)
             {
                 // Scan to find all sets
                 HashTable<label> cSetNames;
diff --git a/applications/utilities/preProcessing/mapFields/mapFieldsDict b/applications/utilities/preProcessing/mapFields/mapFieldsDict
index 96fb5a7f8660e27c2d9f788726a6128c91906529..65e33997d938c53ab2f1b2ff5e35a2961b0ec3ef 100644
--- a/applications/utilities/preProcessing/mapFields/mapFieldsDict
+++ b/applications/utilities/preProcessing/mapFields/mapFieldsDict
@@ -14,6 +14,12 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Specify how to map patches. There are three different options:
+// - patch exists in the source case: specify mapping (patchMap)
+// - patch should be interpolated from internal values in source case
+//   (cuttingPatches)
+// - patch should not be mapped. Default if not in patchMap or cuttingPatches
+
 // List of pairs of target/source patches for mapping
 patchMap
 (
diff --git a/src/ODE/ODESolvers/SIBS/SIMPR.C b/src/ODE/ODESolvers/SIBS/SIMPR.C
index d21a931e1d9decfca0c04ddf97efdf70c28a64bc..58a273b03de49f186e24c4cac39045f4e2b0af52 100644
--- a/src/ODE/ODESolvers/SIBS/SIMPR.C
+++ b/src/ODE/ODESolvers/SIBS/SIMPR.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/basicSymmetry/basicSymmetryPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/basic/basicSymmetry/basicSymmetryPointPatchField.C
index 925719ebcc318ef5e4b1c0edc2ed491ddd613f74..85c718e7030313259849b17962c8b605ee3ee5e5 100644
--- a/src/OpenFOAM/fields/pointPatchFields/basic/basicSymmetry/basicSymmetryPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/basic/basicSymmetry/basicSymmetryPointPatchField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/calculated/calculatedPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/basic/calculated/calculatedPointPatchField.C
index 05b0c876d26361fc71d338d860f87e7ae5874e4c..45f98495dbc6f46cba6f60d46658ab058cb20c2f 100644
--- a/src/OpenFOAM/fields/pointPatchFields/basic/calculated/calculatedPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/basic/calculated/calculatedPointPatchField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/coupled/coupledPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/basic/coupled/coupledPointPatchField.C
index 3ad9e51a79cb9f228bac1e366e0cefc077535875..434b6873f4b77ac71e5389a45587fd0e5abe3290 100644
--- a/src/OpenFOAM/fields/pointPatchFields/basic/coupled/coupledPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/basic/coupled/coupledPointPatchField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/value/valuePointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/basic/value/valuePointPatchField.C
index d8fc9c43a38dd274cc63b909385ebac11dff39ff..6f5c7ef2673d2b87af6a4b7acdfd3756f5599e25 100644
--- a/src/OpenFOAM/fields/pointPatchFields/basic/value/valuePointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/basic/value/valuePointPatchField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C
index d6ce1e4cf1c53632caad45e7afcb2c511fed0991..754ba28d7544661b71e5646b1244c19df5a482e9 100644
--- a/src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/OpenFOAM/fields/pointPatchFields/constraint/empty/emptyPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/constraint/empty/emptyPointPatchField.C
index 8838f0fff5b90fab5408eeda2738ded66e010b21..41f97055aea93a93b87cfa07810d81d5c41dea38 100644
--- a/src/OpenFOAM/fields/pointPatchFields/constraint/empty/emptyPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/constraint/empty/emptyPointPatchField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/OpenFOAM/fields/pointPatchFields/constraint/wedge/wedgePointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/constraint/wedge/wedgePointPatchField.C
index 6e8fa53184f183c3a7bdca0a73bde90a6698c301..447c7771cd3afccbde51721af3f197db3354d611 100644
--- a/src/OpenFOAM/fields/pointPatchFields/constraint/wedge/wedgePointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/constraint/wedge/wedgePointPatchField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C
index 63eb0395652752665ac429cd18aa35bac90a268d..be4b10dfd8f3e36fa3a97885d4fb61a767c838e7 100644
--- a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/OpenFOAM/meshes/Identifiers/patch/patchIdentifier.C b/src/OpenFOAM/meshes/Identifiers/patch/patchIdentifier.C
index 6137986165792382e8a5871dca8258a7d6adda04..b69ffae4f4c7b4695a5ae797806ccd43d80205a9 100644
--- a/src/OpenFOAM/meshes/Identifiers/patch/patchIdentifier.C
+++ b/src/OpenFOAM/meshes/Identifiers/patch/patchIdentifier.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,7 @@ License
 
 #include "patchIdentifier.H"
 #include "dictionary.H"
+#include "ListOps.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -79,6 +80,12 @@ Foam::patchIdentifier::~patchIdentifier()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::patchIdentifier::inGroup(const word& name) const
+{
+    return findIndex(inGroups_, name) != -1;
+}
+
+
 void Foam::patchIdentifier::write(Ostream& os) const
 {
     if (physicalType_.size())
diff --git a/src/OpenFOAM/meshes/Identifiers/patch/patchIdentifier.H b/src/OpenFOAM/meshes/Identifiers/patch/patchIdentifier.H
index cb7d60b7b5b153ba230ccb1b3fe98645ec5d34bd..e54d27673302e5b298e30e7ac8429eded6f92eb6 100644
--- a/src/OpenFOAM/meshes/Identifiers/patch/patchIdentifier.H
+++ b/src/OpenFOAM/meshes/Identifiers/patch/patchIdentifier.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -154,6 +154,8 @@ public:
             return inGroups_;
         }
 
+        //- Test if in group
+        bool inGroup(const word&) const;
 
         //- Write patchIdentifier as a dictionary
         void write(Ostream&) const;
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
index a1935d03684133fe0b5e4d645224389f2f79d430..1e104b4602e43033b178fe6a57c12ef4d8e95c5c 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
@@ -486,6 +486,57 @@ Foam::polyBoundaryMesh::groupPatchIDs() const
 }
 
 
+void Foam::polyBoundaryMesh::setGroup
+(
+    const word& groupName,
+    const labelList& patchIDs
+)
+{
+    groupPatchIDsPtr_.clear();
+
+    polyPatchList& patches = *this;
+
+    boolList donePatch(patches.size(), false);
+
+    // Add to specified patches
+    forAll(patchIDs, i)
+    {
+        label patchI = patchIDs[i];
+        polyPatch& pp = patches[patchI];
+
+        if (!pp.inGroup(groupName))
+        {
+            pp.inGroups().append(groupName);
+        }
+        donePatch[patchI] = true;
+    }
+
+    // Remove from other patches
+    forAll(patches, patchI)
+    {
+        if (!donePatch[patchI])
+        {
+            polyPatch& pp = patches[patchI];
+
+            label newI = 0;
+            if (pp.inGroup(groupName))
+            {
+                wordList& groups = pp.inGroups();
+
+                forAll(groups, i)
+                {
+                    if (groups[i] != groupName)
+                    {
+                        groups[newI++] = groups[i];
+                    }
+                }
+                groups.setSize(newI);
+            }
+        }
+    }
+}
+
+
 Foam::wordList Foam::polyBoundaryMesh::names() const
 {
     const polyPatchList& patches = *this;
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
index 3878c80b36629d47ca53419442460bdd6ad80c55..1906b422ac5587b14286da9693cfb868faad2f1c 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
@@ -186,6 +186,9 @@ public:
         //- Per patch group the patch indices
         const HashTable<labelList, word>& groupPatchIDs() const;
 
+        //- Set/add group with patches
+        void setGroup(const word& groupName, const labelList& patchIDs);
+
         //- Return the set of patch IDs corresponding to the given names
         //  By default warns if given names are not found. Optionally
         //  matches to patchGroups as well as patchNames
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index 7bb72786b8331b5210f47b3a8c11e2beef482cc4..316d376defd8aa15a79fcc8b931c9e16393df682 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -2472,7 +2472,7 @@ void Foam::autoLayerDriver::getLayerCellsFaces
 }
 
 
-bool Foam::autoLayerDriver::writeLayerData
+void Foam::autoLayerDriver::printLayerData
 (
     const fvMesh& mesh,
     const labelList& patchIDs,
@@ -2481,153 +2481,284 @@ bool Foam::autoLayerDriver::writeLayerData
     const scalarField& faceRealThickness
 ) const
 {
-    bool allOk = true;
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 
+    // Find maximum length of a patch name, for a nicer output
+    label maxPatchNameLen = 0;
+    forAll(patchIDs, i)
     {
-        label nAdded = 0;
-        forAll(cellNLayers, cellI)
+        label patchI = patchIDs[i];
+        word patchName = pbm[patchI].name();
+        maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
+    }
+
+    Info<< nl
+        << setf(ios_base::left) << setw(maxPatchNameLen) << "patch"
+        << setw(0) << " faces    layers   overall thickness" << nl
+        << setf(ios_base::left) << setw(maxPatchNameLen) << " "
+        << setw(0) << "                   [m]       [%]" << nl
+        << setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
+        << setw(0) << " -----    ------   ---       ---" << endl;
+
+
+    forAll(patchIDs, i)
+    {
+        label patchI = patchIDs[i];
+        const polyPatch& pp = pbm[patchI];
+
+        label sumSize = pp.size();
+
+        // Number of layers
+        const labelList& faceCells = pp.faceCells();
+        label sumNLayers = 0;
+        forAll(faceCells, i)
         {
-            if (cellNLayers[cellI] > 0)
-            {
-                nAdded++;
-            }
+            sumNLayers += cellNLayers[faceCells[i]];
         }
-        cellSet addedCellSet(mesh, "addedCells", nAdded);
-        forAll(cellNLayers, cellI)
+
+        // Thickness
+        scalarField::subField patchWanted = pbm[patchI].patchSlice
+        (
+            faceWantedThickness
+        );
+        scalarField::subField patchReal = pbm[patchI].patchSlice
+        (
+            faceRealThickness
+        );
+
+        scalar sumRealThickness = sum(patchReal);
+        scalar sumFraction = 0;
+        forAll(patchReal, i)
         {
-            if (cellNLayers[cellI] > 0)
+            if (patchWanted[i] > VSMALL)
             {
-                addedCellSet.insert(cellI);
+                sumFraction += (patchReal[i]/patchWanted[i]);
             }
         }
-        addedCellSet.instance() = meshRefiner_.timeName();
-        Info<< "Writing "
-            << returnReduce(addedCellSet.size(), sumOp<label>())
-            << " added cells to cellSet "
-            << addedCellSet.name() << endl;
-        bool ok = addedCellSet.write();
-        allOk = allOk & ok;
+
+
+        reduce(sumSize, sumOp<label>());
+        reduce(sumNLayers, sumOp<label>());
+        reduce(sumRealThickness, sumOp<scalar>());
+        reduce(sumFraction, sumOp<scalar>());
+
+
+        scalar avgLayers = 0;
+        scalar avgReal = 0;
+        scalar avgFraction = 0;
+        if (sumSize > 0)
+        {
+            avgLayers = scalar(sumNLayers)/sumSize;
+            avgReal = sumRealThickness/sumSize;
+            avgFraction = sumFraction/sumSize;
+        }
+
+        Info<< setf(ios_base::left) << setw(maxPatchNameLen)
+            << pbm[patchI].name() << setprecision(3)
+            << " " << setw(8) << sumSize
+            << " " << setw(8) << avgLayers
+            << " " << setw(8) << avgReal
+            << "  " << setw(8) << 100*avgFraction
+            << endl;
     }
+    Info<< endl;
+}
+
+
+bool Foam::autoLayerDriver::writeLayerData
+(
+    const fvMesh& mesh,
+    const labelList& patchIDs,
+    const labelList& cellNLayers,
+    const scalarField& faceWantedThickness,
+    const scalarField& faceRealThickness
+) const
+{
+    bool allOk = true;
+
+    if (meshRefinement::writeLevel() & meshRefinement::WRITELAYERSETS)
     {
-        label nAdded = 0;
-        forAll(faceRealThickness, faceI)
         {
-            if (faceRealThickness[faceI] > 0)
+            label nAdded = 0;
+            forAll(cellNLayers, cellI)
             {
-                nAdded++;
+                if (cellNLayers[cellI] > 0)
+                {
+                    nAdded++;
+                }
             }
+            cellSet addedCellSet(mesh, "addedCells", nAdded);
+            forAll(cellNLayers, cellI)
+            {
+                if (cellNLayers[cellI] > 0)
+                {
+                    addedCellSet.insert(cellI);
+                }
+            }
+            addedCellSet.instance() = meshRefiner_.timeName();
+            Info<< "Writing "
+                << returnReduce(addedCellSet.size(), sumOp<label>())
+                << " added cells to cellSet "
+                << addedCellSet.name() << endl;
+            bool ok = addedCellSet.write();
+            allOk = allOk & ok;
         }
-
-        faceSet layerFacesSet(mesh, "layerFaces", nAdded);
-        forAll(faceRealThickness, faceI)
         {
-            if (faceRealThickness[faceI] > 0)
+            label nAdded = 0;
+            forAll(faceRealThickness, faceI)
             {
-                layerFacesSet.insert(faceI);
+                if (faceRealThickness[faceI] > 0)
+                {
+                    nAdded++;
+                }
             }
+
+            faceSet layerFacesSet(mesh, "layerFaces", nAdded);
+            forAll(faceRealThickness, faceI)
+            {
+                if (faceRealThickness[faceI] > 0)
+                {
+                    layerFacesSet.insert(faceI);
+                }
+            }
+            layerFacesSet.instance() = meshRefiner_.timeName();
+            Info<< "Writing "
+                << returnReduce(layerFacesSet.size(), sumOp<label>())
+                << " faces inside added layer to faceSet "
+                << layerFacesSet.name() << endl;
+            bool ok = layerFacesSet.write();
+            allOk = allOk & ok;
         }
-        layerFacesSet.instance() = meshRefiner_.timeName();
-        Info<< "Writing "
-            << returnReduce(layerFacesSet.size(), sumOp<label>())
-            << " faces inside added layer to faceSet "
-            << layerFacesSet.name() << endl;
-        bool ok = layerFacesSet.write();
-        allOk = allOk & ok;
     }
+
+    if (meshRefinement::writeLevel() & meshRefinement::WRITELAYERFIELDS)
     {
-        volScalarField fld
-        (
-            IOobject
+        {
+            volScalarField fld
             (
-                "nSurfaceLayers",
-                mesh.time().timeName(),
+                IOobject
+                (
+                    "nSurfaceLayers",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE,
+                    false
+                ),
                 mesh,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE,
-                false
-            ),
-            mesh,
-            dimensionedScalar("zero", dimless, 0),
-            fixedValueFvPatchScalarField::typeName
-        );
-        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
-        forAll(patchIDs, i)
-        {
-            label patchI = patchIDs[i];
-            const polyPatch& pp = pbm[patchI];
-            const labelList& faceCells = pp.faceCells();
-            scalarField pfld(faceCells.size());
-            forAll(faceCells, i)
+                dimensionedScalar("zero", dimless, 0),
+                fixedValueFvPatchScalarField::typeName
+            );
+            const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+            forAll(patchIDs, i)
             {
-                pfld[i] = cellNLayers[faceCells[i]];
+                label patchI = patchIDs[i];
+                const polyPatch& pp = pbm[patchI];
+                const labelList& faceCells = pp.faceCells();
+                scalarField pfld(faceCells.size());
+                forAll(faceCells, i)
+                {
+                    pfld[i] = cellNLayers[faceCells[i]];
+                }
+                fld.boundaryField()[patchI] == pfld;
             }
-            fld.boundaryField()[patchI] == pfld;
+            Info<< "Writing volScalarField " << fld.name()
+                << " with actual number of layers" << endl;
+            bool ok = fld.write();
+            allOk = allOk & ok;
         }
-        Info<< "Writing volScalarField " << fld.name()
-            << " with actual number of layers" << endl;
-        bool ok = fld.write();
-        allOk = allOk & ok;
-    }
-    {
-        volScalarField fld
-        (
-            IOobject
+        {
+            volScalarField fld
             (
-                "wantedThickness",
-                mesh.time().timeName(),
+                IOobject
+                (
+                    "thickness",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE,
+                    false
+                ),
                 mesh,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE,
-                false
-            ),
-            mesh,
-            dimensionedScalar("zero", dimless, 0),
-            fixedValueFvPatchScalarField::typeName
-        );
-        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
-        forAll(patchIDs, i)
+                dimensionedScalar("zero", dimless, 0),
+                fixedValueFvPatchScalarField::typeName
+            );
+            const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+            forAll(patchIDs, i)
+            {
+                label patchI = patchIDs[i];
+                fld.boundaryField()[patchI] == pbm[patchI].patchSlice
+                (
+                    faceRealThickness
+                );
+            }
+            Info<< "Writing volScalarField " << fld.name()
+                << " with overall layer thickness" << endl;
+            bool ok = fld.write();
+            allOk = allOk & ok;
+        }
         {
-            label patchI = patchIDs[i];
-            fld.boundaryField()[patchI] == pbm[patchI].patchSlice
+            volScalarField fld
             (
-                faceWantedThickness
+                IOobject
+                (
+                    "thicknessFraction",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE,
+                    false
+                ),
+                mesh,
+                dimensionedScalar("zero", dimless, 0),
+                fixedValueFvPatchScalarField::typeName
             );
+            const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+            forAll(patchIDs, i)
+            {
+                label patchI = patchIDs[i];
+
+                scalarField::subField patchWanted = pbm[patchI].patchSlice
+                (
+                    faceWantedThickness
+                );
+                scalarField::subField patchReal = pbm[patchI].patchSlice
+                (
+                    faceRealThickness
+                );
+
+                // Convert patchReal to relavtive thickness
+                scalarField pfld(patchReal.size(), 0.0);
+                forAll(patchReal, i)
+                {
+                    if (patchWanted[i] > VSMALL)
+                    {
+                        pfld[i] = patchReal[i]/patchWanted[i];
+                    }
+                }
+
+                fld.boundaryField()[patchI] == pfld;
+            }
+            Info<< "Writing volScalarField " << fld.name()
+                << " with overall layer thickness as fraction"
+                << " of desired thickness" << endl;
+            bool ok = fld.write();
+            allOk = allOk & ok;
         }
-        Info<< "Writing volScalarField " << fld.name()
-            << " with wanted thickness" << endl;
-        bool ok = fld.write();
-        allOk = allOk & ok;
     }
+
+    //if (meshRefinement::outputLevel() & meshRefinement::OUTPUTLAYERINFO)
     {
-        volScalarField fld
+        printLayerData
         (
-            IOobject
-            (
-                "thickness",
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE,
-                false
-            ),
             mesh,
-            dimensionedScalar("zero", dimless, 0),
-            fixedValueFvPatchScalarField::typeName
+            patchIDs,
+            cellNLayers,
+            faceWantedThickness,
+            faceRealThickness
         );
-        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
-        forAll(patchIDs, i)
-        {
-            label patchI = patchIDs[i];
-            fld.boundaryField()[patchI] == pbm[patchI].patchSlice
-            (
-                faceRealThickness
-            );
-        }
-        Info<< "Writing volScalarField " << fld.name()
-            << " with layer thickness" << endl;
-        bool ok = fld.write();
-        allOk = allOk & ok;
     }
+
     return allOk;
 }
 
@@ -2737,7 +2868,12 @@ void Foam::autoLayerDriver::addLayers
             << meshRefiner_.timeName() << endl;
         meshRefiner_.write
         (
-            debug,
+            meshRefinement::debugType(debug),
+            meshRefinement::writeType
+            (
+                meshRefinement::writeLevel()
+              | meshRefinement::WRITEMESH
+            ),
             mesh.time().path()/meshRefiner_.timeName()
         );
     }
@@ -3173,7 +3309,12 @@ void Foam::autoLayerDriver::addLayers
 
             meshRefiner_.write
             (
-                debug,
+                meshRefinement::debugType(debug),
+                meshRefinement::writeType
+                (
+                    meshRefinement::writeLevel()
+                  | meshRefinement::WRITEMESH
+                ),
                 mesh.time().path()/meshRefiner_.timeName()
             );
         }
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
index 4534ebe00d161c586218215e9d0038becf2bcef9..58464ad495419029a2f756b4f7b7b55f0f0adcfb 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
@@ -380,6 +380,16 @@ class autoLayerDriver
                     scalarField& faceRealThickness
                 );
 
+                //- Print layer coverage table
+                void printLayerData
+                (
+                    const fvMesh& mesh,
+                    const labelList& patchIDs,
+                    const labelList& cellNLayers,
+                    const scalarField& faceWantedThickness,
+                    const scalarField& faceRealThickness
+                ) const;
+
                 //- Write cellSet,faceSet for layers
                 bool writeLayerData
                 (
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
index 3e824b1e0287624034d6af5e6e993783a6164e2b..328b4fd86a13dd472583808a434c525f4d98e9fe 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
@@ -1330,7 +1330,12 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
         meshRefiner_.mesh().setInstance(meshRefiner_.timeName());
         meshRefiner_.write
         (
-            debug,
+            meshRefinement::debugType(debug),
+            meshRefinement::writeType
+            (
+                meshRefinement::writeLevel()
+              | meshRefinement::WRITEMESH
+            ),
             mesh.time().path()/meshRefiner_.timeName()
         );
         dispVec.write();
@@ -1775,7 +1780,12 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
 
         meshRefiner_.write
         (
-            debug,
+            meshRefinement::debugType(debug),
+            meshRefinement::writeType
+            (
+                meshRefinement::writeLevel()
+              | meshRefinement::WRITEMESH
+            ),
             mesh.time().path()/meshRefiner_.timeName()
         );
         dispVec.write();
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index a53b6f11eb3bf7b9adf66ec8f04d1c93e7272588..ce97c27104e7b438c650fb4e78203a7538c746bd 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -698,7 +698,16 @@ void Foam::autoRefineDriver::removeInsideCells
     {
         Pout<< "Writing subsetted mesh to time "
             << meshRefiner_.timeName() << '.' << endl;
-        meshRefiner_.write(debug, mesh.time().path()/meshRefiner_.timeName());
+        meshRefiner_.write
+        (
+            meshRefinement::debugType(debug),
+            meshRefinement::writeType
+            (
+                meshRefinement::writeLevel()
+              | meshRefinement::WRITEMESH
+            ),
+            mesh.time().path()/meshRefiner_.timeName()
+        );
         Pout<< "Dumped mesh in = "
             << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
     }
@@ -956,7 +965,12 @@ void Foam::autoRefineDriver::zonify
                 << meshRefiner_.timeName() << '.' << endl;
             meshRefiner_.write
             (
-                debug,
+                meshRefinement::debugType(debug),
+                meshRefinement::writeType
+                (
+                    meshRefinement::writeLevel()
+                  | meshRefinement::WRITEMESH
+                ),
                 mesh.time().path()/meshRefiner_.timeName()
             );
         }
@@ -1071,7 +1085,16 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
     {
         Pout<< "Writing handleProblemCells mesh to time "
             << meshRefiner_.timeName() << '.' << endl;
-        meshRefiner_.write(debug, mesh.time().path()/meshRefiner_.timeName());
+        meshRefiner_.write
+        (
+            meshRefinement::debugType(debug),
+            meshRefinement::writeType
+            (
+                meshRefinement::writeLevel()
+              | meshRefinement::WRITEMESH
+            ),
+            mesh.time().path()/meshRefiner_.timeName()
+        );
     }
 }
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index adb7c3859032a928765510d046141f802df54794..03017d8f8c0be323df2ffb2b90d6fa68d43955e4 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -730,7 +730,12 @@ void Foam::autoSnapDriver::preSmoothPatch
             << meshRefiner.timeName() << '.' << endl;
         meshRefiner.write
         (
-            debug,
+            meshRefinement::debugType(debug),
+            meshRefinement::writeType
+            (
+                meshRefinement::writeLevel()
+              | meshRefinement::WRITEMESH
+            ),
             mesh.time().path()/meshRefiner.timeName()
         );
         Info<< "Dumped mesh in = "
@@ -2027,7 +2032,12 @@ void Foam::autoSnapDriver::smoothDisplacement
 
         meshRefiner_.write
         (
-            debug,
+            meshRefinement::debugType(debug),
+            meshRefinement::writeType
+            (
+                meshRefinement::writeLevel()
+              | meshRefinement::WRITEMESH
+            ),
             mesh.time().path()/meshRefiner_.timeName()
         );
         Info<< "Writing displacement field ..." << endl;
@@ -2497,8 +2507,13 @@ void Foam::autoSnapDriver::doSnap
                     << endl;
                 meshRefiner_.write
                 (
-                    debug, mesh.time().path()
-                   /"duplicatedPoints"
+                    meshRefinement::debugType(debug),
+                    meshRefinement::writeType
+                    (
+                        meshRefinement::writeLevel()
+                      | meshRefinement::WRITEMESH
+                    ),
+                    mesh.time().path()/"duplicatedPoints"
                 );
             }
         }
@@ -2804,7 +2819,12 @@ void Foam::autoSnapDriver::doSnap
                     << meshRefiner_.timeName() << endl;
                 meshRefiner_.write
                 (
-                    debug,
+                    meshRefinement::debugType(debug),
+                    meshRefinement::writeType
+                    (
+                        meshRefinement::writeLevel()
+                      | meshRefinement::WRITEMESH
+                    ),
                     mesh.time().path()/meshRefiner_.timeName()
                 );
                 Info<< "Writing displacement field ..." << endl;
@@ -2866,7 +2886,12 @@ void Foam::autoSnapDriver::doSnap
             << meshRefiner_.timeName() << endl;
         meshRefiner_.write
         (
-            debug,
+            meshRefinement::debugType(debug),
+            meshRefinement::writeType
+            (
+                meshRefinement::writeLevel()
+              | meshRefinement::WRITEMESH
+            ),
             meshRefiner_.timeName()
         );
     }
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index b8282bedfcd8e5f0399be522e5eccbfe2b1bf827..496f980a118c9a52ead5cc4bed659865532638f1 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -131,13 +131,13 @@ class autoSnapDriver
                     const List<pointConstraint>& constraints,
                     vectorField& disp
                 ) const;
-                void smoothAndConstrain2
-                (
-                    const bool applyConstraints,
-                    const indirectPrimitivePatch& pp,
-                    const List<pointConstraint>& constraints,
-                    vectorField& disp
-                ) const;
+                //void smoothAndConstrain2
+                //(
+                //    const bool applyConstraints,
+                //    const indirectPrimitivePatch& pp,
+                //    const List<pointConstraint>& constraints,
+                //    vectorField& disp
+                //) const;
                 void calcNearest
                 (
                     const label iter,
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index 32133ad5f8437c3e754777792bf6385a1260c951..92ef64bf4577b030da9550bf30f6a65bf8351c25 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -25,10 +25,9 @@ License
 
 #include "autoSnapDriver.H"
 #include "polyTopoChange.H"
-#include "OFstream.H"
 #include "syncTools.H"
 #include "fvMesh.H"
-#include "OFstream.H"
+#include "OBJstream.H"
 #include "motionSmoother.H"
 #include "refinementSurfaces.H"
 #include "refinementFeatures.H"
@@ -216,77 +215,77 @@ void Foam::autoSnapDriver::smoothAndConstrain
 }
 //XXXXXX
 //TODO: make proper parallel so coupled edges don't have double influence
-void Foam::autoSnapDriver::smoothAndConstrain2
-(
-    const bool applyConstraints,
-    const indirectPrimitivePatch& pp,
-    const List<pointConstraint>& constraints,
-    vectorField& disp
-) const
-{
-    const fvMesh& mesh = meshRefiner_.mesh();
-
-    for (label avgIter = 0; avgIter < 20; avgIter++)
-    {
-        vectorField dispSum(pp.nPoints(), vector::zero);
-        labelList dispCount(pp.nPoints(), 0);
-
-        const labelListList& pointEdges = pp.pointEdges();
-        const edgeList& edges = pp.edges();
-
-        forAll(pointEdges, pointI)
-        {
-            const labelList& pEdges = pointEdges[pointI];
-
-            forAll(pEdges, i)
-            {
-                label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
-                dispSum[pointI] += disp[nbrPointI];
-                dispCount[pointI]++;
-            }
-        }
-
-        syncTools::syncPointList
-        (
-            mesh,
-            pp.meshPoints(),
-            dispSum,
-            plusEqOp<point>(),
-            vector::zero,
-            mapDistribute::transform()
-        );
-        syncTools::syncPointList
-        (
-            mesh,
-            pp.meshPoints(),
-            dispCount,
-            plusEqOp<label>(),
-            0,
-            mapDistribute::transform()
-        );
-
-        // Constraints
-        forAll(constraints, pointI)
-        {
-            if (dispCount[pointI] > 0)// && constraints[pointI].first() <= 1)
-            {
-                // Mix my displacement with neighbours' displacement
-                disp[pointI] =
-                    0.5
-                   *(disp[pointI] + dispSum[pointI]/dispCount[pointI]);
-
-                if (applyConstraints)
-                {
-                    disp[pointI] = transform
-                    (
-                        constraints[pointI].constraintTransformation(),
-                        disp[pointI]
-                    );
-                }
-            }
-        }
-    }
-}
+//void Foam::autoSnapDriver::smoothAndConstrain2
+//(
+//    const bool applyConstraints,
+//    const indirectPrimitivePatch& pp,
+//    const List<pointConstraint>& constraints,
+//    vectorField& disp
+//) const
+//{
+//    const fvMesh& mesh = meshRefiner_.mesh();
+//
+//    for (label avgIter = 0; avgIter < 20; avgIter++)
+//    {
+//        vectorField dispSum(pp.nPoints(), vector::zero);
+//        labelList dispCount(pp.nPoints(), 0);
+//
+//        const labelListList& pointEdges = pp.pointEdges();
+//        const edgeList& edges = pp.edges();
+//
+//        forAll(pointEdges, pointI)
+//        {
+//            const labelList& pEdges = pointEdges[pointI];
+//
+//            forAll(pEdges, i)
+//            {
+//                label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
+//                dispSum[pointI] += disp[nbrPointI];
+//                dispCount[pointI]++;
+//            }
+//        }
+//
+//        syncTools::syncPointList
+//        (
+//            mesh,
+//            pp.meshPoints(),
+//            dispSum,
+//            plusEqOp<point>(),
+//            vector::zero,
+//            mapDistribute::transform()
+//        );
+//        syncTools::syncPointList
+//        (
+//            mesh,
+//            pp.meshPoints(),
+//            dispCount,
+//            plusEqOp<label>(),
+//            0,
+//            mapDistribute::transform()
+//        );
+//
+//        // Constraints
+//        forAll(constraints, pointI)
+//        {
+//            if (dispCount[pointI] > 0)// && constraints[pointI].first() <= 1)
+//            {
+//                // Mix my displacement with neighbours' displacement
+//                disp[pointI] =
+//                    0.5
+//                   *(disp[pointI] + dispSum[pointI]/dispCount[pointI]);
+//
+//                if (applyConstraints)
+//                {
+//                    disp[pointI] = transform
+//                    (
+//                        constraints[pointI].constraintTransformation(),
+//                        disp[pointI]
+//                    );
+//                }
+//            }
+//        }
+//    }
+//}
 //XXXXXX
 
 
@@ -1082,16 +1081,13 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
     List<pointConstraint>& patchConstraints
 ) const
 {
-    autoPtr<OFstream> feStr;
-    label feVertI = 0;
-    autoPtr<OFstream> fpStr;
-    label fpVertI = 0;
-
+    autoPtr<OBJstream> feStr;
+    autoPtr<OBJstream> fpStr;
     if (debug&meshRefinement::OBJINTERSECTIONS)
     {
         feStr.reset
         (
-            new OFstream
+            new OBJstream
             (
                 meshRefiner_.mesh().time().path()
               / "implicitFeatureEdge_" + name(iter) + ".obj"
@@ -1102,7 +1098,7 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
 
         fpStr.reset
         (
-            new OFstream
+            new OBJstream
             (
                 meshRefiner_.mesh().time().path()
               / "implicitFeaturePoint_" + name(iter) + ".obj"
@@ -1139,7 +1135,10 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
         if
         (
             (constraint.first() > patchConstraints[pointI].first())
-         || (magSqr(attraction) < magSqr(patchAttraction[pointI]))
+         || (
+                (constraint.first() == patchConstraints[pointI].first())
+             && (magSqr(attraction) < magSqr(patchAttraction[pointI]))
+            )
         )
         {
             patchAttraction[pointI] = attraction;
@@ -1149,19 +1148,11 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
 
             if (patchConstraints[pointI].first() == 2 && feStr.valid())
             {
-                meshTools::writeOBJ(feStr(), pt);
-                feVertI++;
-                meshTools::writeOBJ(feStr(), pt+patchAttraction[pointI]);
-                feVertI++;
-                feStr() << "l " << feVertI-1 << ' ' << feVertI << nl;
+                feStr().write(linePointRef(pt, pt+patchAttraction[pointI]));
             }
             else if (patchConstraints[pointI].first() == 3 && fpStr.valid())
             {
-                meshTools::writeOBJ(fpStr(), pt);
-                fpVertI++;
-                meshTools::writeOBJ(fpStr(), pt+patchAttraction[pointI]);
-                fpVertI++;
-                fpStr() << "l " << fpVertI-1 << ' ' << fpVertI << nl;
+                fpStr().write(linePointRef(pt, pt+patchAttraction[pointI]));
             }
         }
     }
@@ -1604,18 +1595,15 @@ void Foam::autoSnapDriver::determineFeatures
     List<pointConstraint>& patchConstraints
 ) const
 {
-    autoPtr<OFstream> featureEdgeStr;
-    label featureEdgeVertI = 0;
-    autoPtr<OFstream> missedEdgeStr;
-    label missedVertI = 0;
-    autoPtr<OFstream> featurePointStr;
-    label featurePointVertI = 0;
+    autoPtr<OBJstream> featureEdgeStr;
+    autoPtr<OBJstream> missedEdgeStr;
+    autoPtr<OBJstream> featurePointStr;
 
     if (debug&meshRefinement::OBJINTERSECTIONS)
     {
         featureEdgeStr.reset
         (
-            new OFstream
+            new OBJstream
             (
                 meshRefiner_.mesh().time().path()
               / "featureEdge_" + name(iter) + ".obj"
@@ -1626,7 +1614,7 @@ void Foam::autoSnapDriver::determineFeatures
 
         missedEdgeStr.reset
         (
-            new OFstream
+            new OBJstream
             (
                 meshRefiner_.mesh().time().path()
               / "missedFeatureEdge_" + name(iter) + ".obj"
@@ -1637,7 +1625,7 @@ void Foam::autoSnapDriver::determineFeatures
 
         featurePointStr.reset
         (
-            new OFstream
+            new OBJstream
             (
                 meshRefiner_.mesh().time().path()
               / "featurePoint_" + name(iter) + ".obj"
@@ -1677,7 +1665,10 @@ void Foam::autoSnapDriver::determineFeatures
         if
         (
             (constraint.first() > patchConstraints[pointI].first())
-         || (magSqr(attraction) < magSqr(patchAttraction[pointI]))
+         || (
+                (constraint.first() == patchConstraints[pointI].first())
+             && (magSqr(attraction) < magSqr(patchAttraction[pointI]))
+            )
         )
         {
             patchAttraction[pointI] = attraction;
@@ -1724,34 +1715,20 @@ void Foam::autoSnapDriver::determineFeatures
                             // Dump
                             if (featureEdgeStr.valid())
                             {
-                                meshTools::writeOBJ(featureEdgeStr(), pt);
-                                featureEdgeVertI++;
-                                meshTools::writeOBJ
+                                featureEdgeStr().write
                                 (
-                                    featureEdgeStr(),
-                                    nearInfo.hitPoint()
+                                    linePointRef(pt, nearInfo.hitPoint())
                                 );
-                                featureEdgeVertI++;
-                                featureEdgeStr()
-                                    << "l " << featureEdgeVertI-1 << ' '
-                                    << featureEdgeVertI << nl;
                             }
                         }
                         else
                         {
                             if (missedEdgeStr.valid())
                             {
-                                meshTools::writeOBJ(missedEdgeStr(), pt);
-                                missedVertI++;
-                                meshTools::writeOBJ
+                                missedEdgeStr().write
                                 (
-                                    missedEdgeStr(),
-                                    nearInfo.missPoint()
+                                    linePointRef(pt, nearInfo.missPoint())
                                 );
-                                missedVertI++;
-                                missedEdgeStr()
-                                    << "l " << missedVertI-1 << ' '
-                                    << missedVertI << nl;
                             }
                         }
                     }
@@ -1788,34 +1765,20 @@ void Foam::autoSnapDriver::determineFeatures
                     // Dump
                     if (featureEdgeStr.valid())
                     {
-                        meshTools::writeOBJ(featureEdgeStr(), pt);
-                        featureEdgeVertI++;
-                        meshTools::writeOBJ
+                        featureEdgeStr().write
                         (
-                            featureEdgeStr(),
-                            nearInfo.hitPoint()
+                            linePointRef(pt, nearInfo.hitPoint())
                         );
-                        featureEdgeVertI++;
-                        featureEdgeStr()
-                            << "l " << featureEdgeVertI-1 << ' '
-                            << featureEdgeVertI << nl;
                     }
                 }
                 else
                 {
                     if (missedEdgeStr.valid())
                     {
-                        meshTools::writeOBJ(missedEdgeStr(), pt);
-                        missedVertI++;
-                        meshTools::writeOBJ
+                        missedEdgeStr().write
                         (
-                            missedEdgeStr(),
-                            nearInfo.missPoint()
+                            linePointRef(pt, nearInfo.missPoint())
                         );
-                        missedVertI++;
-                        missedEdgeStr()
-                            << "l " << missedVertI-1 << ' '
-                            << missedVertI << nl;
                     }
                 }
             }
@@ -1852,13 +1815,7 @@ void Foam::autoSnapDriver::determineFeatures
                         const point& featPt =
                             shapes.points()[nearInfo.second()];
 
-                        meshTools::writeOBJ(featurePointStr(), pt);
-                        featurePointVertI++;
-                        meshTools::writeOBJ(featurePointStr(), featPt);
-                        featurePointVertI++;
-                        featurePointStr()
-                            << "l " << featurePointVertI-1 << ' '
-                            << featurePointVertI << nl;
+                        featurePointStr().write(linePointRef(pt, featPt));
                     }
                 }
             }
@@ -2000,13 +1957,12 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
         const scalar baffleFeatureCos = Foam::cos(degToRad(91));
 
 
-        autoPtr<OFstream> baffleEdgeStr;
-        label baffleEdgeVertI = 0;
+        autoPtr<OBJstream> baffleEdgeStr;
         if (debug&meshRefinement::OBJINTERSECTIONS)
         {
             baffleEdgeStr.reset
             (
-                new OFstream
+                new OBJstream
                 (
                     meshRefiner_.mesh().time().path()
                   / "baffleEdge_" + name(iter) + ".obj"
@@ -2040,12 +1996,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                 {
                     const point& p0 = pp.localPoints()[e[0]];
                     const point& p1 = pp.localPoints()[e[1]];
-                    meshTools::writeOBJ(baffleEdgeStr(), p0);
-                    baffleEdgeVertI++;
-                    meshTools::writeOBJ(baffleEdgeStr(), p1);
-                    baffleEdgeVertI++;
-                    baffleEdgeStr() << "l " << baffleEdgeVertI-1
-                        << ' ' << baffleEdgeVertI << nl;
+                    baffleEdgeStr().write(linePointRef(p0, p1));
                 }
             }
         }
@@ -2226,6 +2177,22 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     {
         const fvMesh& mesh = meshRefiner_.mesh();
 
+        //autoPtr<OBJstream> attractStr;
+        //if (debug&meshRefinement::OBJINTERSECTIONS)
+        //{
+        //    attractStr.reset
+        //    (
+        //        new OBJstream
+        //        (
+        //            meshRefiner_.mesh().time().path()
+        //          / "initAttract_" + name(iter) + ".obj"
+        //        )
+        //    );
+        //    Info<< nl << "Dumping initial attract points to "
+        //        << attractStr().name() << endl;
+        //}
+
+
         boolList isFeatureEdgeOrPoint(pp.nPoints(), false);
         label nFeats = 0;
         forAll(rawPatchConstraints, pointI)
@@ -2234,6 +2201,19 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
             {
                 isFeatureEdgeOrPoint[pointI] = true;
                 nFeats++;
+
+                //if (attractStr.valid())
+                //{
+                //    const point& pt = pp.localPoints()[pointI];
+                //    attractStr().write
+                //    (
+                //        linePointRef
+                //        (
+                //            pt,
+                //            pt+rawPatchAttraction[pointI]
+                //        )
+                //    );
+                //}
             }
         }
 
@@ -2254,6 +2234,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
 
         for (label nGrow = 0; nGrow < 1; nGrow++)
         {
+            boolList newIsFeatureEdgeOrPoint(isFeatureEdgeOrPoint);
+
             forAll(pp.localFaces(), faceI)
             {
                 const face& f = pp.localFaces()[faceI];
@@ -2265,12 +2247,15 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                         // Mark all points on face
                         forAll(f, fp)
                         {
-                            isFeatureEdgeOrPoint[f[fp]] = true;
+                            newIsFeatureEdgeOrPoint[f[fp]] = true;
                         }
                         break;
                     }
                 }
             }
+
+            isFeatureEdgeOrPoint = newIsFeatureEdgeOrPoint;
+
             syncTools::syncPointList
             (
                 mesh,
@@ -2282,12 +2267,40 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
         }
 
 
+        //autoPtr<OBJstream> finalStr;
+        //if (debug&meshRefinement::OBJINTERSECTIONS)
+        //{
+        //    finalStr.reset
+        //    (
+        //        new OBJstream
+        //        (
+        //            meshRefiner_.mesh().time().path()
+        //          / "finalAttract_" + name(iter) + ".obj"
+        //        )
+        //    );
+        //    Info<< nl << "Dumping final attract points to "
+        //        << finalStr().name() << endl;
+        //}
+
         // Collect attractPoints
         forAll(isFeatureEdgeOrPoint, pointI)
         {
             if (isFeatureEdgeOrPoint[pointI])
             {
                 attractPoints.append(pointI);
+
+                //if (finalStr.valid())
+                //{
+                //    const point& pt = pp.localPoints()[pointI];
+                //    finalStr().write
+                //    (
+                //        linePointRef
+                //        (
+                //            pt,
+                //            pt+rawPatchAttraction[pointI]
+                //        )
+                //    );
+                //}
             }
         }
 
@@ -2440,12 +2453,12 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     //     (hopefully)
     if (multiRegionFeatureSnap)
     {
-        autoPtr<OFstream> multiPatchStr;
+        autoPtr<OBJstream> multiPatchStr;
         if (debug&meshRefinement::OBJINTERSECTIONS)
         {
             multiPatchStr.reset
             (
-                new OFstream
+                new OBJstream
                 (
                     meshRefiner_.mesh().time().path()
                   / "multiPatch_" + name(iter) + ".obj"
@@ -2540,11 +2553,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
 
                         if (multiPatchStr.valid())
                         {
-                            meshTools::writeOBJ
-                            (
-                                multiPatchStr(),
-                                pp.localPoints()[pointI]
-                            );
+                            multiPatchStr().write(pp.localPoints()[pointI]);
                         }
                     }
                 }
@@ -2561,53 +2570,34 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     // Dump
     if (debug&meshRefinement::OBJINTERSECTIONS)
     {
-        OFstream featureEdgeStr
+        OBJstream featureEdgeStr
         (
             meshRefiner_.mesh().time().path()
           / "edgeAttractors_" + name(iter) + ".obj"
         );
-        label featureEdgeVertI = 0;
         Pout<< "Dumping feature-edge attraction to "
             << featureEdgeStr.name() << endl;
 
-        OFstream featurePointStr
+        OBJstream featurePointStr
         (
             meshRefiner_.mesh().time().path()
           / "pointAttractors_" + name(iter) + ".obj"
         );
-        label featurePointVertI = 0;
         Pout<< "Dumping feature-point attraction to "
             << featurePointStr.name() << endl;
 
         forAll(patchConstraints, pointI)
         {
             const point& pt = pp.localPoints()[pointI];
+            const vector& attr = patchAttraction[pointI];
 
             if (patchConstraints[pointI].first() == 2)
             {
-                meshTools::writeOBJ(featureEdgeStr, pt);
-                featureEdgeVertI++;
-                meshTools::writeOBJ
-                (
-                    featureEdgeStr,
-                    pt+patchAttraction[pointI]
-                );
-                featureEdgeVertI++;
-                featureEdgeStr << "l " << featureEdgeVertI-1
-                    << ' ' << featureEdgeVertI << nl;
+                featureEdgeStr.write(linePointRef(pt, pt+attr));
             }
             else if (patchConstraints[pointI].first() == 3)
             {
-                meshTools::writeOBJ(featurePointStr, pt);
-                featurePointVertI++;
-                meshTools::writeOBJ
-                (
-                    featurePointStr,
-                    pt+patchAttraction[pointI]
-                );
-                featurePointVertI++;
-                featurePointStr << "l " << featurePointVertI-1
-                    << ' ' << featurePointVertI << nl;
+                featurePointStr.write(linePointRef(pt, pt+attr));
             }
         }
     }
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index 2406b8431021fed46c8b0a747966ef32bf1e0f50..b8609ede3e15e54195d3366756aec614ab419eaf 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -62,8 +62,61 @@ License
 namespace Foam
 {
     defineTypeNameAndDebug(meshRefinement, 0);
+
+    template<>
+    const char* Foam::NamedEnum
+    <
+        Foam::meshRefinement::IOdebugType,
+        4
+    >::names[] =
+    {
+        "mesh",
+        //"scalarLevels",
+        "intersections",
+        "featureSeeds",
+        "layerInfo"
+    };
+
+    template<>
+    const char* Foam::NamedEnum
+    <
+        Foam::meshRefinement::IOoutputType,
+        1
+    >::names[] =
+    {
+        "layerInfo"
+    };
+
+    template<>
+    const char* Foam::NamedEnum
+    <
+        Foam::meshRefinement::IOwriteType,
+        4
+    >::names[] =
+    {
+        "mesh",
+        "scalarLevels",
+        "layerSets",
+        "layerFields"
+    };
+
 }
 
+const Foam::NamedEnum<Foam::meshRefinement::IOdebugType, 4>
+Foam::meshRefinement::IOdebugTypeNames;
+
+const Foam::NamedEnum<Foam::meshRefinement::IOoutputType, 1>
+Foam::meshRefinement::IOoutputTypeNames;
+
+const Foam::NamedEnum<Foam::meshRefinement::IOwriteType, 4>
+Foam::meshRefinement::IOwriteTypeNames;
+
+
+Foam::meshRefinement::writeType Foam::meshRefinement::writeLevel_;
+
+Foam::meshRefinement::outputType Foam::meshRefinement::outputLevel_;
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::meshRefinement::calcNeighbourData
@@ -2728,23 +2781,48 @@ void Foam::meshRefinement::dumpIntersections(const fileName& prefix) const
 
 void Foam::meshRefinement::write
 (
-    const label flag,
+    const debugType debugFlags,
+    const writeType writeFlags,
     const fileName& prefix
 ) const
 {
-    if (flag & MESH)
+    if (writeFlags & WRITEMESH)
     {
         write();
     }
-    if (flag & SCALARLEVELS)
+    if (writeFlags & WRITELEVELS)
     {
         dumpRefinementLevel();
     }
-    if (flag & OBJINTERSECTIONS && prefix.size())
+    if (debugFlags & OBJINTERSECTIONS && prefix.size())
     {
         dumpIntersections(prefix);
     }
 }
 
 
+Foam::meshRefinement::writeType Foam::meshRefinement::writeLevel()
+{
+    return writeLevel_;
+}
+
+
+void Foam::meshRefinement::writeLevel(const writeType flags)
+{
+    writeLevel_ = flags;
+}
+
+
+Foam::meshRefinement::outputType Foam::meshRefinement::outputLevel()
+{
+    return outputLevel_;
+}
+
+
+void Foam::meshRefinement::outputLevel(const outputType flags)
+{
+    outputLevel_ = flags;
+}
+
+
 // ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 873167bc4f0227cfb7302e7f79ef7382c6e3e8cb..168cd487fedd9db93711049f4f90e81c92067bbd 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -85,16 +85,52 @@ public:
 
     // Public data types
 
-        //- Enumeration for debug dumping
-        enum writeFlag
+        //- Enumeration for what to debug
+        enum IOdebugType
         {
-            MESH = 1,
-            SCALARLEVELS = 2,
-            OBJINTERSECTIONS = 4,
-            FEATURESEEDS = 8,
-            LAYERINFO = 16
+            IOMESH,
+            //IOSCALARLEVELS,
+            IOOBJINTERSECTIONS,
+            IOFEATURESEEDS,
+            IOLAYERINFO
+        };
+        static const NamedEnum<IOdebugType, 4> IOdebugTypeNames;
+        enum debugType
+        {
+            MESH = 1<<IOMESH,
+            //SCALARLEVELS = 1<<IOSCALARLEVELS,
+            OBJINTERSECTIONS = 1<<IOOBJINTERSECTIONS,
+            FEATURESEEDS = 1<<IOFEATURESEEDS,
+            LAYERINFO = 1<<IOLAYERINFO
         };
 
+        //- Enumeration for what to output
+        enum IOoutputType
+        {
+            IOOUTPUTLAYERINFO
+        };
+        static const NamedEnum<IOoutputType, 1> IOoutputTypeNames;
+        enum outputType
+        {
+            OUTPUTLAYERINFO = 1<<IOOUTPUTLAYERINFO
+        };
+
+        //- Enumeration for what to write
+        enum IOwriteType
+        {
+            IOWRITEMESH,
+            IOWRITELEVELS,
+            IOWRITELAYERSETS,
+            IOWRITELAYERFIELDS
+        };
+        static const NamedEnum<IOwriteType, 4> IOwriteTypeNames;
+        enum writeType
+        {
+            WRITEMESH = 1<<IOWRITEMESH,
+            WRITELEVELS = 1<<IOWRITELEVELS,
+            WRITELAYERSETS = 1<<IOWRITELAYERSETS,
+            WRITELAYERFIELDS = 1<<IOWRITELAYERFIELDS
+        };
 
         //- Enumeration for how the userdata is to be mapped upon refinement.
         enum mapType
@@ -107,6 +143,15 @@ public:
 
 private:
 
+    // Static data members
+
+        //- Control of writing level
+        static writeType writeLevel_;
+
+        //- Control of output/log level
+        static outputType outputLevel_;
+
+
     // Private data
 
         //- Reference to mesh
@@ -1007,9 +1052,13 @@ public:
             //- Debug: Write intersection information to OBJ format
             void dumpIntersections(const fileName& prefix) const;
 
-            //- Do any one of above IO functions. flag is combination of
-            //  writeFlag values.
-            void write(const label flag, const fileName&) const;
+            //- Do any one of above IO functions
+            void write
+            (
+                const debugType debugFlags,
+                const writeType writeFlags,
+                const fileName&
+            ) const;
 
             //- Helper: calculate average
             template<class T>
@@ -1030,6 +1079,21 @@ public:
                 const UList<T>& values
             );
 
+
+            //- Get/set write level
+            static writeType writeLevel();
+            static void writeLevel(const writeType);
+
+            //- Get/set output level
+            static outputType outputLevel();
+            static void outputLevel(const outputType);
+
+
+            //- Helper: convert wordList into bit pattern using provided
+            //  NamedEnum
+            template<class Enum>
+            static int readFlags(const Enum& namedEnum, const wordList&);
+
 };
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 491b9f49d42b901a4235d1b4b5e23f5cfacb81cd..8f4c504c00d65fa7e0fed3f4aca40f5d479ec0d5 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -685,7 +685,7 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::getDuplicateFaces
         << " pairs of duplicate faces." << nl << endl;
 
 
-    if (debug&meshRefinement::MESH)
+    if (debug&MESH)
     {
         faceSet duplicateFaceSet(mesh_, "duplicateFaces", 2*dupI);
 
@@ -790,12 +790,17 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createZoneBaffles
                     << abort(FatalError);
             }
 
-            if (debug&meshRefinement::MESH)
+            if (debug&MESH)
             {
                 const_cast<Time&>(mesh_.time())++;
                 Pout<< "Writing zone-baffled mesh to time " << timeName()
                     << endl;
-                write(debug, mesh_.time().path()/"baffles");
+                write
+                (
+                    debugType(debug),
+                    writeType(writeLevel() | WRITEMESH),
+                    mesh_.time().path()/"baffles"
+                );
             }
         }
         Info<< "Created " << nZoneFaces << " baffles in = "
@@ -1920,9 +1925,9 @@ void Foam::meshRefinement::handleSnapProblems
     Info<< "Analyzed problem cells in = "
         << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
 
-    if (debug&meshRefinement::MESH)
+    if (debug&MESH)
     {
-        faceSet problemFaces(mesh_, "problemFaces", 100);
+        faceSet problemFaces(mesh_, "problemFaces", mesh_.nFaces()/100);
 
         forAll(facePatch, faceI)
         {
@@ -1957,11 +1962,16 @@ void Foam::meshRefinement::handleSnapProblems
 
     printMeshInfo(debug, "After introducing baffles");
 
-    if (debug&meshRefinement::MESH)
+    if (debug&MESH)
     {
         Pout<< "Writing extra baffled mesh to time "
             << timeName() << endl;
-        write(debug, runTime.path()/"extraBaffles");
+        write
+        (
+            debugType(debug),
+            writeType(writeLevel() | WRITEMESH),
+            runTime.path()/"extraBaffles"
+        );
         Pout<< "Dumped debug data in = "
             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
     }
@@ -2519,11 +2529,16 @@ void Foam::meshRefinement::baffleAndSplitMesh
 
     printMeshInfo(debug, "After introducing baffles");
 
-    if (debug&meshRefinement::MESH)
+    if (debug&MESH)
     {
         Pout<< "Writing baffled mesh to time " << timeName()
             << endl;
-        write(debug, runTime.path()/"baffles");
+        write
+        (
+            debugType(debug),
+            writeType(writeLevel() | WRITEMESH),
+            runTime.path()/"baffles"
+        );
         Pout<< "Dumped debug data in = "
             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
     }
@@ -2575,11 +2590,16 @@ void Foam::meshRefinement::baffleAndSplitMesh
 
     printMeshInfo(debug, "After subsetting");
 
-    if (debug&meshRefinement::MESH)
+    if (debug&MESH)
     {
         Pout<< "Writing subsetted mesh to time " << timeName()
             << endl;
-        write(debug, runTime.path()/timeName());
+        write
+        (
+            debugType(debug),
+            writeType(writeLevel() | WRITEMESH),
+            runTime.path()/timeName()
+        );
         Pout<< "Dumped debug data in = "
             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
     }
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
index a1c1ed0f58ccf9128142c3494afa9be2fdd62b3d..f49ad65a32c335f272c8ffb9d3d2937220e1c40b 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
@@ -593,7 +593,12 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
             mesh_.movePoints(newPoints);
             Pout<< "Writing newPoints mesh to time " << timeName()
                 << endl;
-            write(debug, mesh_.time().path()/"newPoints");
+            write
+            (
+                debugType(debug),
+                writeType(writeLevel() | WRITEMESH),
+                mesh_.time().path()/"newPoints"
+            );
             mesh_.movePoints(oldPoints);
         }
     }
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index d2556c4a9d0d92144bd354e37689d79472679fc4..83ed367e89abdf1e6fbac3b549d97c03d24c2a99 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -2253,9 +2253,9 @@ Foam::meshRefinement::refineAndBalance
             << " mesh to time " << timeName() << endl;
         write
         (
-            debug,
-            mesh_.time().path()
-           /timeName()
+            debugType(debug),
+            writeType(writeLevel() | WRITEMESH),
+            mesh_.time().path()/timeName()
         );
         Pout<< "Dumped debug data in = "
             << mesh_.time().cpuTimeIncrement() << " s" << endl;
@@ -2317,7 +2317,8 @@ Foam::meshRefinement::refineAndBalance
                     << " mesh to time " << timeName() << endl;
                 write
                 (
-                    debug,
+                    debugType(debug),
+                    writeType(writeLevel() | WRITEMESH),
                     mesh_.time().path()/timeName()
                 );
                 Pout<< "Dumped debug data in = "
@@ -2439,7 +2440,8 @@ Foam::meshRefinement::balanceAndRefine
                 << " mesh to time " << timeName() << endl;
             write
             (
-                debug,
+                debugType(debug),
+                writeType(writeLevel() | WRITEMESH),
                 mesh_.time().path()/timeName()
             );
             Pout<< "Dumped debug data in = "
@@ -2462,9 +2464,9 @@ Foam::meshRefinement::balanceAndRefine
             << " mesh to time " << timeName() << endl;
         write
         (
-            debug,
-            mesh_.time().path()
-           /timeName()
+            debugType(debug),
+            writeType(writeLevel() | WRITEMESH),
+            mesh_.time().path()/timeName()
         );
         Pout<< "Dumped debug data in = "
             << mesh_.time().cpuTimeIncrement() << " s" << endl;
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
index 83e3389a5cdab446efd255076c14773dba2230dd..45fcee4d5d793861002b045680ef7763618c8005 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
@@ -315,6 +315,21 @@ void meshRefinement::reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
 }
 
 
+template<class Enum>
+int meshRefinement::readFlags(const Enum& namedEnum, const wordList& words)
+{
+    int flags = 0;
+
+    forAll(words, i)
+    {
+        int index = namedEnum[words[i]];
+        int val = 1<<index;
+        flags |= val;
+    }
+    return flags;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/parallel/reconstruct/reconstruct/processorMeshes.C b/src/parallel/reconstruct/reconstruct/processorMeshes.C
index 875e96c38e005fb3f8c694026915c225ee9e5179..2ddbd7fe8a257fa224e64120f90394c290cb9c7e 100644
--- a/src/parallel/reconstruct/reconstruct/processorMeshes.C
+++ b/src/parallel/reconstruct/reconstruct/processorMeshes.C
@@ -31,6 +31,8 @@ License
 
 void Foam::processorMeshes::read()
 {
+    // Make sure to clear (and hence unregister) any previously loaded meshes
+    // and fields
     forAll(databases_, procI)
     {
         meshes_.set(procI, NULL);
diff --git a/src/postProcessing/functionObjects/field/Make/options b/src/postProcessing/functionObjects/field/Make/options
index de6d19ef29068e561374d756e9bb1d53718325b3..0cb1bb8f0921cd7a87be8a84e55cd06b375a2ebe 100644
--- a/src/postProcessing/functionObjects/field/Make/options
+++ b/src/postProcessing/functionObjects/field/Make/options
@@ -3,6 +3,7 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/turbulenceModels \
@@ -15,6 +16,7 @@ LIB_LIBS = \
     -lsurfMesh \
     -llagrangian \
     -lfileFormats \
+    -lsurfMesh \
     -lsampling \
     -lincompressibleTransportModels \
     -lcompressibleTurbulenceModel \
diff --git a/src/postProcessing/functionObjects/field/readFields/readFields.C b/src/postProcessing/functionObjects/field/readFields/readFields.C
index 3b28809035648096ce35aa8f321d3c0df9f80ec7..7897bfd92b1661d68cdcc04048bb1e8667be8130 100644
--- a/src/postProcessing/functionObjects/field/readFields/readFields.C
+++ b/src/postProcessing/functionObjects/field/readFields/readFields.C
@@ -120,7 +120,7 @@ void Foam::readFields::execute()
 
 void Foam::readFields::end()
 {
-    // Do nothing
+    execute();
 }