From 202d7c1ad3d0b40d4484c9b0a6a5f73027584a67 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 12 Nov 2013 09:00:04 +0000
Subject: [PATCH] ENH: snappyHexMesh: fix growing of attraction. Split off
 debug/writing/output

---
 .../generation/snappyHexMesh/snappyHexMesh.C  | 126 ++++--
 .../snappyHexMesh/snappyHexMeshDict           |  24 +-
 .../autoHexMeshDriver/autoLayerDriver.C       | 387 ++++++++++++------
 .../autoHexMeshDriver/autoLayerDriver.H       |  10 +
 .../autoHexMeshDriver/autoLayerDriverShrink.C |  14 +-
 .../autoHexMeshDriver/autoRefineDriver.C      |  29 +-
 .../autoHexMeshDriver/autoSnapDriver.C        |  37 +-
 .../autoHexMeshDriver/autoSnapDriver.H        |  14 +-
 .../autoHexMeshDriver/autoSnapDriverFeature.C | 360 ++++++++--------
 .../meshRefinement/meshRefinement.C           |  86 +++-
 .../meshRefinement/meshRefinement.H           |  84 +++-
 .../meshRefinement/meshRefinementBaffles.C    |  42 +-
 .../meshRefinementProblemCells.C              |   7 +-
 .../meshRefinement/meshRefinementRefine.C     |  18 +-
 .../meshRefinement/meshRefinementTemplates.C  |  15 +
 15 files changed, 854 insertions(+), 399 deletions(-)

diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index c349740ea7..b0cf1879c1 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;
     }
 
-    const bool writeLevel = meshDict.lookupOrDefault<bool>("writeLevel", false);
+    // Set file writing level
+    {
+        wordList flags;
+        if (meshDict.readIfPresent("writeFlags", flags))
+        {
+            meshRefinement::writeLevel
+            (
+                meshRefinement::writeType
+                (
+                    meshRefinement::readFlags
+                    (
+                        meshRefinement::IOwriteTypeNames,
+                        flags
+                    )
+                )
+            );
+        }
+    }
+
+    // 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 b850033a44..507891d356 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/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index 7bb72786b8..316d376def 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -2472,6 +2472,99 @@ void Foam::autoLayerDriver::getLayerCellsFaces
 }
 
 
+void Foam::autoLayerDriver::printLayerData
+(
+    const fvMesh& mesh,
+    const labelList& patchIDs,
+    const labelList& cellNLayers,
+    const scalarField& faceWantedThickness,
+    const scalarField& faceRealThickness
+) const
+{
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+    // Find maximum length of a patch name, for a nicer output
+    label maxPatchNameLen = 0;
+    forAll(patchIDs, i)
+    {
+        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)
+        {
+            sumNLayers += cellNLayers[faceCells[i]];
+        }
+
+        // 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 (patchWanted[i] > VSMALL)
+            {
+                sumFraction += (patchReal[i]/patchWanted[i]);
+            }
+        }
+
+
+        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,
@@ -2483,151 +2576,189 @@ bool Foam::autoLayerDriver::writeLayerData
 {
     bool allOk = true;
 
+    if (meshRefinement::writeLevel() & meshRefinement::WRITELAYERSETS)
     {
-        label nAdded = 0;
-        forAll(cellNLayers, cellI)
         {
-            if (cellNLayers[cellI] > 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;
         }
-        cellSet addedCellSet(mesh, "addedCells", nAdded);
-        forAll(cellNLayers, cellI)
         {
-            if (cellNLayers[cellI] > 0)
+            label nAdded = 0;
+            forAll(faceRealThickness, faceI)
             {
-                addedCellSet.insert(cellI);
+                if (faceRealThickness[faceI] > 0)
+                {
+                    nAdded++;
+                }
             }
-        }
-        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;
-    }
-    {
-        label nAdded = 0;
-        forAll(faceRealThickness, faceI)
-        {
-            if (faceRealThickness[faceI] > 0)
-            {
-                nAdded++;
-            }
-        }
 
-        faceSet layerFacesSet(mesh, "layerFaces", nAdded);
-        forAll(faceRealThickness, faceI)
-        {
-            if (faceRealThickness[faceI] > 0)
+            faceSet layerFacesSet(mesh, "layerFaces", nAdded);
+            forAll(faceRealThickness, faceI)
             {
-                layerFacesSet.insert(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
-            (
-                "nSurfaceLayers",
-                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];
-            const polyPatch& pp = pbm[patchI];
-            const labelList& faceCells = pp.faceCells();
-            scalarField pfld(faceCells.size());
-            forAll(faceCells, i)
+            volScalarField fld
+            (
+                IOobject
+                (
+                    "nSurfaceLayers",
+                    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)
             {
-                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
-            (
-                "wantedThickness",
-                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];
-            fld.boundaryField()[patchI] == pbm[patchI].patchSlice
+            volScalarField fld
             (
-                faceWantedThickness
-            );
-        }
-        Info<< "Writing volScalarField " << fld.name()
-            << " with wanted thickness" << endl;
-        bool ok = fld.write();
-        allOk = allOk & ok;
-    }
-    {
-        volScalarField fld
-        (
-            IOobject
-            (
-                "thickness",
-                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)
-        {
-            label patchI = patchIDs[i];
-            fld.boundaryField()[patchI] == pbm[patchI].patchSlice
-            (
-                faceRealThickness
+                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;
+        }
+        {
+            volScalarField fld
+            (
+                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 layer thickness" << endl;
-        bool ok = fld.write();
-        allOk = allOk & ok;
     }
+
+    //if (meshRefinement::outputLevel() & meshRefinement::OUTPUTLAYERINFO)
+    {
+        printLayerData
+        (
+            mesh,
+            patchIDs,
+            cellNLayers,
+            faceWantedThickness,
+            faceRealThickness
+        );
+    }
+
     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 4534ebe00d..58464ad495 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 3e824b1e02..328b4fd86a 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 a53b6f11eb..ce97c27104 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 adb7c38590..03017d8f8c 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 b8282bedfc..496f980a11 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 32133ad5f8..92ef64bf45 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 2406b84310..b8609ede3e 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 873167bc4f..168cd487fe 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 491b9f49d4..8f4c504c00 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 a1c1ed0f58..f49ad65a32 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 d2556c4a9d..83ed367e89 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 83e3389a5c..45fcee4d5d 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