diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
index 28930c573f..39e7375a67 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
@@ -7,6 +7,7 @@
 #include "pointSet.H"
 #include "IOmanip.H"
 #include "emptyPolyPatch.H"
+#include "processorPolyPatch.H"
 
 Foam::label Foam::checkTopology
 (
@@ -289,10 +290,18 @@ Foam::label Foam::checkTopology
         }
     }
 
-    if (!Pstream::parRun())
+
     {
-        Pout<< "\nChecking patch topology for multiply connected surfaces ..."
-            << endl;
+        if (!Pstream::parRun())
+        {
+            Info<< "\nChecking patch topology for multiply connected"
+                << " surfaces..." << endl;
+        }
+        else
+        {
+            Info<< "\nChecking basic patch addressing..." << endl;
+        }
+
 
         const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
@@ -301,90 +310,101 @@ Foam::label Foam::checkTopology
         (
             mesh,
             "nonManifoldPoints",
-            mesh.nPoints()/100
+            mesh.nPoints()/1000
         );
 
         Pout.setf(ios_base::left);
 
-        Pout<< "    "
+        Info<< "    "
             << setw(20) << "Patch"
             << setw(9) << "Faces"
-            << setw(9) << "Points"
-            << setw(34) << "Surface topology";
+            << setw(9) << "Points";
+        if (!Pstream::parRun())
+        {
+            Info<< setw(34) << "Surface topology";
+        }
         if (allGeometry)
         {
-            Pout<< " Bounding box";
+            Info<< " Bounding box";
         }
-        Pout<< endl;
+        Info<< endl;
 
         forAll(patches, patchI)
         {
             const polyPatch& pp = patches[patchI];
 
-                Pout<< "    "
+            if (!isA<processorPolyPatch>(pp))
+            {
+                Info<< "    "
                     << setw(20) << pp.name()
-                    << setw(9) << pp.size()
-                    << setw(9) << pp.nPoints();
+                    << setw(9) << returnReduce(pp.size(), sumOp<label>())
+                    << setw(9) << returnReduce(pp.nPoints(), sumOp<label>());
 
-
-            primitivePatch::surfaceTopo pTyp = pp.surfaceType();
-
-            if (pp.empty())
-            {
-                Pout<< setw(34) << "ok (empty)";
-            }
-            else if (pTyp == primitivePatch::MANIFOLD)
-            {
-                if (pp.checkPointManifold(true, &points))
+                if (!Pstream::parRun())
                 {
-                    Pout<< setw(34) << "multiply connected (shared point)";
-                }
-                else
-                {
-                    Pout<< setw(34) << "ok (closed singly connected)";
-                }
+                    primitivePatch::surfaceTopo pTyp = pp.surfaceType();
 
-                // Add points on non-manifold edges to make set complete
-                pp.checkTopology(false, &points);
-            }
-            else
-            {
-                pp.checkTopology(false, &points);
-
-                if (pTyp == primitivePatch::OPEN)
-                {
-                    Pout<< setw(34) << "ok (non-closed singly connected)";
-                }
-                else
-                {
-                    Pout<< setw(34) << "multiply connected (shared edge)";
-                }
-            }
-
-            if (allGeometry)
-            {
-                const pointField& pts = pp.points();
-                const labelList& mp = pp.meshPoints();
-
-                if (returnReduce(mp.size(), sumOp<label>()) > 0)
-                {
-                    boundBox bb(point::max, point::min);
-                    forAll (mp, i)
+                    if (pp.empty())
                     {
-                        bb.min() = min(bb.min(), pts[mp[i]]);
-                        bb.max() = max(bb.max(), pts[mp[i]]);
+                        Info<< setw(34) << "ok (empty)";
+                    }
+                    else if (pTyp == primitivePatch::MANIFOLD)
+                    {
+                        if (pp.checkPointManifold(true, &points))
+                        {
+                            Info<< setw(34)
+                                << "multiply connected (shared point)";
+                        }
+                        else
+                        {
+                            Info<< setw(34) << "ok (closed singly connected)";
+                        }
+
+                        // Add points on non-manifold edges to make set complete
+                        pp.checkTopology(false, &points);
+                    }
+                    else
+                    {
+                        pp.checkTopology(false, &points);
+
+                        if (pTyp == primitivePatch::OPEN)
+                        {
+                            Info<< setw(34)
+                                << "ok (non-closed singly connected)";
+                        }
+                        else
+                        {
+                            Info<< setw(34)
+                                << "multiply connected (shared edge)";
+                        }
                     }
-                    reduce(bb.min(), minOp<vector>());
-                    reduce(bb.max(), maxOp<vector>());
-                    Pout<< ' ' << bb;
                 }
+
+                if (allGeometry)
+                {
+                    const pointField& pts = pp.points();
+                    const labelList& mp = pp.meshPoints();
+
+                    if (returnReduce(mp.size(), sumOp<label>()) > 0)
+                    {
+                        boundBox bb(point::max, point::min);
+                        forAll (mp, i)
+                        {
+                            bb.min() = min(bb.min(), pts[mp[i]]);
+                            bb.max() = max(bb.max(), pts[mp[i]]);
+                        }
+                        reduce(bb.min(), minOp<vector>());
+                        reduce(bb.max(), maxOp<vector>());
+                        Info<< ' ' << bb;
+                    }
+                }
+                Info<< endl;
             }
-            Pout<< endl;
         }
 
         if (points.size())
         {
-            Pout<< "  <<Writing " << points.size()
+            Info<< "  <<Writing " << returnReduce(points.size(), sumOp<label>())
                 << " conflicting points to set "
                 << points.name() << endl;
 
@@ -392,7 +412,7 @@ Foam::label Foam::checkTopology
             points.write();
         }
 
-        //Pout.setf(ios_base::right);
+        //Info.setf(ios_base::right);
     }
 
     // Force creation of all addressing if requested.
diff --git a/applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C b/applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C
index 3841e16bf3..98304abf41 100644
--- a/applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C
+++ b/applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,6 +115,7 @@ int main(int argc, char *argv[])
             runTimeMaster
         )
     );
+    const word oldInstance = masterMesh.pointsInstance();
 
 
     Info<< "Reading mesh to add for time = " << runTimeToAdd.timeName() << nl;
@@ -139,7 +140,13 @@ int main(int argc, char *argv[])
 
     masterMesh.addMesh(meshToAdd);
     masterMesh.merge();
-    masterMesh.polyMesh::write();
+
+    if (overwrite)
+    {
+        masterMesh.setInstance(oldInstance);
+    }
+
+    masterMesh.write();
 
     Info<< "\nEnd\n" << endl;
 
diff --git a/src/OSspecific/POSIX/printStack.C b/src/OSspecific/POSIX/printStack.C
index df31b35f7b..4a3859dcae 100644
--- a/src/OSspecific/POSIX/printStack.C
+++ b/src/OSspecific/POSIX/printStack.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -167,6 +167,26 @@ void getSymbolForRaw
 }
 
 
+void error::safePrintStack(std::ostream& os)
+{
+    // Get raw stack symbols
+    void *array[100];
+    size_t size = backtrace(array, 100);
+    char **strings = backtrace_symbols(array, size);
+
+    // See if they contain function between () e.g. "(__libc_start_main+0xd0)"
+    // and see if cplus_demangle can make sense of part before +
+    for (size_t i = 0; i < size; i++)
+    {
+        string msg(strings[i]);
+        fileName programFile;
+        word address;
+
+        os  << '#' << label(i) << '\t' << msg << std::endl;
+    }
+}
+
+
 void error::printStack(Ostream& os)
 {
     // Reads the starting addresses for the dynamically linked libraries
diff --git a/src/OpenFOAM/db/IOstreams/Sstreams/ISstream.H b/src/OpenFOAM/db/IOstreams/Sstreams/ISstream.H
index a45a784c60..9e1bb29610 100644
--- a/src/OpenFOAM/db/IOstreams/Sstreams/ISstream.H
+++ b/src/OpenFOAM/db/IOstreams/Sstreams/ISstream.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,8 +29,7 @@ Description
 
 SourceFiles
     ISstreamI.H
-    ISread.C
-    ISreadToken.C
+    ISstream.C
 
 \*---------------------------------------------------------------------------*/
 
diff --git a/src/OpenFOAM/db/error/error.H b/src/OpenFOAM/db/error/error.H
index ea599b6550..8216f5f71f 100644
--- a/src/OpenFOAM/db/error/error.H
+++ b/src/OpenFOAM/db/error/error.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -162,6 +162,10 @@ public:
         operator dictionary() const;
 
 
+        //- Helper function to print a stack (if OpenFOAM IO not yet
+        //  initialised)
+        static void safePrintStack(std::ostream&);
+
         //- Helper function to print a stack
         static void printStack(Ostream&);
 
diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/fixedValue/fixedValuePointPatchField.H b/src/OpenFOAM/fields/pointPatchFields/basic/fixedValue/fixedValuePointPatchField.H
index 5c8f1466bd..d739bc239f 100644
--- a/src/OpenFOAM/fields/pointPatchFields/basic/fixedValue/fixedValuePointPatchField.H
+++ b/src/OpenFOAM/fields/pointPatchFields/basic/fixedValue/fixedValuePointPatchField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -121,6 +121,17 @@ public:
         }
 
 
+    // Member functions
+
+        // Access
+
+            //- Return true if this patch field fixes a value
+            virtual bool fixesValue() const
+            {
+                return true;
+            }
+
+
     // Member operators
 
         // Disable assignment operators
diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H
index 5c3b19d8e7..39b1c9e551 100644
--- a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H
+++ b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -278,6 +278,12 @@ public:
                 return internalField_;
             }
 
+            //- Return true if this patch field fixes a value
+            virtual bool fixesValue() const
+            {
+                return false;
+            }
+
             //- Return true if this patch field is coupled
             virtual bool coupled() const
             {
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index 26b4fb0187..464fdd2b8f 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -710,7 +710,8 @@ void Foam::autoLayerDriver::setNumLayers
     const indirectPrimitivePatch& pp,
     pointField& patchDisp,
     labelList& patchNLayers,
-    List<extrudeMode>& extrudeStatus
+    List<extrudeMode>& extrudeStatus,
+    label& nAddedCells
 ) const
 {
     const fvMesh& mesh = meshRefiner_.mesh();
@@ -797,6 +798,24 @@ void Foam::autoLayerDriver::setNumLayers
         }
     }
 
+
+    // Calculate number of cells to create
+    nAddedCells = 0;
+    forAll(pp.localFaces(), faceI)
+    {
+        const face& f = pp.localFaces()[faceI];
+
+        // Get max of extrusion per point
+        label nCells = 0;
+        forAll(f, fp)
+        {
+            nCells = max(nCells, patchNLayers[f[fp]]);
+        }
+
+        nAddedCells += nCells;
+    }
+    reduce(nAddedCells, sumOp<label>());
+
     //reduce(nConflicts, sumOp<label>());
     //
     //Info<< "Set displacement to zero for " << nConflicts
@@ -2491,9 +2510,13 @@ void Foam::autoLayerDriver::addLayers
     // extrudeStatus = EXTRUDE.
     labelList patchNLayers(pp().nPoints(), 0);
 
+    // Ideal number of cells added
+    label nIdealAddedCells = 0;
+
     // Whether to add edge for all pp.localPoints.
     List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
 
+
     {
         // Get number of layer per point from number of layers per patch
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -2506,7 +2529,8 @@ void Foam::autoLayerDriver::addLayers
 
             patchDisp,
             patchNLayers,
-            extrudeStatus
+            extrudeStatus,
+            nIdealAddedCells
         );
 
         // Precalculate mesh edge labels for patch edges
@@ -2776,6 +2800,7 @@ void Foam::autoLayerDriver::addLayers
         layerParams.nSmoothNormals(),
         layerParams.nSmoothSurfaceNormals(),
         layerParams.minMedianAxisAngleCos(),
+        layerParams.featureAngle(),
 
         dispVec,
         medialRatio,
@@ -3101,10 +3126,24 @@ void Foam::autoLayerDriver::addLayers
 
         label nExtruded = countExtrusion(pp, extrudeStatus);
         label nTotFaces = returnReduce(pp().size(), sumOp<label>());
+        label nAddedCells = 0;
+        {
+            forAll(flaggedCells, cellI)
+            {
+                if (flaggedCells[cellI])
+                {
+                    nAddedCells++;
+                }
+            }
+            reduce(nAddedCells, sumOp<label>());
+        }
         Info<< "Extruding " << nExtruded
             << " out of " << nTotFaces
             << " faces (" << 100.0*nExtruded/nTotFaces << "%)."
             << " Removed extrusion at " << nTotChanged << " faces."
+            << endl
+            << "Added " << nAddedCells << " out of " << nIdealAddedCells
+            << " cells (" << 100.0*nAddedCells/nIdealAddedCells << "%)."
             << endl;
 
         if (nTotChanged == 0)
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
index d36d15a4cb..f0a2f5d377 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
@@ -192,7 +192,8 @@ class autoLayerDriver
                     const indirectPrimitivePatch& pp,
                     pointField& patchDisp,
                     labelList& patchNLayers,
-                    List<extrudeMode>& extrudeStatus
+                    List<extrudeMode>& extrudeStatus,
+                    label& nIdealAddedCells
                 ) const;
 
                 //- Helper function to make a pointVectorField with correct
@@ -469,6 +470,7 @@ class autoLayerDriver
                     const label nSmoothNormals,
                     const label nSmoothSurfaceNormals,
                     const scalar minMedianAxisAngleCos,
+                    const scalar featureAngle,
 
                     pointVectorField& dispVec,
                     pointScalarField& medialRatio,
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
index db246ba1f2..ba0b77effe 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
@@ -33,9 +33,10 @@ Description
 #include "motionSmoother.H"
 #include "pointData.H"
 #include "PointEdgeWave.H"
-#include "OFstream.H"
+#include "OBJstream.H"
 #include "meshTools.H"
 #include "PatchTools.H"
+#include "unitConversion.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -775,6 +776,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
     const label nSmoothNormals,
     const label nSmoothSurfaceNormals,
     const scalar minMedianAxisAngleCos,
+    const scalar featureAngle,
 
     pointVectorField& dispVec,
     pointScalarField& medialRatio,
@@ -904,7 +906,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
                 const edge& e = edges[edgeI];
                 // Approximate medial axis location on edge.
                 //const point medialAxisPt = e.centre(points);
-                vector eVec = e.vec(mesh.points());
+                vector eVec = e.vec(points);
                 scalar eMag = mag(eVec);
                 if (eMag > VSMALL)
                 {
@@ -961,42 +963,111 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
 
         labelHashSet adaptPatches(meshMover.adaptPatchIDs());
 
+
         forAll(patches, patchI)
         {
             const polyPatch& pp = patches[patchI];
-            //const pointPatchVectorField& pvf =
-            //    meshMover.displacement().boundaryField()[patchI];
+            const pointPatchVectorField& pvf =
+                meshMover.displacement().boundaryField()[patchI];
 
             if
             (
                 !pp.coupled()
              && !isA<emptyPolyPatch>(pp)
-             //&&  pvf.constraintType() != word::null //exclude fixedValue
              && !adaptPatches.found(patchI)
             )
             {
-                Info<< "Inserting points on patch " << pp.name() << endl;
-
                 const labelList& meshPoints = pp.meshPoints();
 
-                forAll(meshPoints, i)
+                // Check the type of the patchField. The types are
+                //  - fixedValue (0 or more layers) but the >0 layers have
+                //    already been handled in the adaptPatches loop
+                //  - constraint (but not coupled) types, e.g. symmetryPlane,
+                //    slip.
+                if (pvf.fixesValue())
                 {
-                    label pointI = meshPoints[i];
+                    // Disable all movement on fixedValue patchFields
+                    Info<< "Inserting all points on patch " << pp.name()
+                        << endl;
 
-                    if (!pointMedialDist[pointI].valid(dummyTrackData))
+                    forAll(meshPoints, i)
                     {
-                        maxPoints.append(pointI);
-                        maxInfo.append
-                        (
-                            pointData
+                        label pointI = meshPoints[i];
+                        if (!pointMedialDist[pointI].valid(dummyTrackData))
+                        {
+                            maxPoints.append(pointI);
+                            maxInfo.append
                             (
-                                points[pointI],
-                                0.0,
-                                pointI,         // passive data
-                                vector::zero    // passive data
-                            )
-                        );
-                        pointMedialDist[pointI] = maxInfo.last();
+                                pointData
+                                (
+                                    points[pointI],
+                                    0.0,
+                                    pointI,         // passive data
+                                    vector::zero    // passive data
+                                )
+                            );
+                            pointMedialDist[pointI] = maxInfo.last();
+                        }
+                    }
+                }
+                else
+                {
+                    // Based on geometry: analyse angle w.r.t. nearest moving
+                    // point. In the pointWallDist we transported the
+                    // normal as the passive vector. Note that this points
+                    // out of the originating wall so inside of the domain
+                    // on this patch.
+                    Info<< "Inserting points on patch " << pp.name()
+                        << " if angle to nearest layer patch > "
+                        << featureAngle << " degrees." << endl;
+
+                    scalar featureAngleCos = Foam::cos(degToRad(featureAngle));
+                    pointField pointNormals
+                    (
+                        PatchTools::pointNormals
+                        (
+                            mesh,
+                            pp,
+                            identity(pp.size())+pp.start()
+                        )
+                    );
+
+                    forAll(meshPoints, i)
+                    {
+                        label pointI = meshPoints[i];
+
+                        if (!pointMedialDist[pointI].valid(dummyTrackData))
+                        {
+                            // Check if angle not too large.
+                            scalar cosAngle =
+                            (
+                               -pointWallDist[pointI].v()
+                              & pointNormals[i]
+                            );
+                            if (cosAngle > featureAngleCos)
+                            {
+                                // Extrusion direction practically perpendicular
+                                // to the patch. Disable movement at the patch.
+
+                                maxPoints.append(pointI);
+                                maxInfo.append
+                                (
+                                    pointData
+                                    (
+                                        points[pointI],
+                                        0.0,
+                                        pointI,         // passive data
+                                        vector::zero    // passive data
+                                    )
+                                );
+                                pointMedialDist[pointI] = maxInfo.last();
+                            }
+                            else
+                            {
+                                // Extrusion direction makes angle with patch
+                                // so allow sliding - don't insert zero points
+                            }
+                        }
                     }
                 }
             }
@@ -1141,13 +1212,12 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
     // reduce thickness where thickness/medial axis distance large
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    autoPtr<OFstream> str;
-    label vertI = 0;
+    autoPtr<OBJstream> str;
     if (debug)
     {
         str.reset
         (
-            new OFstream
+            new OBJstream
             (
                 mesh.time().path()
               / "thicknessRatioExcludePoints_"
@@ -1159,13 +1229,12 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
             << str().name() << endl;
     }
 
-    autoPtr<OFstream> medialVecStr;
-    label medialVertI = 0;
+    autoPtr<OBJstream> medialVecStr;
     if (debug)
     {
         medialVecStr.reset
         (
-            new OFstream
+            new OBJstream
             (
                 mesh.time().path()
               / "thicknessRatioExcludeMedialVec_"
@@ -1227,21 +1296,19 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
                 if (str.valid())
                 {
                     const point& pt = mesh.points()[pointI];
-                    meshTools::writeOBJ(str(), pt);
-                    vertI++;
-                    meshTools::writeOBJ(str(), pt+patchDisp[patchPointI]);
-                    vertI++;
-                    str()<< "l " << vertI-1 << ' ' << vertI << nl;
+                    str().write(linePointRef(pt, pt+patchDisp[patchPointI]));
                 }
                 if (medialVecStr.valid())
                 {
                     const point& pt = mesh.points()[pointI];
-                    meshTools::writeOBJ(medialVecStr(), pt);
-                    medialVertI++;
-                    meshTools::writeOBJ(medialVecStr(), medialVec[pointI]);
-                    medialVertI++;
-                    medialVecStr()<< "l " << medialVertI-1
-                        << ' ' << medialVertI << nl;
+                    medialVecStr().write
+                    (
+                        linePointRef
+                        (
+                            pt,
+                            medialVec[pointI]
+                        )
+                    );
                 }
             }
         }
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index 9431f9b751..bd5dd86188 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -178,6 +178,24 @@ class autoSnapDriver
                     vector& edgeOffset  // offset from pt to point on edge
                 ) const;
 
+
+                //- For any reverse (so from feature back to mesh) attraction:
+                //  add attraction if diagonal points on face attracted
+                void stringFeatureEdges
+                (
+                    const label iter,
+                    const scalar featureCos,
+
+                    const indirectPrimitivePatch& pp,
+                    const scalarField& snapDist,
+
+                    const vectorField& rawPatchAttraction,
+                    const List<pointConstraint>& rawPatchConstraints,
+
+                    vectorField& patchAttraction,
+                    List<pointConstraint>& patchConstraints
+                ) const;
+
                 //- Return hit if on multiple points
                 pointIndexHit findMultiPatchPoint
                 (
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index 88ec55b1c8..7124875465 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -848,17 +848,17 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
         vector d = r.refPoint()-pt;
         d -= (d&n)*n;
 
-        // Correct for attraction to non-dominant face
-        correctAttraction
-        (
-            surfacePoints,
-            surfaceCount,
-            r.refPoint(),
-            n,                  // normalised normal
-            pt,
-
-            d                   // perpendicular offset vector
-        );
+        //// Correct for attraction to non-dominant face
+        //correctAttraction
+        //(
+        //    surfacePoints,
+        //    surfaceCount,
+        //    r.refPoint(),
+        //    n,                  // normalised normal
+        //    pt,
+        //
+        //    d                   // perpendicular offset vector
+        //);
 
         // Trim to snap distance
         if (magSqr(d) > sqr(snapDist[pointI]))
@@ -893,6 +893,15 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
         patchConstraint.applyConstraint(surfaceNormals[0]);
         patchConstraint.applyConstraint(surfaceNormals[1]);
         patchConstraint.applyConstraint(surfaceNormals[2]);
+
+        //Pout<< "# Feature point " << pt << nl;
+        //meshTools::writeOBJ(Pout, pt);
+        //meshTools::writeOBJ(Pout, surfacePoints[0]);
+        //meshTools::writeOBJ(Pout, surfacePoints[1]);
+        //meshTools::writeOBJ(Pout, surfacePoints[2]);
+        //Pout<< "l 1 2" << nl
+        //    << "l 1 3" << nl
+        //    << "l 1 4" << nl;
     }
 }
 
@@ -1001,6 +1010,195 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
 }
 
 
+void Foam::autoSnapDriver::stringFeatureEdges
+(
+    const label iter,
+    const scalar featureCos,
+
+    const indirectPrimitivePatch& pp,
+    const scalarField& snapDist,
+
+    const vectorField& rawPatchAttraction,
+    const List<pointConstraint>& rawPatchConstraints,
+
+    vectorField& patchAttraction,
+    List<pointConstraint>& patchConstraints
+) const
+{
+    // Snap edges to feature edges
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Walk existing edges and snap remaining ones (that are marked as
+    // feature edges in rawPatchConstraints)
+
+    // What this does is fill in any faces where not all points
+    // on the face are being attracted:
+    /*
+           +
+          / \
+         /   \
+      ---+    +---
+         \   /
+          \ /
+           +
+    */
+    // so the top and bottom will never get attracted since the nearest
+    // back from the feature edge will always be one of the left or right
+    // points since the face is diamond like. So here we walk the feature edges
+    // and add any non-attracted points.
+
+
+    while (true)
+    {
+        label nChanged = 0;
+
+        const labelListList& pointEdges = pp.pointEdges();
+        forAll(pointEdges, pointI)
+        {
+            if (patchConstraints[pointI].first() == 2)
+            {
+                const point& pt = pp.localPoints()[pointI];
+                const labelList& pEdges = pointEdges[pointI];
+                const vector& featVec = patchConstraints[pointI].second();
+
+                // Detect whether there are edges in both directions.
+                // (direction along the feature edge that is)
+                bool hasPos = false;
+                bool hasNeg = false;
+
+                forAll(pEdges, pEdgeI)
+                {
+                    const edge& e = pp.edges()[pEdges[pEdgeI]];
+                    label nbrPointI = e.otherVertex(pointI);
+
+                    if (patchConstraints[nbrPointI].first() > 1)
+                    {
+                        const point& nbrPt = pp.localPoints()[nbrPointI];
+                        const point featPt =
+                            nbrPt + patchAttraction[nbrPointI];
+                        const scalar cosAngle = (featVec & (featPt-pt));
+
+                        if (cosAngle > 0)
+                        {
+                            hasPos = true;
+                        }
+                        else
+                        {
+                            hasNeg = true;
+                        }
+                    }
+                }
+
+                if (!hasPos || !hasNeg)
+                {
+                    //Pout<< "**Detected feature string end at  "
+                    //    << pp.localPoints()[pointI] << endl;
+
+                    // No string. Assign best choice on either side
+                    label bestPosPointI = -1;
+                    scalar minPosDistSqr = GREAT;
+                    label bestNegPointI = -1;
+                    scalar minNegDistSqr = GREAT;
+
+                    forAll(pEdges, pEdgeI)
+                    {
+                        const edge& e = pp.edges()[pEdges[pEdgeI]];
+                        label nbrPointI = e.otherVertex(pointI);
+
+                        if
+                        (
+                            patchConstraints[nbrPointI].first() <= 1
+                         && rawPatchConstraints[nbrPointI].first() > 1
+                        )
+                        {
+                            const vector& nbrFeatVec =
+                                rawPatchConstraints[pointI].second();
+
+                            if (mag(featVec&nbrFeatVec) > featureCos)
+                            {
+                                // nbrPointI attracted to sameish feature
+                                // Note: also check on position.
+
+                                scalar d2 = magSqr
+                                (
+                                    rawPatchAttraction[nbrPointI]
+                                );
+
+                                const point featPt =
+                                    pp.localPoints()[nbrPointI]
+                                  + rawPatchAttraction[nbrPointI];
+                                const scalar cosAngle =
+                                    (featVec & (featPt-pt));
+
+                                if (cosAngle > 0)
+                                {
+                                    if (!hasPos && d2 < minPosDistSqr)
+                                    {
+                                        minPosDistSqr = d2;
+                                        bestPosPointI = nbrPointI;
+                                    }
+                                }
+                                else
+                                {
+                                    if (!hasNeg && d2 < minNegDistSqr)
+                                    {
+                                        minNegDistSqr = d2;
+                                        bestNegPointI = nbrPointI;
+                                    }
+                                }
+                            }
+                        }
+                    }
+
+                    if (bestPosPointI != -1)
+                    {
+                        // Use reconstructed-feature attraction. Use only
+                        // part of it since not sure...
+                        //const point& bestPt =
+                        //    pp.localPoints()[bestPosPointI];
+                        //Pout<< "**Overriding point " << bestPt
+                        //    << " on reconstructed feature edge at "
+                        //    << rawPatchAttraction[bestPosPointI]+bestPt
+                        //    << " to attracted-to-feature-edge." << endl;
+                        patchAttraction[bestPosPointI] =
+                            0.5*rawPatchAttraction[bestPosPointI];
+                        patchConstraints[bestPosPointI] =
+                            rawPatchConstraints[bestPosPointI];
+
+                        nChanged++;
+                    }
+                    if (bestNegPointI != -1)
+                    {
+                        // Use reconstructed-feature attraction. Use only
+                        // part of it since not sure...
+                        //const point& bestPt =
+                        //    pp.localPoints()[bestNegPointI];
+                        //Pout<< "**Overriding point " << bestPt
+                        //    << " on reconstructed feature edge at "
+                        //    << rawPatchAttraction[bestNegPointI]+bestPt
+                        //    << " to attracted-to-feature-edge." << endl;
+                        patchAttraction[bestNegPointI] =
+                            0.5*rawPatchAttraction[bestNegPointI];
+                        patchConstraints[bestNegPointI] =
+                            rawPatchConstraints[bestNegPointI];
+
+                        nChanged++;
+                    }
+                }
+            }
+        }
+
+
+        reduce(nChanged, sumOp<label>());
+        Info<< "Stringing feature edges : changed " << nChanged << " points"
+            << endl;
+        if (nChanged == 0)
+        {
+            break;
+        }
+    }
+}
+
+
 Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
 (
     const indirectPrimitivePatch& pp,
@@ -2113,156 +2311,20 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     // Walk existing edges and snap remaining ones (that are marked as
     // feature edges in allPatchConstraints)
 
-    while (true)
-    {
-        label nChanged = 0;
+    stringFeatureEdges
+    (
+        iter,
+        featureCos,
 
-        const labelListList& pointEdges = pp.pointEdges();
-        forAll(pointEdges, pointI)
-        {
-            if (patchConstraints[pointI].first() == 2)
-            {
-                const point& pt = pp.localPoints()[pointI];
-                const labelList& pEdges = pointEdges[pointI];
-                const vector& featVec = patchConstraints[pointI].second();
+        pp,
+        snapDist,
 
-                // Detect whether there are edges in both directions.
-                // (direction along the feature edge that is)
-                bool hasPos = false;
-                bool hasNeg = false;
-
-                forAll(pEdges, pEdgeI)
-                {
-                    const edge& e = pp.edges()[pEdges[pEdgeI]];
-                    label nbrPointI = e.otherVertex(pointI);
-
-                    if (patchConstraints[nbrPointI].first() > 1)
-                    {
-                        const point& nbrPt = pp.localPoints()[nbrPointI];
-                        const point featPt =
-                            nbrPt + patchAttraction[nbrPointI];
-                        const scalar cosAngle = (featVec & (featPt-pt));
-
-                        if (cosAngle > 0)
-                        {
-                            hasPos = true;
-                        }
-                        else
-                        {
-                            hasNeg = true;
-                        }
-                    }
-                }
-
-                if (!hasPos || !hasNeg)
-                {
-                    //Pout<< "**Detected feature string end at  "
-                    //    << pp.localPoints()[pointI] << endl;
-
-                    // No string. Assign best choice on either side
-                    label bestPosPointI = -1;
-                    scalar minPosDistSqr = GREAT;
-                    label bestNegPointI = -1;
-                    scalar minNegDistSqr = GREAT;
-
-                    forAll(pEdges, pEdgeI)
-                    {
-                        const edge& e = pp.edges()[pEdges[pEdgeI]];
-                        label nbrPointI = e.otherVertex(pointI);
-
-                        if
-                        (
-                            patchConstraints[nbrPointI].first() <= 1
-                         && allPatchConstraints[nbrPointI].first() > 1
-                        )
-                        {
-                            const vector& nbrFeatVec =
-                                allPatchConstraints[pointI].second();
-
-                            if (mag(featVec&nbrFeatVec) > featureCos)
-                            {
-                                // nbrPointI attracted to sameish feature
-                                // Note: also check on position.
-
-                                scalar d2 = magSqr
-                                (
-                                    allPatchAttraction[nbrPointI]
-                                );
-
-                                const point featPt =
-                                    pp.localPoints()[nbrPointI]
-                                  + allPatchAttraction[nbrPointI];
-                                const scalar cosAngle =
-                                    (featVec & (featPt-pt));
-
-                                if (cosAngle > 0)
-                                {
-                                    if (!hasPos && d2 < minPosDistSqr)
-                                    {
-                                        minPosDistSqr = d2;
-                                        bestPosPointI = nbrPointI;
-                                    }
-                                }
-                                else
-                                {
-                                    if (!hasNeg && d2 < minNegDistSqr)
-                                    {
-                                        minNegDistSqr = d2;
-                                        bestNegPointI = nbrPointI;
-                                    }
-                                }
-                            }
-                        }
-                    }
-
-                    if (bestPosPointI != -1)
-                    {
-                        // Use reconstructed-feature attraction. Use only
-                        // part of it since not sure...
-                        //const point& bestPt =
-                        //    pp.localPoints()[bestPosPointI];
-                        //Pout<< "**Overriding point " << bestPt
-                        //    << " on reconstructed feature edge at "
-                        //    << allPatchAttraction[bestPosPointI]+bestPt
-                        //    << " to attracted-to-feature-edge." << endl;
-                        patchAttraction[bestPosPointI] =
-                            0.5*allPatchAttraction[bestPosPointI];
-                        patchConstraints[bestPosPointI] =
-                            allPatchConstraints[bestPosPointI];
-
-                        nChanged++;
-                    }
-                    if (bestNegPointI != -1)
-                    {
-                        // Use reconstructed-feature attraction. Use only
-                        // part of it since not sure...
-                        //const point& bestPt =
-                        //    pp.localPoints()[bestNegPointI];
-                        //Pout<< "**Overriding point " << bestPt
-                        //    << " on reconstructed feature edge at "
-                        //    << allPatchAttraction[bestNegPointI]+bestPt
-                        //    << " to attracted-to-feature-edge." << endl;
-                        patchAttraction[bestNegPointI] =
-                            0.5*allPatchAttraction[bestNegPointI];
-                        patchConstraints[bestNegPointI] =
-                            allPatchConstraints[bestNegPointI];
-
-                        nChanged++;
-                    }
-                }
-            }
-        }
-
-
-        reduce(nChanged, sumOp<label>());
-        Info<< "Stringing feature edges : changed " << nChanged << " points"
-            << endl;
-        if (nChanged == 0)
-        {
-            break;
-        }
-    }
+        allPatchAttraction,
+        allPatchConstraints,
 
+        patchAttraction,
+        patchConstraints
+    );
 
 
     // Avoid diagonal attraction
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index e1b76c44e1..8cc7ba12f2 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -312,7 +312,8 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
                     {
                         Pout<< "Adding particle from point:" << pointI
                             << " coord:" << featureMesh.points()[pointI]
-                            << " pEdges:" << pointEdges[pointI]
+                            << " since number of emanating edges:"
+                            << pointEdges[pointI].size()
                             << endl;
                     }
 
diff --git a/src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermo.C b/src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermo.C
index 00f75ac7cc..b0f985f2f3 100644
--- a/src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermo.C
+++ b/src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermo.C
@@ -31,7 +31,7 @@ License
 namespace Foam
 {
     defineTypeNameAndDebug(solidReactionThermo, 0);
-    defineRunTimeSelectionTable(solidReactionThermo, mesh);
+    defineRunTimeSelectionTable(solidReactionThermo, fvMesh);
     defineRunTimeSelectionTable(solidReactionThermo, dictionary);
 }
 
diff --git a/src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermo.H b/src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermo.H
index 6aa8a76047..e9112602e1 100644
--- a/src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermo.H
+++ b/src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermo.H
@@ -66,7 +66,7 @@ public:
     (
         autoPtr,
         solidReactionThermo,
-        mesh,
+        fvMesh,
         (const fvMesh& mesh),
         (mesh)
     );
diff --git a/src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermoNew.C b/src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermoNew.C
index 6ef34f0e76..38d044c9f1 100644
--- a/src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermoNew.C
+++ b/src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermoNew.C
@@ -53,16 +53,16 @@ Foam::autoPtr<Foam::solidReactionThermo> Foam::solidReactionThermo::New
 
     Info<< "Selecting thermodynamics package " << modelType << endl;
 
-    meshConstructorTable::iterator cstrIter =
-        meshConstructorTablePtr_->find(modelType);
+    fvMeshConstructorTable::iterator cstrIter =
+        fvMeshConstructorTablePtr_->find(modelType);
 
-    if (cstrIter == meshConstructorTablePtr_->end())
+    if (cstrIter == fvMeshConstructorTablePtr_->end())
     {
         FatalErrorIn("solidReactionThermo::New(const fvMesh&)")
             << "Unknown solidReactionThermo type "
             << modelType << nl << nl
             << "Valid solidReactionThermo types:" << nl
-            << meshConstructorTablePtr_->sortedToc() << nl
+            << fvMeshConstructorTablePtr_->sortedToc() << nl
             << exit(FatalError);
     }