From 3e71574d8afd4686dedd5153527dd484a7d36467 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 5 Apr 2011 11:56:57 +0100
Subject: [PATCH] ENH: snappyHexMesh: initial feature-line support

---
 .../generation/snappyHexMesh/snappyHexMesh.C  |   26 +-
 .../snappyHexMesh/snappyHexMeshDict           |    8 +-
 src/mesh/autoMesh/Make/files                  |    2 +
 .../autoHexMeshDriver/autoLayerDriver.C       |  638 +----
 .../autoHexMeshDriver/autoLayerDriver.H       |   40 -
 .../autoHexMeshDriver/autoRefineDriver.C      |  101 +-
 .../autoHexMeshDriver/autoRefineDriver.H      |   15 +-
 .../autoHexMeshDriver/autoSnapDriver.C        |  188 +-
 .../autoHexMeshDriver/autoSnapDriver.H        |  197 +-
 .../autoHexMeshDriver/autoSnapDriverFeature.C | 2146 +++++++++++++++++
 .../snapParameters/snapParameters.C           |   15 +-
 .../snapParameters/snapParameters.H           |   15 +-
 .../meshRefinement/meshRefinement.C           |   55 +-
 .../meshRefinement/meshRefinement.H           |   72 +-
 .../meshRefinement/meshRefinementBaffles.C    |   34 +-
 .../meshRefinement/meshRefinementMerge.C      |  895 ++++++-
 .../meshRefinementProblemCells.C              |   14 +-
 .../meshRefinement/meshRefinementRefine.C     |   56 +-
 .../refinementFeatures/refinementFeatures.C   |  233 ++
 .../refinementFeatures/refinementFeatures.H   |  139 ++
 20 files changed, 3835 insertions(+), 1054 deletions(-)
 create mode 100644 src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
 create mode 100644 src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
 create mode 100644 src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H

diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index 5691f01034..5b36e99301 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,6 +37,7 @@ Description
 #include "autoLayerDriver.H"
 #include "searchableSurfaces.H"
 #include "refinementSurfaces.H"
+#include "refinementFeatures.H"
 #include "shellSurfaces.H"
 #include "decompositionMethod.H"
 #include "fvMeshDistribute.H"
@@ -250,6 +251,20 @@ int main(int argc, char *argv[])
         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
 
 
+    // Read feature meshes
+    // ~~~~~~~~~~~~~~~~~~~
+
+    Info<< "Reading features." << endl;
+    refinementFeatures features
+    (
+        mesh,
+        refineDict.lookup("features")
+    );
+    Info<< "Read features in = "
+        << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
+
+
+
     // Refinement engine
     // ~~~~~~~~~~~~~~~~~
 
@@ -265,6 +280,7 @@ int main(int argc, char *argv[])
         mergeDist,          // tolerance used in sorting coordinates
         overwrite,          // overwrite mesh files?
         surfaces,           // for surface intersection refinement
+        features,           // for feature edges/point based refinement
         shells              // for volume (inside/outside) refinement
     );
     Info<< "Calculated surface intersections in = "
@@ -429,13 +445,19 @@ int main(int argc, char *argv[])
 
         // Snap parameters
         snapParameters snapParams(snapDict);
+        // Temporary hack to get access to resolveFeatureAngle
+        scalar curvature;
+        {
+            refinementParameters refineParams(refineDict);
+            curvature = refineParams.curvature();
+        }
 
         if (!overwrite)
         {
             const_cast<Time&>(mesh.time())++;
         }
 
-        snapDriver.doSnap(snapDict, motionDict, snapParams);
+        snapDriver.doSnap(snapDict, motionDict, curvature, snapParams);
 
         writeMesh
         (
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index 8ee2c4b74c..da7a8ac63c 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -217,10 +217,10 @@ snapControls
     //  to surface
     nSmoothPatch 3;
 
-    //- Relative distance for points to be attracted by surface feature point
-    //  or edge. True distance is this factor times local
-    //  maximum edge length.
-    tolerance 4.0;
+    //- Maximum relative distance for points to be attracted by surface.
+    //  True distance is this factor times local maximum edge length.
+    // Note: changed(corrected) w.r.t 17x! (17x used 2* tolerance)
+    tolerance 2.0;
 
     //- Number of mesh displacement relaxation iterations.
     nSolveIter 30;
diff --git a/src/mesh/autoMesh/Make/files b/src/mesh/autoMesh/Make/files
index ec877c5f94..66cb3c6202 100644
--- a/src/mesh/autoMesh/Make/files
+++ b/src/mesh/autoMesh/Make/files
@@ -4,6 +4,7 @@ autoHexMeshDriver       = $(autoHexMesh)/autoHexMeshDriver
 $(autoHexMeshDriver)/autoLayerDriver.C
 $(autoHexMeshDriver)/autoLayerDriverShrink.C
 $(autoHexMeshDriver)/autoSnapDriver.C
+$(autoHexMeshDriver)/autoSnapDriverFeature.C
 $(autoHexMeshDriver)/autoRefineDriver.C
 
 $(autoHexMeshDriver)/layerParameters/layerParameters.C
@@ -16,6 +17,7 @@ $(autoHexMesh)/meshRefinement/meshRefinement.C
 $(autoHexMesh)/meshRefinement/meshRefinementMerge.C
 $(autoHexMesh)/meshRefinement/meshRefinementProblemCells.C
 $(autoHexMesh)/meshRefinement/meshRefinementRefine.C
+$(autoHexMesh)/refinementFeatures/refinementFeatures.C
 $(autoHexMesh)/refinementSurfaces/refinementSurfaces.C
 $(autoHexMesh)/shellSurfaces/shellSurfaces.C
 $(autoHexMesh)/trackedParticle/trackedParticle.C
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index ca38821229..cdaa9f83c7 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -59,622 +59,6 @@ defineTypeNameAndDebug(autoLayerDriver, 0);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::label Foam::autoLayerDriver::mergePatchFacesUndo
-(
-    const scalar minCos,
-    const scalar concaveCos,
-    const dictionary& motionDict
-)
-{
-    fvMesh& mesh = meshRefiner_.mesh();
-
-    // Patch face merging engine
-    combineFaces faceCombiner(mesh, true);
-
-    // Pick up all candidate cells on boundary
-    labelHashSet boundaryCells(mesh.nFaces()-mesh.nInternalFaces());
-
-    {
-        labelList patchIDs(meshRefiner_.meshedPatches());
-
-        const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-        forAll(patchIDs, i)
-        {
-            label patchI = patchIDs[i];
-
-            const polyPatch& patch = patches[patchI];
-
-            if (!patch.coupled())
-            {
-                forAll(patch, i)
-                {
-                    boundaryCells.insert(mesh.faceOwner()[patch.start()+i]);
-                }
-            }
-        }
-    }
-
-    // Get all sets of faces that can be merged
-    labelListList allFaceSets
-    (
-        faceCombiner.getMergeSets
-        (
-            minCos,
-            concaveCos,
-            boundaryCells
-        )
-    );
-
-    label nFaceSets = returnReduce(allFaceSets.size(), sumOp<label>());
-
-    Info<< "Merging " << nFaceSets << " sets of faces." << nl << endl;
-
-    if (nFaceSets > 0)
-    {
-        if (debug)
-        {
-            faceSet allSets(mesh, "allFaceSets", allFaceSets.size());
-            forAll(allFaceSets, setI)
-            {
-                forAll(allFaceSets[setI], i)
-                {
-                    allSets.insert(allFaceSets[setI][i]);
-                }
-            }
-            Pout<< "Writing all faces to be merged to set "
-                << allSets.objectPath() << endl;
-            allSets.instance() = mesh.time().timeName();
-            allSets.write();
-        }
-
-
-        // Topology changes container
-        polyTopoChange meshMod(mesh);
-
-        // Merge all faces of a set into the first face of the set.
-        faceCombiner.setRefinement(allFaceSets, meshMod);
-
-        // Experimental: store data for all the points that have been deleted
-        meshRefiner_.storeData
-        (
-            faceCombiner.savedPointLabels(),    // points to store
-            labelList(0),                       // faces to store
-            labelList(0)                        // cells to store
-        );
-
-        // Change the mesh (no inflation)
-        autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
-
-        // Update fields
-        mesh.updateMesh(map);
-
-        // Move mesh (since morphing does not do this)
-        if (map().hasMotionPoints())
-        {
-            mesh.movePoints(map().preMotionPoints());
-        }
-        else
-        {
-            // Delete mesh volumes.
-            mesh.clearOut();
-        }
-
-        faceCombiner.updateMesh(map);
-
-        meshRefiner_.updateMesh(map, labelList(0));
-
-
-        for (label iteration = 0; iteration < 100; iteration++)
-        {
-            Info<< nl
-                << "Undo iteration " << iteration << nl
-                << "----------------" << endl;
-
-
-            // Check mesh for errors
-            // ~~~~~~~~~~~~~~~~~~~~~
-
-            faceSet errorFaces
-            (
-                mesh,
-                "errorFaces",
-                mesh.nFaces()-mesh.nInternalFaces()
-            );
-            bool hasErrors = motionSmoother::checkMesh
-            (
-                false,  // report
-                mesh,
-                motionDict,
-                errorFaces
-            );
-
-            //if (checkEdgeConnectivity)
-            //{
-            //    Info<< "Checking edge-face connectivity (duplicate faces"
-            //        << " or non-consecutive shared vertices)" << endl;
-            //
-            //    label nOldSize = errorFaces.size();
-            //
-            //    hasErrors =
-            //        mesh.checkFaceFaces
-            //        (
-            //            false,
-            //            &errorFaces
-            //        )
-            //     || hasErrors;
-            //
-            //    Info<< "Detected additional "
-            //        <<  returnReduce
-            //            (
-            //                errorFaces.size() - nOldSize,
-            //                sumOp<label>()
-            //            )
-            //        << " faces with illegal face-face connectivity"
-            //        << endl;
-            //}
-
-            if (!hasErrors)
-            {
-                break;
-            }
-
-
-            if (debug)
-            {
-                errorFaces.instance() = mesh.time().timeName();
-                Pout<< "Writing all faces in error to faceSet "
-                    << errorFaces.objectPath() << nl << endl;
-                errorFaces.write();
-            }
-
-
-            // Check any master cells for using any of the error faces
-            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-            DynamicList<label> mastersToRestore(allFaceSets.size());
-
-            forAll(allFaceSets, setI)
-            {
-                label masterFaceI = faceCombiner.masterFace()[setI];
-
-                if (masterFaceI != -1)
-                {
-                    label masterCellII = mesh.faceOwner()[masterFaceI];
-
-                    const cell& cFaces = mesh.cells()[masterCellII];
-
-                    forAll(cFaces, i)
-                    {
-                        if (errorFaces.found(cFaces[i]))
-                        {
-                            mastersToRestore.append(masterFaceI);
-                            break;
-                        }
-                    }
-                }
-            }
-            mastersToRestore.shrink();
-
-            label nRestore = returnReduce
-            (
-                mastersToRestore.size(),
-                sumOp<label>()
-            );
-
-            Info<< "Masters that need to be restored:"
-                << nRestore << endl;
-
-            if (debug)
-            {
-                faceSet restoreSet(mesh, "mastersToRestore", mastersToRestore);
-                restoreSet.instance() = mesh.time().timeName();
-                Pout<< "Writing all " << mastersToRestore.size()
-                    << " masterfaces to be restored to set "
-                    << restoreSet.objectPath() << endl;
-                restoreSet.write();
-            }
-
-
-            if (nRestore == 0)
-            {
-                break;
-            }
-
-
-            // Undo
-            // ~~~~
-
-            // Topology changes container
-            polyTopoChange meshMod(mesh);
-
-            // Merge all faces of a set into the first face of the set.
-            // Experimental:mark all points/faces/cells that have been restored.
-            Map<label> restoredPoints(0);
-            Map<label> restoredFaces(0);
-            Map<label> restoredCells(0);
-
-            faceCombiner.setUnrefinement
-            (
-                mastersToRestore,
-                meshMod,
-                restoredPoints,
-                restoredFaces,
-                restoredCells
-            );
-
-            // Change the mesh (no inflation)
-            autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
-
-            // Update fields
-            mesh.updateMesh(map);
-
-            // Move mesh (since morphing does not do this)
-            if (map().hasMotionPoints())
-            {
-                mesh.movePoints(map().preMotionPoints());
-            }
-            else
-            {
-                // Delete mesh volumes.
-                mesh.clearOut();
-            }
-
-            faceCombiner.updateMesh(map);
-
-            // Renumber restore maps
-            inplaceMapKey(map().reversePointMap(), restoredPoints);
-            inplaceMapKey(map().reverseFaceMap(), restoredFaces);
-            inplaceMapKey(map().reverseCellMap(), restoredCells);
-
-            // Experimental:restore all points/face/cells in maps
-            meshRefiner_.updateMesh
-            (
-                map,
-                labelList(0),       // changedFaces
-                restoredPoints,
-                restoredFaces,
-                restoredCells
-            );
-
-            Info<< endl;
-        }
-
-        if (debug)
-        {
-            Pout<< "Writing merged-faces mesh to time "
-                << meshRefiner_.timeName() << nl << endl;
-            mesh.write();
-        }
-    }
-    else
-    {
-        Info<< "No faces merged ..." << endl;
-    }
-
-    return nFaceSets;
-}
-
-
-// Remove points. pointCanBeDeleted is parallel synchronised.
-Foam::autoPtr<Foam::mapPolyMesh> Foam::autoLayerDriver::doRemovePoints
-(
-    removePoints& pointRemover,
-    const boolList& pointCanBeDeleted
-)
-{
-    fvMesh& mesh = meshRefiner_.mesh();
-
-    // Topology changes container
-    polyTopoChange meshMod(mesh);
-
-    pointRemover.setRefinement(pointCanBeDeleted, meshMod);
-
-    // Change the mesh (no inflation)
-    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
-
-    // Update fields
-    mesh.updateMesh(map);
-
-    // Move mesh (since morphing does not do this)
-    if (map().hasMotionPoints())
-    {
-        mesh.movePoints(map().preMotionPoints());
-    }
-    else
-    {
-        // Delete mesh volumes.
-        mesh.clearOut();
-    }
-
-    pointRemover.updateMesh(map);
-    meshRefiner_.updateMesh(map, labelList(0));
-
-    return map;
-}
-
-
-// Restore faces (which contain removed points)
-Foam::autoPtr<Foam::mapPolyMesh> Foam::autoLayerDriver::doRestorePoints
-(
-    removePoints& pointRemover,
-    const labelList& facesToRestore
-)
-{
-    fvMesh& mesh = meshRefiner_.mesh();
-
-    // Topology changes container
-    polyTopoChange meshMod(mesh);
-
-    // Determine sets of points and faces to restore
-    labelList localFaces, localPoints;
-    pointRemover.getUnrefimentSet
-    (
-        facesToRestore,
-        localFaces,
-        localPoints
-    );
-
-    // Undo the changes on the faces that are in error.
-    pointRemover.setUnrefinement
-    (
-        localFaces,
-        localPoints,
-        meshMod
-    );
-
-    // Change the mesh (no inflation)
-    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
-
-    // Update fields
-    mesh.updateMesh(map);
-
-    // Move mesh (since morphing does not do this)
-    if (map().hasMotionPoints())
-    {
-        mesh.movePoints(map().preMotionPoints());
-    }
-    else
-    {
-        // Delete mesh volumes.
-        mesh.clearOut();
-    }
-
-    pointRemover.updateMesh(map);
-    meshRefiner_.updateMesh(map, labelList(0));
-
-    return map;
-}
-
-
-// Collect all faces that are both in candidateFaces and in the set.
-// If coupled face also collects the coupled face.
-Foam::labelList Foam::autoLayerDriver::collectFaces
-(
-    const labelList& candidateFaces,
-    const labelHashSet& set
-) const
-{
-    const fvMesh& mesh = meshRefiner_.mesh();
-
-    // Has face been selected?
-    boolList selected(mesh.nFaces(), false);
-
-    forAll(candidateFaces, i)
-    {
-        label faceI = candidateFaces[i];
-
-        if (set.found(faceI))
-        {
-            selected[faceI] = true;
-        }
-    }
-    syncTools::syncFaceList
-    (
-        mesh,
-        selected,
-        orEqOp<bool>()      // combine operator
-    );
-
-    labelList selectedFaces(findIndices(selected, true));
-
-    return selectedFaces;
-}
-
-
-// Pick up faces of cells of faces in set.
-Foam::labelList Foam::autoLayerDriver::growFaceCellFace
-(
-    const labelHashSet& set
-) const
-{
-    const fvMesh& mesh = meshRefiner_.mesh();
-
-    boolList selected(mesh.nFaces(), false);
-
-    forAllConstIter(faceSet, set, iter)
-    {
-        label faceI = iter.key();
-
-        label own = mesh.faceOwner()[faceI];
-
-        const cell& ownFaces = mesh.cells()[own];
-        forAll(ownFaces, ownFaceI)
-        {
-            selected[ownFaces[ownFaceI]] = true;
-        }
-
-        if (mesh.isInternalFace(faceI))
-        {
-            label nbr = mesh.faceNeighbour()[faceI];
-
-            const cell& nbrFaces = mesh.cells()[nbr];
-            forAll(nbrFaces, nbrFaceI)
-            {
-                selected[nbrFaces[nbrFaceI]] = true;
-            }
-        }
-    }
-    syncTools::syncFaceList
-    (
-        mesh,
-        selected,
-        orEqOp<bool>()      // combine operator
-    );
-    return findIndices(selected, true);
-}
-
-
-// Remove points not used by any face or points used by only two faces where
-// the edges are in line
-Foam::label Foam::autoLayerDriver::mergeEdgesUndo
-(
-    const scalar minCos,
-    const dictionary& motionDict
-)
-{
-    fvMesh& mesh = meshRefiner_.mesh();
-
-    Info<< nl
-        << "Merging all points on surface that" << nl
-        << "- are used by only two boundary faces and" << nl
-        << "- make an angle with a cosine of more than " << minCos
-        << "." << nl << endl;
-
-    // Point removal analysis engine with undo
-    removePoints pointRemover(mesh, true);
-
-    // Count usage of points
-    boolList pointCanBeDeleted;
-    label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
-
-    if (nRemove > 0)
-    {
-        Info<< "Removing " << nRemove
-            << " straight edge points ..." << nl << endl;
-
-        // Remove points
-        // ~~~~~~~~~~~~~
-
-        doRemovePoints(pointRemover, pointCanBeDeleted);
-
-
-        for (label iteration = 0; iteration < 100; iteration++)
-        {
-            Info<< nl
-                << "Undo iteration " << iteration << nl
-                << "----------------" << endl;
-
-
-            // Check mesh for errors
-            // ~~~~~~~~~~~~~~~~~~~~~
-
-            faceSet errorFaces
-            (
-                mesh,
-                "errorFaces",
-                mesh.nFaces()-mesh.nInternalFaces()
-            );
-            bool hasErrors = motionSmoother::checkMesh
-            (
-                false,  // report
-                mesh,
-                motionDict,
-                errorFaces
-            );
-            //if (checkEdgeConnectivity)
-            //{
-            //    Info<< "Checking edge-face connectivity (duplicate faces"
-            //        << " or non-consecutive shared vertices)" << endl;
-            //
-            //    label nOldSize = errorFaces.size();
-            //
-            //    hasErrors =
-            //        mesh.checkFaceFaces
-            //        (
-            //            false,
-            //            &errorFaces
-            //        )
-            //     || hasErrors;
-            //
-            //    Info<< "Detected additional "
-            //        << returnReduce(errorFaces.size()-nOldSize,sumOp<label>())
-            //        << " faces with illegal face-face connectivity"
-            //        << endl;
-            //}
-
-            if (!hasErrors)
-            {
-                break;
-            }
-
-            if (debug)
-            {
-                errorFaces.instance() = mesh.time().timeName();
-                Pout<< "**Writing all faces in error to faceSet "
-                    << errorFaces.objectPath() << nl << endl;
-                errorFaces.write();
-            }
-
-            labelList masterErrorFaces
-            (
-                collectFaces
-                (
-                    pointRemover.savedFaceLabels(),
-                    errorFaces
-                )
-            );
-
-            label n = returnReduce(masterErrorFaces.size(), sumOp<label>());
-
-            Info<< "Detected " << n
-                << " error faces on boundaries that have been merged."
-                << " These will be restored to their original faces." << nl
-                << endl;
-
-            if (n == 0)
-            {
-                if (hasErrors)
-                {
-                    Info<< "Detected "
-                        << returnReduce(errorFaces.size(), sumOp<label>())
-                        << " error faces in mesh."
-                        << " Restoring neighbours of faces in error." << nl
-                        << endl;
-
-                    labelList expandedErrorFaces
-                    (
-                        growFaceCellFace
-                        (
-                            errorFaces
-                        )
-                    );
-
-                    doRestorePoints(pointRemover, expandedErrorFaces);
-                }
-
-                break;
-            }
-
-            doRestorePoints(pointRemover, masterErrorFaces);
-        }
-
-        if (debug)
-        {
-            Pout<< "Writing merged-edges mesh to time "
-                << meshRefiner_.timeName() << nl << endl;
-            mesh.write();
-        }
-    }
-    else
-    {
-        Info<< "No straight edges simplified and no points removed ..." << endl;
-    }
-
-    return nRemove;
-}
-
-
 // For debugging: Dump displacement to .obj files
 void Foam::autoLayerDriver::dumpDisplacement
 (
@@ -793,7 +177,7 @@ void Foam::autoLayerDriver::checkMeshManifold() const
             << " points where this happens to pointSet "
             << nonManifoldPoints.name() << endl;
 
-        nonManifoldPoints.instance() = mesh.time().timeName();
+        nonManifoldPoints.instance() = meshRefiner_.timeName();
         nonManifoldPoints.write();
     }
     Info<< endl;
@@ -1080,7 +464,6 @@ void Foam::autoLayerDriver::handleWarpedFaces
 }
 
 
-
 //// No extrusion on cells with multiple patch faces. There ususally is a reason
 //// why combinePatchFaces hasn't succeeded.
 //void Foam::autoLayerDriver::handleMultiplePatchFaces
@@ -1133,7 +516,7 @@ void Foam::autoLayerDriver::handleWarpedFaces
 //
 //    if (nMultiPatchCells > 0)
 //    {
-//        multiPatchCells.instance() = mesh.time().timeName();
+//        multiPatchCells.instance() = meshRefiner_.timeName();
 //        Info<< "Writing " << nMultiPatchCells
 //            << " cells with multiple (connected) patch faces to cellSet "
 //            << multiPatchCells.objectPath() << endl;
@@ -2547,9 +1930,15 @@ void Foam::autoLayerDriver::mergePatchFacesUndo
         << "      (0=straight, 180=fully concave)" << nl
         << endl;
 
-    label nChanged = mergePatchFacesUndo(minCos, concaveCos, motionDict);
+    label nChanged = meshRefiner_.mergePatchFacesUndo
+    (
+        minCos,
+        concaveCos,
+        meshRefiner_.meshedPatches(),
+        motionDict
+    );
 
-    nChanged += mergeEdgesUndo(minCos, motionDict);
+    nChanged += meshRefiner_.mergeEdgesUndo(minCos, motionDict);
 }
 
 
@@ -3191,6 +2580,7 @@ void Foam::autoLayerDriver::addLayers
                 "addedCells",
                 findIndices(flaggedCells, true)
             );
+            addedCellSet.instance() = meshRefiner_.timeName();
             Info<< "Writing "
                 << returnReduce(addedCellSet.size(), sumOp<label>())
                 << " added cells to cellSet "
@@ -3203,6 +2593,7 @@ void Foam::autoLayerDriver::addLayers
                 "layerFaces",
                 findIndices(flaggedCells, true)
             );
+            layerFacesSet.instance() = meshRefiner_.timeName();
             Info<< "Writing "
                 << returnReduce(layerFacesSet.size(), sumOp<label>())
                 << " faces inside added layer to faceSet "
@@ -3279,6 +2670,9 @@ void Foam::autoLayerDriver::addLayers
         mesh.clearOut();
     }
 
+    // Reset the instance for if in overwrite mode
+    mesh.setInstance(meshRefiner_.timeName());
+
     meshRefiner_.updateMesh(map, labelList(0));
 
 
@@ -3345,6 +2739,7 @@ void Foam::autoLayerDriver::addLayers
     // ~~~~~~~~~~
 
     cellSet addedCellSet(mesh, "addedCells", findIndices(flaggedCells, true));
+    addedCellSet.instance() = meshRefiner_.timeName();
     Info<< "Writing "
         << returnReduce(addedCellSet.size(), sumOp<label>())
         << " added cells to cellSet "
@@ -3352,6 +2747,7 @@ void Foam::autoLayerDriver::addLayers
     addedCellSet.write();
 
     faceSet layerFacesSet(mesh, "layerFaces", findIndices(flaggedFaces, true));
+    layerFacesSet.instance() = meshRefiner_.timeName();
     Info<< "Writing "
         << returnReduce(layerFacesSet.size(), sumOp<label>())
         << " faces inside added layer to faceSet "
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
index fd5a1f3aa8..ff3586c422 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
@@ -105,46 +105,6 @@ class autoLayerDriver
 
     // Private Member Functions
 
-
-        // Face merging
-
-            //- Merge patch faces. Undo until no checkMesh errors.
-            label mergePatchFacesUndo
-            (
-                const scalar minCos,
-                const scalar concaveCos,
-                const dictionary&
-            );
-
-            //- Remove points.
-            autoPtr<mapPolyMesh> doRemovePoints
-            (
-                removePoints& pointRemover,
-                const boolList& pointCanBeDeleted
-            );
-
-            //- Restore faces (which contain removed points)
-            autoPtr<mapPolyMesh> doRestorePoints
-            (
-                removePoints& pointRemover,
-                const labelList& facesToRestore
-            );
-
-            //- Return candidateFaces that are also in set.
-            labelList collectFaces
-            (
-                const labelList& candidateFaces,
-                const labelHashSet& set
-            ) const;
-
-            //- Pick up faces of cells of faces in set.
-            labelList growFaceCellFace(const labelHashSet&) const;
-
-            //- Remove points not used by any face or points used by only
-            //  two faces where the edges are in line
-            label mergeEdgesUndo(const scalar minCos, const dictionary&);
-
-
         // Layers
 
             //- For debugging: Dump displacement to .obj files
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index e05db743c2..03c29aefd3 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,8 +30,8 @@ License
 #include "cellSet.H"
 #include "syncTools.H"
 #include "refinementParameters.H"
-#include "featureEdgeMesh.H"
 #include "refinementSurfaces.H"
+#include "refinementFeatures.H"
 #include "shellSurfaces.H"
 #include "mapDistributePolyMesh.H"
 #include "unitConversion.H"
@@ -48,62 +48,6 @@ defineTypeNameAndDebug(autoRefineDriver, 0);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Read explicit feature edges
-Foam::label Foam::autoRefineDriver::readFeatureEdges
-(
-    const PtrList<dictionary>& featDicts,
-    PtrList<featureEdgeMesh>& featureMeshes,
-    labelList& featureLevels
-) const
-{
-    Info<< "Reading external feature lines." << endl;
-
-    const fvMesh& mesh = meshRefiner_.mesh();
-
-    featureMeshes.setSize(featDicts.size());
-    featureLevels.setSize(featDicts.size());
-
-    forAll(featDicts, i)
-    {
-        const dictionary& dict = featDicts[i];
-
-        fileName featFileName(dict.lookup("file"));
-
-        featureMeshes.set
-        (
-            i,
-            new featureEdgeMesh
-            (
-                IOobject
-                (
-                    featFileName,                       // name
-                    //mesh.time().findInstance("triSurface", featFileName),
-                    //                                    // instance
-                    mesh.time().constant(),             // instance
-                    "triSurface",                       // local
-                    mesh.time(),                        // registry
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE,
-                    false
-                )
-            )
-        );
-
-        featureMeshes[i].mergePoints(meshRefiner_.mergeDistance());
-        featureLevels[i] = readLabel(dict.lookup("level"));
-
-        Info<< "Refinement level " << featureLevels[i]
-            << " for all cells crossed by feature " << featFileName
-            << " (" << featureMeshes[i].points().size() << " points, "
-            << featureMeshes[i].edges().size() << " edges)." << endl;
-    }
-
-    Info<< "Read feature lines in = "
-        << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
-
-    return featureMeshes.size();
-}
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -128,23 +72,15 @@ Foam::autoRefineDriver::autoRefineDriver
 Foam::label Foam::autoRefineDriver::featureEdgeRefine
 (
     const refinementParameters& refineParams,
-    const PtrList<dictionary>& featDicts,
     const label maxIter,
     const label minRefine
 )
 {
     const fvMesh& mesh = meshRefiner_.mesh();
 
-    // Read explicit feature edges
-    PtrList<featureEdgeMesh> featureMeshes;
-    // Per feature the refinement level
-    labelList featureLevels;
-    readFeatureEdges(featDicts, featureMeshes, featureLevels);
-
-
     label iter = 0;
 
-    if (featureMeshes.size() && maxIter > 0)
+    if (meshRefiner_.features().size() && maxIter > 0)
     {
         for (; iter < maxIter; iter++)
         {
@@ -160,9 +96,6 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
                     refineParams.keepPoints()[0],    // For now only use one.
                     refineParams.curvature(),
 
-                    featureMeshes,
-                    featureLevels,
-
                     true,               // featureRefinement
                     false,              // internalRefinement
                     false,              // surfaceRefinement
@@ -266,8 +199,6 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
         // Only look at surface intersections (minLevel and surface curvature),
         // do not do internal refinement (refinementShells)
 
-        const PtrList<featureEdgeMesh> dummyFeatures;
-
         labelList candidateCells
         (
             meshRefiner_.refineCandidates
@@ -275,9 +206,6 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
                 refineParams.keepPoints()[0],
                 refineParams.curvature(),
 
-                dummyFeatures,      // dummy featureMeshes;
-                labelList(0),       // dummy featureLevels;
-
                 false,              // featureRefinement
                 false,              // internalRefinement
                 true,               // surfaceRefinement
@@ -432,8 +360,6 @@ Foam::label Foam::autoRefineDriver::shellRefine
             << "----------------------------" << nl
             << endl;
 
-        const PtrList<featureEdgeMesh> dummyFeatures;
-
         labelList candidateCells
         (
             meshRefiner_.refineCandidates
@@ -441,9 +367,6 @@ Foam::label Foam::autoRefineDriver::shellRefine
                 refineParams.keepPoints()[0],
                 refineParams.curvature(),
 
-                dummyFeatures,      // dummy featureMeshes;
-                labelList(0),       // dummy featureLevels;
-
                 false,              // featureRefinement
                 true,               // internalRefinement
                 false,              // surfaceRefinement
@@ -459,7 +382,7 @@ Foam::label Foam::autoRefineDriver::shellRefine
                 << " cells to cellSet candidateCellsFromShells." << endl;
 
             cellSet c(mesh, "candidateCellsFromShells", candidateCells);
-            c.instance() = mesh.time().timeName();
+            c.instance() = meshRefiner_.timeName();
             c.write();
         }
 
@@ -677,6 +600,7 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
         handleSnapProblems,                 // remove perp edge connected cells
         perpAngle,                          // perp angle
         false,                              // merge free standing baffles?
+        //true,                               // merge free standing baffles?
         motionDict,
         const_cast<Time&>(mesh.time()),
         globalToPatch_,
@@ -735,7 +659,8 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
 
 void Foam::autoRefineDriver::mergePatchFaces
 (
-    const refinementParameters& refineParams
+    const refinementParameters& refineParams,
+    const dictionary& motionDict
 )
 {
     const fvMesh& mesh = meshRefiner_.mesh();
@@ -750,11 +675,12 @@ void Foam::autoRefineDriver::mergePatchFaces
         const_cast<Time&>(mesh.time())++;
     }
 
-    meshRefiner_.mergePatchFaces
+    meshRefiner_.mergePatchFacesUndo
     (
         Foam::cos(degToRad(45.0)),
         Foam::cos(degToRad(45.0)),
-        meshRefiner_.meshedPatches()
+        meshRefiner_.meshedPatches(),
+        motionDict
     );
 
     if (debug)
@@ -762,7 +688,7 @@ void Foam::autoRefineDriver::mergePatchFaces
         meshRefiner_.checkData();
     }
 
-    meshRefiner_.mergeEdges(Foam::cos(degToRad(45.0)));
+    meshRefiner_.mergeEdgesUndo(Foam::cos(degToRad(45.0)), motionDict);
 
     if (debug)
     {
@@ -789,13 +715,10 @@ void Foam::autoRefineDriver::doRefine
     // Check that all the keep points are inside the mesh.
     refineParams.findCells(mesh);
 
-    PtrList<dictionary> featDicts(refineDict.lookup("features"));
-
     // Refine around feature edges
     featureEdgeRefine
     (
         refineParams,
-        featDicts,
         100,    // maxIter
         0       // min cells to refine
     );
@@ -833,7 +756,7 @@ void Foam::autoRefineDriver::doRefine
     // Do something about cells with refined faces on the boundary
     if (prepareForSnapping)
     {
-        mergePatchFaces(refineParams);
+        mergePatchFaces(refineParams, motionDict);
     }
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
index a47d40082e..abe8a53be6 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,7 +42,6 @@ namespace Foam
 {
 
 // Forward declaration of classes
-class featureEdgeMesh;
 class refinementParameters;
 class meshRefinement;
 class decompositionMethod;
@@ -71,19 +70,10 @@ class autoRefineDriver
 
     // Private Member Functions
 
-        //- Read explicit feature edges
-        label readFeatureEdges
-        (
-            const PtrList<dictionary>& featDicts,
-            PtrList<featureEdgeMesh>& featureMeshes,
-            labelList& featureLevel
-        ) const;
-
         //- Refine all cells pierced by explicit feature edges
         label featureEdgeRefine
         (
             const refinementParameters& refineParams,
-            const PtrList<dictionary>& featDicts,
             const label maxIter,
             const label minRefine
         );
@@ -130,7 +120,8 @@ class autoRefineDriver
         //- Merge refined boundary faces (from exposing coarser cell)
         void mergePatchFaces
         (
-            const refinementParameters& refineParams
+            const refinementParameters& refineParams,
+            const dictionary& motionDict
         );
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index c6455cdc10..3787156a53 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,7 +27,6 @@ Description
 \*----------------------------------------------------------------------------*/
 
 #include "autoSnapDriver.H"
-#include "Time.H"
 #include "motionSmoother.H"
 #include "polyTopoChange.H"
 #include "OFstream.H"
@@ -36,13 +35,12 @@ Description
 #include "Time.H"
 #include "OFstream.H"
 #include "mapPolyMesh.H"
-#include "motionSmoother.H"
 #include "pointEdgePoint.H"
 #include "PointEdgeWave.H"
 #include "mergePoints.H"
 #include "snapParameters.H"
 #include "refinementSurfaces.H"
-
+#include "unitConversion.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -509,9 +507,7 @@ void Foam::autoSnapDriver::dumpMove
 )
 {
     // Dump direction of growth into file
-    Pout<< nl << "Dumping move direction to " << fName << nl
-        << "View this Lightwave-OBJ file with e.g. javaview" << nl
-        << endl;
+    Pout<< nl << "Dumping move direction to " << fName << endl;
 
     OFstream nearestStream(fName);
 
@@ -652,8 +648,7 @@ Foam::scalarField Foam::autoSnapDriver::calcSnapDistance
         -GREAT              // null value
     );
 
-    tmp<scalarField> tfld = snapParams.snapTol()*maxEdgeLen;
-    return tfld();
+    return snapParams.snapTol()*maxEdgeLen;
 }
 
 
@@ -688,6 +683,7 @@ void Foam::autoSnapDriver::preSmoothPatch
 
         // The current mesh is the starting mesh to smooth from.
         meshMover.setDisplacement(patchDisp);
+
         meshMover.correct();
 
         scalar oldErrorReduction = -1;
@@ -828,7 +824,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
             (
                 unzonedSurfaces,
                 localPoints,
-                sqr(4*snapDist),        // sqr of attract distance
+                sqr(snapDist),        // sqr of attract distance
                 hitSurface,
                 hitInfo
             );
@@ -880,7 +876,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
             (
                 labelList(1, zoneSurfI),
                 pointField(localPoints, zonePointIndices),
-                sqr(4*scalarField(minSnapDist, zonePointIndices)),
+                sqr(scalarField(minSnapDist, zonePointIndices)),
                 hitSurface,
                 hitInfo
             );
@@ -956,8 +952,8 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
         mesh,
         pp.meshPoints(),
         patchDisp,
-        minMagEqOp(),                   // combine op
-        vector(GREAT, GREAT, GREAT)     // null value
+        minMagSqrEqOp<point>(),         // combine op
+        vector(GREAT, GREAT, GREAT)     // null value (note: cannot use VGREAT)
     );
 
 
@@ -990,13 +986,12 @@ void Foam::autoSnapDriver::smoothDisplacement
 ) const
 {
     const fvMesh& mesh = meshRefiner_.mesh();
-    const pointMesh& pMesh = meshMover.pMesh();
     const indirectPrimitivePatch& pp = meshMover.patch();
 
     Info<< "Smoothing displacement ..." << endl;
 
     // Set edge diffusivity as inverse of distance to patch
-    scalarField edgeGamma(1.0/(edgePatchDist(pMesh, pp) + SMALL));
+    scalarField edgeGamma(1.0/(edgePatchDist(meshMover.pMesh(), pp) + SMALL));
     //scalarField edgeGamma(mesh.nEdges(), 1.0);
     //scalarField edgeGamma(wallGamma(mesh, pp, 10, 1));
 
@@ -1010,7 +1005,6 @@ void Foam::autoSnapDriver::smoothDisplacement
             Info<< "Iteration " << iter << endl;
         }
         pointVectorField oldDisp(disp);
-
         meshMover.smooth(oldDisp, edgeGamma, disp);
     }
     Info<< "Displacement smoothed in = "
@@ -1045,7 +1039,7 @@ void Foam::autoSnapDriver::smoothDisplacement
 }
 
 
-void Foam::autoSnapDriver::scaleMesh
+bool Foam::autoSnapDriver::scaleMesh
 (
     const snapParameters& snapParams,
     const label nInitErrors,
@@ -1061,6 +1055,8 @@ void Foam::autoSnapDriver::scaleMesh
 
     scalar oldErrorReduction = -1;
 
+    bool meshOk = false;
+
     Info<< "Moving mesh ..." << endl;
     for (label iter = 0; iter < 2*snapParams.nSnap(); iter++)
     {
@@ -1072,10 +1068,11 @@ void Foam::autoSnapDriver::scaleMesh
             oldErrorReduction = meshMover.setErrorReduction(0.0);
         }
 
-        if (meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors))
+        meshOk = meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors);
+
+        if (meshOk)
         {
             Info<< "Successfully moved mesh" << endl;
-
             break;
         }
         if (debug)
@@ -1098,6 +1095,8 @@ void Foam::autoSnapDriver::scaleMesh
     }
     Info<< "Moved mesh in = "
         << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
+
+    return meshOk;
 }
 
 
@@ -1191,7 +1190,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
         (
             unzonedSurfaces,
             localFaceCentres,
-            sqr(4*faceSnapDist),    // sqr of attract distance
+            sqr(faceSnapDist),    // sqr of attract distance
             hitSurface,
             hitRegion
         );
@@ -1260,6 +1259,7 @@ void Foam::autoSnapDriver::doSnap
 (
     const dictionary& snapDict,
     const dictionary& motionDict,
+    const scalar featureCos,
     const snapParameters& snapParams
 )
 {
@@ -1278,6 +1278,21 @@ void Foam::autoSnapDriver::doSnap
     List<labelPair> baffles;
     meshRefiner_.createZoneBaffles(globalToPatch_, baffles);
 
+
+    bool doFeatures = false;
+    label nFeatIter = 1;
+    if (snapParams.nFeatureSnap() > 0)
+    {
+        doFeatures = true;
+        nFeatIter = snapParams.nFeatureSnap();
+
+        Info<< "Snapping to features in " << nFeatIter
+            << " iterations ..." << endl;
+    }
+
+
+    bool meshOk = false;
+
     {
         autoPtr<indirectPrimitivePatch> ppPtr
         (
@@ -1287,10 +1302,9 @@ void Foam::autoSnapDriver::doSnap
                 adaptPatchIDs
             )
         );
-        indirectPrimitivePatch& pp = ppPtr();
 
         // Distance to attract to nearest feature on surface
-        const scalarField snapDist(calcSnapDistance(snapParams, pp));
+        const scalarField snapDist(calcSnapDistance(snapParams, ppPtr()));
 
 
         // Construct iterative mesh mover.
@@ -1302,7 +1316,7 @@ void Foam::autoSnapDriver::doSnap
         motionSmoother meshMover
         (
             mesh,
-            pp,
+            ppPtr(),
             adaptPatchIDs,
             meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs),
             motionDict
@@ -1330,19 +1344,97 @@ void Foam::autoSnapDriver::doSnap
         // Pre-smooth patch vertices (so before determining nearest)
         preSmoothPatch(snapParams, nInitErrors, baffles, meshMover);
 
-        // Calculate displacement at every patch point. Insert into
-        // meshMover.
-        calcNearestSurface(snapDist, meshMover);
+        if (debug)
+        {
+            Pout<< "Writing patch smoothed mesh to time "
+                << meshRefiner_.timeName() << '.' << endl;
+            meshRefiner_.write
+            (
+                debug,
+                mesh.time().path()/meshRefiner_.timeName()
+            );
+            Pout<< "Dumped mesh in = "
+                << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
+        }
 
-        //// Get smoothly varying internal displacement field.
-        //- 2009-12-16 : was not found to be beneficial. Keeping internal
-        // fields fixed slightly increases skewness (on boundary)
-        // but lowers non-orthogonality quite a bit (e.g. 65->59 degrees).
-        // Maybe if better smoother?
-        //smoothDisplacement(snapParams, meshMover);
 
-        // Apply internal displacement to mesh.
-        scaleMesh(snapParams, nInitErrors, baffles, meshMover);
+        for (label iter = 0; iter < nFeatIter; iter++)
+        {
+            Info<< nl
+                << "Morph iteration " << iter << nl
+                << "-----------------" << endl;
+
+            // Calculate displacement at every patch point. Insert into
+            // meshMover.
+            vectorField disp = calcNearestSurface(snapDist, meshMover);
+
+            // Override displacement with feature edge attempt
+            if (doFeatures)
+            {
+                disp = calcNearestSurfaceFeature
+                (
+                    iter,
+                    featureCos,
+                    scalar(iter+1)/nFeatIter,
+                    snapDist,
+                    disp,
+                    meshMover
+                );
+            }
+
+            if (debug&meshRefinement::OBJINTERSECTIONS)
+            {
+                dumpMove
+                (
+                    mesh.time().path()
+                  / "patchDisplacement_" + name(iter) + ".obj",
+                    ppPtr().localPoints(),
+                    ppPtr().localPoints() + disp
+                );
+            }
+
+            // Get smoothly varying internal displacement field.
+            smoothDisplacement(snapParams, meshMover);
+
+            // Apply internal displacement to mesh.
+            meshOk = scaleMesh
+            (
+                snapParams,
+                nInitErrors,
+                baffles,
+                meshMover
+            );
+
+            if (!meshOk)
+            {
+                Info<< "Did not succesfully snap mesh. Giving up."
+                    << nl << endl;
+
+                // Use current mesh as base mesh
+                meshMover.correct();
+
+                break;
+            }
+
+            if (debug)
+            {
+                const_cast<Time&>(mesh.time())++;
+                Pout<< "Writing scaled mesh to time "
+                    << meshRefiner_.timeName() << endl;
+                meshRefiner_.write
+                (
+                    debug,
+                    mesh.time().path()/meshRefiner_.timeName()
+                );
+                Pout<< "Writing displacement field ..." << endl;
+                meshMover.displacement().write();
+                tmp<pointScalarField> magDisp(mag(meshMover.displacement()));
+                magDisp().write();
+            }
+
+            // Use current mesh as base mesh
+            meshMover.correct();
+        }
     }
 
     // Merge any introduced baffles.
@@ -1350,6 +1442,34 @@ void Foam::autoSnapDriver::doSnap
 
     // Repatch faces according to nearest.
     repatchToSurface(snapParams, adaptPatchIDs);
+
+    // Repatching might have caused faces to be on same patch and hence
+    // mergeable so try again to merge coplanar faces
+    label nChanged = meshRefiner_.mergePatchFacesUndo
+    (
+        featureCos,  // minCos
+        featureCos,  // concaveCos
+        meshRefiner_.meshedPatches(),
+        motionDict
+    );
+
+    nChanged += meshRefiner_.mergeEdgesUndo
+    (
+        featureCos,
+        motionDict
+    );
+
+    if (nChanged > 0 && debug)
+    {
+        const_cast<Time&>(mesh.time())++;
+        Pout<< "Writing patchFace merged mesh to time "
+            << meshRefiner_.timeName() << endl;
+        meshRefiner_.write
+        (
+            debug,
+            mesh.time().path()/meshRefiner_.timeName()
+        );
+    }
 }
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index 6264d4c8cc..4b0e1a30fd 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,6 +29,7 @@ Description
 
 SourceFiles
     autoSnapDriver.C
+    autoSnapDriverFeature.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -45,6 +46,7 @@ namespace Foam
 // Forward declaration of classes
 class motionSmoother;
 class snapParameters;
+class pointConstraint;
 
 /*---------------------------------------------------------------------------*\
                            Class autoSnapDriver Declaration
@@ -52,23 +54,6 @@ class snapParameters;
 
 class autoSnapDriver
 {
-    // Private classes
-
-        //- Combine operator class for equalizing displacements.
-        class minMagEqOp
-        {
-        public:
-
-            void operator()(vector& x, const vector& y) const
-            {
-                if (magSqr(y) < magSqr(x))
-                {
-                    x = y;
-                }
-            }
-        };
-
-
     // Private data
 
         //- Mesh+surface
@@ -124,6 +109,174 @@ class autoSnapDriver
                 const vectorField&
             );
 
+            // Feature line snapping
+
+                void smoothAndConstrain
+                (
+                    const indirectPrimitivePatch& pp,
+                    const List<pointConstraint>& constraints,
+                    vectorField& disp
+                ) const;
+                void calcNearest
+                (
+                    const pointField& points,
+                    vectorField& disp,
+                    vectorField& surfaceNormal
+                ) const;
+                void calcNearestFace
+                (
+                    const label iter,
+                    const indirectPrimitivePatch& pp,
+                    vectorField& faceDisp,
+                    vectorField& faceSurfaceNormal,
+                    vectorField& faceRotation
+                ) const;
+                void interpolateFaceToPoint
+                (
+                    const label iter,
+                    const indirectPrimitivePatch& pp,
+                    const vectorField& faceSurfaceNormal,
+
+                    const vectorField& faceDisp,
+                    const vectorField& faceRotation,
+
+                    vectorField& patchDisp,
+                    vectorField& patchRotationDisp
+                ) const;
+                void correctAttraction
+                (
+                    const DynamicList<point>& surfacePoints,
+                    const DynamicList<label>& surfaceCount,
+                    const point& edgePt,
+                    const vector& edgeNormal,   // normalised normal
+                    const point& pt,
+                    vector& edgeOffset  // offset from pt to point on edge
+                ) const;
+                void binFeatureFace
+                (
+                    const label iter,
+                    const scalar featureCos,
+
+                    const indirectPrimitivePatch& pp,
+                    const scalar snapDist,
+
+                    const point& fc,
+                    const vector& faceSurfaceNormal,
+                    const vector& faceDisp,
+
+                    DynamicList<point>& surfacePoints,
+                    DynamicList<vector>& surfaceNormals,
+                    DynamicList<label>& surfaceCount
+                ) const;
+                void binFeatureFaces
+                (
+                    const label iter,
+                    const scalar featureCos,
+
+                    const indirectPrimitivePatch& pp,
+                    const scalarField& snapDist,
+                    const label pointI,
+
+                    const List<List<point> >& pointFaceNormals,
+                    const List<List<point> >& pointFaceDisp,
+                    const List<List<point> >& pointFaceCentres,
+
+                    DynamicList<point>& surfacePoints,
+                    DynamicList<vector>& surfaceNormals,
+                    DynamicList<label>& surfaceCount
+                ) const;
+
+                void featureAttractionUsingReconstruction
+                (
+                    const label iter,
+                    const scalar featureCos,
+                    const indirectPrimitivePatch& pp,
+                    const scalarField& snapDist,
+
+                    const List<List<point> >& pointFaceNormals,
+                    const List<List<point> >& pointFaceDisp,
+                    const List<List<point> >& pointFaceCentres,
+
+                    vectorField& patchAttraction,
+                    List<pointConstraint>& patchConstraints
+                ) const;
+
+                void determineAllFeatures
+                (
+                    const label iter,
+                    const scalar featureCos,
+
+                    const indirectPrimitivePatch&,
+                    const scalarField&,
+
+                    const List<List<point> >& pointFaceNormals,
+                    const List<List<point> >& pointFaceDisp,
+                    const List<List<point> >& pointFaceCentres,
+
+                    List<labelList>& pointAttractor,
+                    List<List<pointConstraint> >& pointConstraints,
+                    // Feature-edge to pp point
+                    List<List<DynamicList<point> > >& edgeAttractors,
+                    List<List<DynamicList<pointConstraint> > >& edgeConstraints,
+                    vectorField& patchAttraction,
+                    List<pointConstraint>& patchConstraints
+                ) const;
+                void determineFeatures
+                (
+                    const label iter,
+                    const scalar featureCos,
+
+                    const indirectPrimitivePatch&,
+                    const scalarField&,
+
+                    const List<List<point> >& pointFaceNormals,
+                    const List<List<point> >& pointFaceDisp,
+                    const List<List<point> >& pointFaceCentres,
+
+                    List<labelList>& pointAttractor,
+                    List<List<pointConstraint> >& pointConstraints,
+                    // Feature-edge to pp point
+                    List<List<DynamicList<point> > >& edgeAttractors,
+                    List<List<DynamicList<pointConstraint> > >& edgeConstraints,
+                    vectorField& patchAttraction,
+                    List<pointConstraint>& patchConstraints
+                ) const;
+
+                void featureAttractionUsingFeatureEdges
+                (
+                    const label iter,
+                    const scalar featureCos,
+                    const indirectPrimitivePatch& pp,
+                    const scalarField& snapDist,
+
+                    const List<List<point> >& pointFaceNormals,
+                    const List<List<point> >& pointFaceDisp,
+                    const List<List<point> >& pointFaceCentres,
+
+                    vectorField& patchAttraction,
+                    List<pointConstraint>& patchConstraints
+                ) const;
+                void preventFaceSqueeze
+                (
+                    const label iter,
+                    const scalar featureCos,
+                    const indirectPrimitivePatch& pp,
+                    const scalarField& snapDist,
+
+                    vectorField& patchAttraction,
+                    List<pointConstraint>& patchConstraints
+                ) const;
+
+                vectorField calcNearestSurfaceFeature
+                (
+                    const label iter,
+                    const scalar featureCos,
+                    const scalar featureAttract,
+                    const scalarField& snapDist,
+                    const vectorField& nearestDisp,
+                    motionSmoother& meshMover
+                ) const;
+
 
         //- Disallow default bitwise copy construct
         autoSnapDriver(const autoSnapDriver&);
@@ -188,6 +341,7 @@ public:
                 motionSmoother& meshMover
             ) const;
 
+
             //- Smooth the displacement field to the internal.
             void smoothDisplacement
             (
@@ -196,8 +350,9 @@ public:
             ) const;
 
             //- Do the hard work: move the mesh according to displacement,
-            //  locally relax the displacement.
-            void scaleMesh
+            //  locally relax the displacement. Return true if ended up with
+            //  correct mesh, false if not.
+            bool scaleMesh
             (
                 const snapParameters& snapParams,
                 const label nInitErrors,
@@ -212,11 +367,11 @@ public:
                 const labelList& adaptPatchIDs
             );
 
-
             void doSnap
             (
                 const dictionary& snapDict,
                 const dictionary& motionDict,
+                const scalar featureCos,
                 const snapParameters& snapParams
             );
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
new file mode 100644
index 0000000000..cc73efbe51
--- /dev/null
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -0,0 +1,2146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "autoSnapDriver.H"
+#include "polyTopoChange.H"
+#include "OFstream.H"
+#include "syncTools.H"
+#include "fvMesh.H"
+#include "OFstream.H"
+#include "motionSmoother.H"
+#include "refinementSurfaces.H"
+#include "refinementFeatures.H"
+#include "unitConversion.H"
+#include "plane.H"
+#include "featureEdgeMesh.H"
+#include "treeDataPoint.H"
+#include "indexedOctree.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    class listTransform
+    {
+    public:
+
+        void operator()
+        (
+            const vectorTensorTransform& vt,
+            const bool forward,
+            List<List<point> >& fld
+        ) const
+        {
+            const tensor T
+            (
+                forward
+              ? vt.R()
+              : vt.R().T()
+            );
+
+            forAll(fld, i)
+            {
+                List<point>& elems = fld[i];
+                forAll(elems, elemI)
+                {
+                    elems[elemI] = transform(T, elems[elemI]);
+                }
+            }
+        }
+    };
+
+    class listPlusEqOp
+    {
+    public:
+
+        void operator()
+        (
+            List<point>& x,
+            const List<point>& y
+        ) const
+        {
+            label sz = x.size();
+            x.setSize(sz+y.size());
+            forAll(y, i)
+            {
+                x[sz++] = y[i];
+            }
+        }
+    };
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::autoSnapDriver::smoothAndConstrain
+(
+    const indirectPrimitivePatch& pp,
+    const List<pointConstraint>& constraints,
+    vectorField& disp
+) const
+{
+    for (label avgIter = 0; avgIter < 20; avgIter++)
+    {
+        // Calculate average displacement of neighbours
+        vectorField dispAvg(pp.nPoints(), vector::zero);
+
+        const labelListList& pointEdges = pp.pointEdges();
+
+        forAll(pointEdges, pointI)
+        {
+            const labelList& pEdges = pointEdges[pointI];
+
+            forAll(pEdges, i)
+            {
+                const edge& e = pp.edges()[pEdges[i]];
+                label nbrPointI = e.otherVertex(pointI);
+                dispAvg[pointI] += disp[nbrPointI];
+            }
+            dispAvg[pointI] /= pEdges.size();
+        }
+
+        // Constraints
+        forAll(constraints, pointI)
+        {
+            if (constraints[pointI].first() == 0)
+            {
+                // Mix my displacement with neighbours' displacement
+                disp[pointI] = 0.5*disp[pointI] + 0.5*dispAvg[pointI];
+            }
+        }
+    }
+}
+
+
+void Foam::autoSnapDriver::calcNearest
+(
+    const pointField& points,
+    vectorField& disp,
+    vectorField& surfaceNormal
+) const
+{
+    const refinementSurfaces& surfaces = meshRefiner_.surfaces();
+
+    // Displacement and orientation per pp face.
+    disp.setSize(points.size());
+    disp = vector::zero;
+    surfaceNormal.setSize(points.size());
+    surfaceNormal = vector::zero;
+
+    {
+        // Divide surfaces into zoned and unzoned
+        labelList zonedSurfaces =
+            meshRefiner_.surfaces().getNamedSurfaces();
+        labelList unzonedSurfaces =
+            meshRefiner_.surfaces().getUnnamedSurfaces();
+
+        {
+            List<pointIndexHit> hitInfo;
+            labelList hitSurface;
+            labelList hitRegion;
+            surfaces.findNearestRegion
+            (
+                unzonedSurfaces,
+                points,
+                sqr(scalarField(points.size(), GREAT)),// sqr of attract dist
+                hitSurface,
+                hitInfo,
+                hitRegion,
+                surfaceNormal
+            );
+
+            forAll(hitInfo, i)
+            {
+                if (hitInfo[i].hit())
+                {
+                    disp[i] =
+                        hitInfo[i].hitPoint()
+                      - points[i];
+                }
+                else
+                {
+                    WarningIn("calcNearest(..)")
+                        << "Did not hit anything from face:" << i
+                        << " at:" << points[i] << endl;
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::autoSnapDriver::calcNearestFace
+(
+    const label iter,
+    const indirectPrimitivePatch& pp,
+    vectorField& faceDisp,
+    vectorField& faceSurfaceNormal,
+    vectorField& faceRotation
+) const
+{
+    const fvMesh& mesh = meshRefiner_.mesh();
+
+    // Displacement and orientation per pp face.
+    faceDisp.setSize(pp.size());
+    faceDisp = vector::zero;
+    faceSurfaceNormal.setSize(pp.size());
+    faceSurfaceNormal = vector::zero;
+
+    calcNearest(pp.faceCentres(), faceDisp, faceSurfaceNormal);
+
+    // Determine rotation axis
+    faceRotation.setSize(pp.size());
+    faceRotation = vector::zero;
+
+    forAll(faceRotation, faceI)
+    {
+        // Note: extend to >180 degrees checking
+        faceRotation[faceI] =
+            pp.faceNormals()[faceI]
+          ^ -faceSurfaceNormal[faceI];
+    }
+
+    if (debug&meshRefinement::OBJINTERSECTIONS)
+    {
+        dumpMove
+        (
+            mesh.time().path()
+          / "faceDisp_" + name(iter) + ".obj",
+            pp.faceCentres(),
+            pp.faceCentres() + faceDisp
+        );
+        dumpMove
+        (
+            mesh.time().path()
+          / "faceRotation_" + name(iter) + ".obj",
+            pp.faceCentres(),
+            pp.faceCentres() + faceRotation
+        );
+    }
+}
+
+
+void Foam::autoSnapDriver::interpolateFaceToPoint
+(
+    const label iter,
+    const indirectPrimitivePatch& pp,
+    const vectorField& faceSurfaceNormal,
+    const vectorField& faceDisp,
+    const vectorField& faceRotation,
+
+    vectorField& patchDisp,
+    vectorField& patchRotationDisp
+) const
+{
+    const fvMesh& mesh = meshRefiner_.mesh();
+
+    // Displacement due to normal snapping
+    patchDisp.setSize(pp.nPoints());
+    patchDisp = vector::zero;
+    // Displacement due to rotation
+    patchRotationDisp.setSize(pp.nPoints());
+    patchRotationDisp = vector::zero;
+
+
+    // Unweighted interpolation
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelList nPatchDisp(pp.nPoints(), 0.0);
+
+    forAll(pp, faceI)
+    {
+        const face& f = pp.localFaces()[faceI];
+        const point& fc = pp.faceCentres()[faceI];
+
+        forAll(f, fp)
+        {
+            label pointI = f[fp];
+            const point& pt = pp.localPoints()[pointI];
+
+            patchDisp[pointI] += faceDisp[faceI];
+            patchRotationDisp[pointI] += (faceRotation[faceI]^(pt-fc));
+            nPatchDisp[pointI]++;
+        }
+    }
+
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        patchDisp,
+        plusEqOp<vector>(),
+        vector::zero,
+        mapDistribute::transform()
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        patchRotationDisp,
+        plusEqOp<vector>(),
+        vector::zero,
+        mapDistribute::transform()
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        nPatchDisp,
+        plusEqOp<label>(),
+        0,
+        mapDistribute::transform()
+    );
+
+    forAll(patchDisp, pointI)
+    {
+        patchDisp[pointI] /= nPatchDisp[pointI];
+        patchRotationDisp[pointI] /= nPatchDisp[pointI];
+    }
+
+
+    if (debug&meshRefinement::OBJINTERSECTIONS)
+    {
+        dumpMove
+        (
+            mesh.time().path()
+          / "patchDisp_" + name(iter) + ".obj",
+            pp.localPoints(),
+            pp.localPoints() + patchDisp
+        );
+        dumpMove
+        (
+            mesh.time().path()
+          / "patchRotationDisp_" + name(iter) + ".obj",
+            pp.localPoints(),
+            pp.localPoints() + patchRotationDisp
+        );
+    }
+}
+
+
+// Gets passed in offset to nearest point on feature edge. Calculates
+// if the point has a different number of faces on either side of the feature
+// and if so attracts the point to that non-dominant plane.
+void Foam::autoSnapDriver::correctAttraction
+(
+    const DynamicList<point>& surfacePoints,
+    const DynamicList<label>& surfaceCount,
+    const point& edgePt,
+    const vector& edgeNormal,       // normalised normal
+    const point& pt,
+
+    vector& edgeOffset              // offset from pt to point on edge
+) const
+{
+    // Tangential component along edge
+    scalar tang = ((pt-edgePt)&edgeNormal);
+
+    labelList order;
+    Foam::sortedOrder(surfaceCount, order);
+
+    if (order[0] < order[1])
+    {
+        // There is a non-dominant plane. Use the point on the plane to
+        // attract to.
+        vector attractD = surfacePoints[order[0]]-edgePt;
+        // Tangential component along edge
+        scalar tang2 = (attractD&edgeNormal);
+        // Normal component
+        attractD -= tang2*edgeNormal;
+        // Calculate fraction of normal distances
+        scalar magAttractD = mag(attractD);
+        scalar fraction = magAttractD/(magAttractD+mag(edgeOffset));
+
+        point linePt =
+            edgePt
+          + ((1.0-fraction)*tang2 + fraction*tang)*edgeNormal;
+        edgeOffset = linePt-pt;
+    }
+}
+
+
+void Foam::autoSnapDriver::binFeatureFace
+(
+    const label iter,
+    const scalar featureCos,
+
+    const indirectPrimitivePatch& pp,
+    const scalar snapDist,
+
+    const point& fc,
+    const vector& faceSurfaceNormal,
+    const vector& faceDisp,
+
+    DynamicList<point>& surfacePoints,
+    DynamicList<vector>& surfaceNormals,
+    DynamicList<label>& surfaceCount
+) const
+{
+    // What to do with very far attraction? For now just ignore the face
+    if (magSqr(faceDisp) < sqr(snapDist))
+    {
+        const point pt = fc + faceDisp;
+
+        bool same = false;
+        forAll(surfaceNormals, j)
+        {
+            scalar cosAngle = (faceSurfaceNormal&surfaceNormals[j]);
+
+            if
+            (
+                (cosAngle >= featureCos)
+             || (cosAngle < (-1+0.001)) // triangle baffles
+            )
+            {
+                same = true;
+                surfaceCount[j]++;
+                break;
+            }
+        }
+
+        if (!same)
+        {
+            // Now check if the planes go through the same edge or point
+
+            if (surfacePoints.size() <= 1)
+            {
+                surfacePoints.append(pt);
+                surfaceNormals.append(faceSurfaceNormal);
+                surfaceCount.append(1);
+            }
+            else if (surfacePoints.size() == 2)
+            {
+                plane pl0(surfacePoints[0], surfaceNormals[0]);
+                plane pl1(surfacePoints[1], surfaceNormals[1]);
+                plane::ray r(pl0.planeIntersect(pl1));
+                vector featureNormal = r.dir() / mag(r.dir());
+
+                if (mag(faceSurfaceNormal&featureNormal) >= 0.001)
+                {
+                    // Definitely makes a feature point
+                    surfacePoints.append(pt);
+                    surfaceNormals.append(faceSurfaceNormal);
+                    surfaceCount.append(1);
+                }
+            }
+            else if (surfacePoints.size() == 3)
+            {
+                // Have already feature point. See if this new plane is the
+                // same point or not.
+                plane pl0(surfacePoints[0], surfaceNormals[0]);
+                plane pl1(surfacePoints[1], surfaceNormals[1]);
+                plane pl2(surfacePoints[2], surfaceNormals[2]);
+                point p012(pl0.planePlaneIntersect(pl1, pl2));
+
+                plane::ray r(pl0.planeIntersect(pl1));
+                vector featureNormal = r.dir() / mag(r.dir());
+                if (mag(faceSurfaceNormal&featureNormal) >= 0.001)
+                {
+                    plane pl3(pt, faceSurfaceNormal);
+                    point p013(pl0.planePlaneIntersect(pl1, pl3));
+
+                    if (mag(p012-p013) > 0.0001)    //TBD
+                    {
+                        // Different feature point
+                        //Pout<< "** differing feature point :" << p012
+                        //    << " and " << p013 << endl;
+                        // Mark point as illegal (for now by setting size
+                        // to 4)
+                        surfacePoints.append(pt);
+                        surfaceNormals.append(faceSurfaceNormal);
+                        surfaceCount.append(1);
+                    }
+                }
+            }
+        }
+    }
+    //else
+    //{
+    //    Pout<< "binFeatureFace : for point:" << pp.localPoints()[pointI]
+    //        << " found face-nearest point:"
+    //        << fc + faceDisp
+    //        << " at distance:" << mag(faceDisp)
+    //        << " further away than:" << snapDist
+    //        << endl;
+    //}
+}
+
+
+// Check the faces surrounding a point. Bin them according to normal.
+void Foam::autoSnapDriver::binFeatureFaces
+(
+    const label iter,
+    const scalar featureCos,
+
+    const indirectPrimitivePatch& pp,
+    const scalarField& snapDist,
+
+    const label pointI,
+
+//    const vectorField& faceSurfaceNormal,
+//    const vectorField& faceDisp,
+    const List<List<point> >& pointFaceNormals,
+    const List<List<point> >& pointFaceDisp,
+    const List<List<point> >& pointFaceCentres,
+
+    DynamicList<point>& surfacePoints,
+    DynamicList<vector>& surfaceNormals,
+    DynamicList<label>& surfaceCount
+) const
+{
+    const List<point>& pfNormals = pointFaceNormals[pointI];
+    const List<point>& pfDisp = pointFaceDisp[pointI];
+    const List<point>& pfCentres = pointFaceCentres[pointI];
+
+    // Collect all different directions
+    forAll(pfNormals, i)
+    {
+        binFeatureFace
+        (
+            iter,
+            featureCos,
+
+            pp,
+            snapDist[pointI],
+
+            pfCentres[i],
+            pfNormals[i],
+            pfDisp[i],
+
+            surfacePoints,
+            surfaceNormals,
+            surfaceCount
+        );
+    }
+}
+
+
+// Special version that calculates attraction in one go
+void Foam::autoSnapDriver::featureAttractionUsingReconstruction
+(
+    const label iter,
+    const scalar featureCos,
+
+    const indirectPrimitivePatch& pp,
+    const scalarField& snapDist,
+
+    //const vectorField& faceSurfaceNormal,
+    //const vectorField& faceDisp,
+    //const vectorField& faceRotation,
+    const List<List<point> >& pointFaceNormals,
+    const List<List<point> >& pointFaceDisp,
+    const List<List<point> >& pointFaceCentres,
+
+    vectorField& patchAttraction,
+    List<pointConstraint>& patchConstraints
+) const
+{
+    autoPtr<OFstream> featureEdgeStr;
+    label featureEdgeVertI = 0;
+    autoPtr<OFstream> featurePointStr;
+    label featurePointVertI = 0;
+
+    if (debug&meshRefinement::OBJINTERSECTIONS)
+    {
+        featureEdgeStr.reset
+        (
+            new OFstream
+            (
+                meshRefiner_.mesh().time().path()
+              / "featureEdge_" + name(iter) + ".obj"
+            )
+        );
+        Pout<< nl << "Dumping feature-edge direction to "
+            << featureEdgeStr().name() << endl;
+
+        featurePointStr.reset
+        (
+            new OFstream
+            (
+                meshRefiner_.mesh().time().path()
+              / "featurePoint_" + name(iter) + ".obj"
+            )
+        );
+        Pout<< nl << "Dumping feature-point direction to "
+            << featurePointStr().name() << endl;
+    }
+
+
+    // For points on multiple normals calculate intersection
+    patchAttraction.setSize(pp.nPoints());
+    patchAttraction = vector::zero;
+
+    forAll(pp.localPoints(), pointI)
+    {
+        // Collect all different directions
+        DynamicList<point> surfacePoints;
+        DynamicList<vector> surfaceNormals;
+        DynamicList<label> surfaceCount;
+
+        binFeatureFaces
+        (
+            iter,
+            featureCos,
+
+            pp,
+            snapDist,
+            pointI,
+
+            //faceSurfaceNormal,
+            //faceDisp,
+            pointFaceNormals,
+            pointFaceDisp,
+            pointFaceCentres,
+
+            surfacePoints,
+            surfaceNormals,
+            surfaceCount
+        );
+
+        const point& pt = pp.localPoints()[pointI];
+
+        // Check the number of directions
+        if (surfaceNormals.size() == 2)
+        {
+            plane pl0(surfacePoints[0], surfaceNormals[0]);
+            plane pl1(surfacePoints[1], surfaceNormals[1]);
+            plane::ray r(pl0.planeIntersect(pl1));
+            vector n = r.dir() / mag(r.dir());
+
+            // Get nearest point on infinite ray
+            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
+            );
+
+            // Trim to snap distance
+            if (magSqr(d) > sqr(snapDist[pointI]))
+            {
+                d *= Foam::sqrt(sqr(snapDist[pointI])/magSqr(d));
+            }
+
+            patchAttraction[pointI] = d;
+
+            // Store constraints
+            patchConstraints[pointI].applyConstraint(surfaceNormals[0]);
+            patchConstraints[pointI].applyConstraint(surfaceNormals[1]);
+
+
+            // Dump vector from point to points on faces
+            if (featureEdgeStr.valid())
+            {
+                meshTools::writeOBJ(featureEdgeStr(), pt);
+                featureEdgeVertI++;
+                meshTools::writeOBJ(featureEdgeStr(), surfacePoints[0]);
+                featureEdgeVertI++;
+                meshTools::writeOBJ(featureEdgeStr(), surfacePoints[1]);
+                featureEdgeVertI++;
+                featureEdgeStr()
+                    << "l " << featureEdgeVertI-2 << ' '
+                    << featureEdgeVertI-1 << nl
+                    << "l " << featureEdgeVertI-2 << ' '
+                    << featureEdgeVertI << nl;
+            }
+        }
+        else if (surfaceNormals.size() == 3)
+        {
+            // Calculate point from the faces.
+            plane pl0(surfacePoints[0], surfaceNormals[0]);
+            plane pl1(surfacePoints[1], surfaceNormals[1]);
+            plane pl2(surfacePoints[2], surfaceNormals[2]);
+            point cornerPt(pl0.planePlaneIntersect(pl1, pl2));
+
+            vector d = cornerPt - pt;
+
+            if (magSqr(d) > sqr(snapDist[pointI]))
+            {
+                d *= Foam::sqrt(sqr(snapDist[pointI])/magSqr(d));
+            }
+
+            patchAttraction[pointI] = d;
+
+            // Store constraints
+            patchConstraints[pointI].applyConstraint(surfaceNormals[0]);
+            patchConstraints[pointI].applyConstraint(surfaceNormals[1]);
+            patchConstraints[pointI].applyConstraint(surfaceNormals[2]);
+
+            // Dump vector from point to points on faces
+            if (featurePointStr.valid())
+            {
+                meshTools::writeOBJ(featurePointStr(), pt);
+                featurePointVertI++;
+                meshTools::writeOBJ(featurePointStr(), surfacePoints[0]);
+                featurePointVertI++;
+                meshTools::writeOBJ(featurePointStr(), surfacePoints[1]);
+                featurePointVertI++;
+                meshTools::writeOBJ(featurePointStr(), surfacePoints[2]);
+                featurePointVertI++;
+                featurePointStr()
+                    << "l " << featurePointVertI-3 << ' '
+                    << featurePointVertI-2 << nl
+                    << "l " << featurePointVertI-3 << ' '
+                    << featurePointVertI-1 << nl
+                    << "l " << featurePointVertI-3 << ' '
+                    << featurePointVertI << nl;
+            }
+        }
+    }
+
+    if (debug&meshRefinement::OBJINTERSECTIONS)
+    {
+        dumpMove
+        (
+            meshRefiner_.mesh().time().path()
+          / "patchAttraction_" + name(iter) + ".obj",
+            pp.localPoints(),
+            pp.localPoints() + patchAttraction
+        );
+    }
+}
+
+
+// Finds nearest feature (within snapDist) for all points of pp.
+void Foam::autoSnapDriver::determineAllFeatures
+(
+    const label iter,
+    const scalar featureCos,
+
+    const indirectPrimitivePatch& pp,
+    const scalarField& snapDist,
+
+    //const vectorField& faceSurfaceNormal,
+    //const vectorField& faceDisp,
+    //const vectorField& faceRotation,
+    const List<List<point> >& pointFaceNormals,
+    const List<List<point> >& pointFaceDisp,
+    const List<List<point> >& pointFaceCentres,
+
+    List<labelList>& pointAttractor,
+    List<List<pointConstraint> >& pointConstraints,
+    // Feature-edge to pp point
+    List<List<DynamicList<point> > >& edgeAttractors,
+    List<List<DynamicList<pointConstraint> > >& edgeConstraints,
+    // pp point to nearest feature
+    vectorField& patchAttraction,
+    List<pointConstraint>& patchConstraints
+) const
+{
+    autoPtr<OFstream> featureEdgeStr;
+    label featureEdgeVertI = 0;
+    autoPtr<OFstream> featurePointStr;
+    label featurePointVertI = 0;
+
+    if (debug&meshRefinement::OBJINTERSECTIONS)
+    {
+        featureEdgeStr.reset
+        (
+            new OFstream
+            (
+                meshRefiner_.mesh().time().path()
+              / "featureEdge_" + name(iter) + ".obj"
+            )
+        );
+        Pout<< nl << "Dumping feature-edge direction to "
+            << featureEdgeStr().name() << endl;
+
+        featurePointStr.reset
+        (
+            new OFstream
+            (
+                meshRefiner_.mesh().time().path()
+              / "featurePoint_" + name(iter) + ".obj"
+            )
+        );
+        Pout<< nl << "Dumping feature-point direction to "
+            << featurePointStr().name() << endl;
+    }
+
+    const refinementFeatures& features = meshRefiner_.features();
+
+    // Look at near feature edges
+    labelList nearEdgeFeat;
+    List<pointIndexHit> nearEdgeInfo;
+
+    labelList nearPointFeat;
+    labelList nearPointIndex;
+    {
+        scalarField snapDistSqr = sqr(snapDist);
+        features.findNearestEdge
+        (
+            pp.localPoints(),
+            snapDistSqr,
+            nearEdgeFeat,
+            nearEdgeInfo
+        );
+
+        // Look at near feature points
+        features.findNearestPoint
+        (
+            pp.localPoints(),
+            snapDistSqr,
+            nearPointFeat,
+            nearPointIndex
+        );
+    }
+
+    forAll(pp.localPoints(), pointI)
+    {
+        const point& pt = pp.localPoints()[pointI];
+        const label featI = nearEdgeFeat[pointI];
+        const pointIndexHit& nearEdge = nearEdgeInfo[pointI];
+
+        // Mark point on the nearest feature edge.
+        if (nearEdge.hit())
+        {
+            // So we have a point on the feature edge. Use this instead
+            // of our estimate from planes.
+
+            label featEdgeI = nearEdge.index();
+            edgeAttractors[featI][featEdgeI].append(nearEdge.hitPoint());
+            const featureEdgeMesh& eMesh = features[nearEdgeFeat[pointI]];
+            const edge& e = eMesh.edges()[featEdgeI];
+            vector eVec = e.vec(eMesh.points());
+            eVec /= mag(eVec);
+            pointConstraint c;
+            c.first() = 2;
+            c.second() = eVec;
+            edgeConstraints[featI][featEdgeI].append(c);
+
+            // Store for later use
+            patchAttraction[pointI] = nearEdge.hitPoint()-pt;
+            patchConstraints[pointI] = c;
+
+            // Dump
+            if (featureEdgeStr.valid())
+            {
+                meshTools::writeOBJ(featureEdgeStr(), pt);
+                featureEdgeVertI++;
+                meshTools::writeOBJ(featureEdgeStr(), nearEdge.hitPoint());
+                featureEdgeVertI++;
+                featureEdgeStr()
+                    << "l " << featureEdgeVertI-1 << ' '
+                    << featureEdgeVertI << nl;
+            }
+        }
+
+
+        // Mark point on the nearest feature point.
+        if (nearPointFeat[pointI] != -1)
+        {
+            label featI = nearPointFeat[pointI];
+            label index = nearPointIndex[pointI];
+            const treeDataPoint& shapes = features.pointTrees()[featI].shapes();
+            label featPointI = shapes.pointLabels()[index];
+            const point& featPt = shapes.points()[featPointI];
+
+            pointAttractor[featI][featPointI] = pointI;
+            pointConstraints[featI][featPointI].first() = 3;
+            pointConstraints[featI][featPointI].second() = vector::zero;
+
+            // Dump
+            if (featurePointStr.valid())
+            {
+                meshTools::writeOBJ(featurePointStr(), pt);
+                featurePointVertI++;
+                meshTools::writeOBJ(featurePointStr(), featPt);
+                featurePointVertI++;
+                featurePointStr()
+                    << "l " << featurePointVertI-1 << ' '
+                    << featurePointVertI << nl;
+            }
+        }
+    }
+}
+
+
+void Foam::autoSnapDriver::determineFeatures
+(
+    const label iter,
+    const scalar featureCos,
+
+    const indirectPrimitivePatch& pp,
+    const scalarField& snapDist,
+
+    //const vectorField& faceSurfaceNormal,
+    //const vectorField& faceDisp,
+    //const vectorField& faceRotation,
+    const List<List<point> >& pointFaceNormals,
+    const List<List<point> >& pointFaceDisp,
+    const List<List<point> >& pointFaceCentres,
+
+    // Feature-point to pp point
+    List<labelList>& pointAttractor,
+    List<List<pointConstraint> >& pointConstraints,
+    // Feature-edge to pp point
+    List<List<DynamicList<point> > >& edgeAttractors,
+    List<List<DynamicList<pointConstraint> > >& edgeConstraints,
+    // pp point to nearest feature
+    vectorField& patchAttraction,
+    List<pointConstraint>& patchConstraints
+) const
+{
+    autoPtr<OFstream> featureEdgeStr;
+    label featureEdgeVertI = 0;
+    autoPtr<OFstream> missedEdgeStr;
+    label missedVertI = 0;
+    autoPtr<OFstream> featurePointStr;
+    label featurePointVertI = 0;
+
+    if (debug&meshRefinement::OBJINTERSECTIONS)
+    {
+        featureEdgeStr.reset
+        (
+            new OFstream
+            (
+                meshRefiner_.mesh().time().path()
+              / "featureEdge_" + name(iter) + ".obj"
+            )
+        );
+        Pout<< nl << "Dumping feature-edge sampling to "
+            << featureEdgeStr().name() << endl;
+
+        missedEdgeStr.reset
+        (
+            new OFstream
+            (
+                meshRefiner_.mesh().time().path()
+              / "missedFeatureEdge_" + name(iter) + ".obj"
+            )
+        );
+        Pout<< nl << "Dumping feature-edges that are too far away to "
+            << missedEdgeStr().name() << endl;
+
+        featurePointStr.reset
+        (
+            new OFstream
+            (
+                meshRefiner_.mesh().time().path()
+              / "featurePoint_" + name(iter) + ".obj"
+            )
+        );
+        Pout<< nl << "Dumping feature-point sampling to "
+            << featurePointStr().name() << endl;
+    }
+
+    const refinementFeatures& features = meshRefiner_.features();
+
+    forAll(pp.localPoints(), pointI)
+    {
+        if (patchConstraints[pointI].first() == 0)
+        {
+            const point& pt = pp.localPoints()[pointI];
+
+            // Collect all different directions
+            DynamicList<point> surfacePoints;
+            DynamicList<vector> surfaceNormals;
+            DynamicList<label> surfaceCount;
+
+            binFeatureFaces
+            (
+                iter,
+                featureCos,
+
+                pp,
+                snapDist,
+                pointI,
+
+                //faceSurfaceNormal,
+                //faceDisp,
+                pointFaceNormals,
+                pointFaceDisp,
+                pointFaceCentres,
+
+                surfacePoints,
+                surfaceNormals,
+                surfaceCount
+            );
+
+            // Check the number of directions
+            if (surfaceNormals.size() == 2)
+            {
+                plane pl0(surfacePoints[0], surfaceNormals[0]);
+                plane pl1(surfacePoints[1], surfaceNormals[1]);
+                plane::ray r(pl0.planeIntersect(pl1));
+                const vector n = r.dir() / mag(r.dir());
+
+                // Get nearest point on infinite ray
+                vector d = r.refPoint()-pt;
+                d -= (d&n)*n;
+
+                // Trim to snap distance
+                if (magSqr(d) > sqr(snapDist[pointI]))
+                {
+                    d *= Foam::sqrt(sqr(snapDist[pointI])/magSqr(d));
+                }
+
+
+                // Mark point on the nearest feature edge. Note that we
+                // only search within the surrounding since the plane
+                // reconstruction might find a feature where there isn't one.
+                const point estimatedPt(pt + d);
+
+                labelList nearEdgeFeat;
+                List<pointIndexHit> nearEdgeInfo;
+                features.findNearestEdge
+                (
+                    pointField(1, estimatedPt),
+                    scalarField(1, sqr(snapDist[pointI])),
+                    nearEdgeFeat,
+                    nearEdgeInfo
+                );
+
+                const pointIndexHit& nearInfo = nearEdgeInfo[0];
+                label featI = nearEdgeFeat[0];
+
+                if (nearInfo.hit())
+                {
+                    // So we have a point on the feature edge. Use this instead
+                    // of our estimate from planes.
+
+                    edgeAttractors[featI][nearInfo.index()].append
+                    (
+                        nearInfo.hitPoint()
+                    );
+                    pointConstraint c;
+                    c.applyConstraint(surfaceNormals[0]);
+                    c.applyConstraint(surfaceNormals[1]);
+                    edgeConstraints[featI][nearInfo.index()].append(c);
+
+                    // Store for later use
+                    patchAttraction[pointI] = nearInfo.hitPoint()-pt;
+                    patchConstraints[pointI] = c;
+
+                    // Dump
+                    if (featureEdgeStr.valid())
+                    {
+                        meshTools::writeOBJ(featureEdgeStr(), pt);
+                        featureEdgeVertI++;
+                        meshTools::writeOBJ
+                        (
+                            featureEdgeStr(),
+                            nearInfo.hitPoint()
+                        );
+                        featureEdgeVertI++;
+                        featureEdgeStr()
+                            << "l " << featureEdgeVertI-1 << ' '
+                            << featureEdgeVertI << nl;
+                    }
+                }
+                else
+                {
+                    if (missedEdgeStr.valid())
+                    {
+                        meshTools::writeOBJ(missedEdgeStr(), pt);
+                        missedVertI++;
+                        meshTools::writeOBJ
+                        (
+                            missedEdgeStr(),
+                            nearInfo.missPoint()
+                        );
+                        missedVertI++;
+                        missedEdgeStr()
+                            << "l " << missedVertI-1 << ' '
+                            << missedVertI << nl;
+                    }
+                }
+            }
+            else if (surfaceNormals.size() == 3)
+            {
+                // Calculate point from the faces.
+                plane pl0(surfacePoints[0], surfaceNormals[0]);
+                plane pl1(surfacePoints[1], surfaceNormals[1]);
+                plane pl2(surfacePoints[2], surfaceNormals[2]);
+                point cornerPt(pl0.planePlaneIntersect(pl1, pl2));
+                vector d = cornerPt - pt;
+
+                // Trim to snap distance
+                if (magSqr(d) > sqr(snapDist[pointI]))
+                {
+                    d *= Foam::sqrt(sqr(snapDist[pointI])/magSqr(d));
+                }
+
+
+                // Mark point on the nearest feature point.
+                const point estimatedPt(pt + d);
+
+                labelList nearPointFeat;
+                labelList nearPointIndex;
+                features.findNearestPoint
+                (
+                    pointField(1, estimatedPt),
+                    scalarField(1, sqr(snapDist[pointI])),
+                    nearPointFeat,
+                    nearPointIndex
+                );
+
+
+                if (nearPointIndex[0] != -1)
+                {
+                    label featI = nearPointFeat[0];
+                    label index = nearPointIndex[0];
+
+                    const treeDataPoint& shapes =
+                        features.pointTrees()[featI].shapes();
+                    label featPointI = shapes.pointLabels()[index];
+                    const point& featPt = shapes.points()[featPointI];
+
+                    pointAttractor[featI][featPointI] = pointI;
+                    pointConstraint& c = pointConstraints[featI][featPointI];
+                    c.applyConstraint(surfaceNormals[0]);
+                    c.applyConstraint(surfaceNormals[1]);
+                    c.applyConstraint(surfaceNormals[2]);
+
+                    // Store for later use
+                    patchAttraction[pointI] = featPt-pt;
+                    patchConstraints[pointI] = c;
+
+                    // Dump
+                    if (featurePointStr.valid())
+                    {
+                        meshTools::writeOBJ(featurePointStr(), pt);
+                        featurePointVertI++;
+                        meshTools::writeOBJ(featurePointStr(), featPt);
+                        featurePointVertI++;
+                        featurePointStr()
+                            << "l " << featurePointVertI-1 << ' '
+                            << featurePointVertI << nl;
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
+(
+    const label iter,
+    const scalar featureCos,
+
+    const indirectPrimitivePatch& pp,
+    const scalarField& snapDist,
+
+    //const vectorField& faceSurfaceNormal,
+    //const vectorField& faceDisp,
+    //const vectorField& faceRotation,
+    const List<List<point> >& pointFaceNormals,
+    const List<List<point> >& pointFaceDisp,
+    const List<List<point> >& pointFaceCentres,
+
+    vectorField& patchAttraction,
+    List<pointConstraint>& patchConstraints
+) const
+{
+    const refinementFeatures& features = meshRefiner_.features();
+
+    // Collect ordered attractions on feature edges
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Per feature, per feature-edge a list of attraction points and their
+    // originating vertex.
+    List<List<DynamicList<point> > > edgeAttractors(features.size());
+    List<List<DynamicList<pointConstraint> > > edgeConstraints(features.size());
+    forAll(features, featI)
+    {
+        label nFeatEdges = features[featI].edges().size();
+        edgeAttractors[featI].setSize(nFeatEdges);
+        edgeConstraints[featI].setSize(nFeatEdges);
+    }
+
+    // Per feature, per feature-point the pp point that is attracted to it.
+    // This list is only used to subset the feature-points that are actually
+    // used.
+    List<labelList> pointAttractor(features.size());
+    List<List<pointConstraint> > pointConstraints(features.size());
+    forAll(features, featI)
+    {
+        label nFeatPoints = features[featI].points().size();
+        pointAttractor[featI].setSize(nFeatPoints, -1);
+        pointConstraints[featI].setSize(nFeatPoints);
+    }
+
+    // Reverse: from pp point to nearest feature
+    vectorField allPatchAttraction(pp.nPoints(), vector::zero);
+    List<pointConstraint> allPatchConstraints(pp.nPoints());
+
+    determineFeatures
+    //determineAllFeatures
+    (
+        iter,
+        featureCos,
+
+        pp,
+        snapDist,
+
+        pointFaceNormals,
+        pointFaceDisp,
+        pointFaceCentres,
+
+        // Feature-point to pp point
+        pointAttractor,
+        pointConstraints,
+        // Feature-edge to pp point
+        edgeAttractors,
+        edgeConstraints,
+        // pp point to nearest feature
+        allPatchAttraction,
+        allPatchConstraints
+    );
+
+
+    // Find nearest mesh point to feature edge
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Reverse lookup : go through all edgeAttractors and find the
+    // nearest point on pp
+
+    indexedOctree<treeDataPoint> ppTree
+    (
+        treeDataPoint(pp.localPoints()),
+        treeBoundBox(pp.localPoints()), // overall search domain
+        8,                              // maxLevel
+        10,                             // leafsize
+        3.0                             // duplicity
+    );
+
+    // Per mesh point the point on nearest feature edge.
+    patchAttraction.setSize(pp.nPoints());
+    patchAttraction = vector::zero;
+    patchConstraints.setSize(pp.nPoints());
+    patchConstraints = pointConstraint();
+
+    forAll(edgeAttractors, featI)
+    {
+        const List<DynamicList<point> >& edgeAttr = edgeAttractors[featI];
+        const List<DynamicList<pointConstraint> >& edgeConstr =
+            edgeConstraints[featI];
+
+        forAll(edgeAttr, featEdgeI)
+        {
+            const DynamicList<point>& attr = edgeAttr[featEdgeI];
+            forAll(attr, i)
+            {
+                // Find nearest pp point
+
+                const point& featPt = attr[i];
+                pointIndexHit nearInfo = ppTree.findNearest(featPt, sqr(GREAT));
+
+                if (nearInfo.hit())
+                {
+                    label pointI = nearInfo.index();
+                    const point attraction = featPt-pp.localPoints()[pointI];
+
+                    // Check if this point is already being attracted. If so
+                    // override it only if nearer.
+                    if
+                    (
+                        patchConstraints[pointI].first() == 0
+                     || magSqr(attraction) < magSqr(patchAttraction[pointI])
+                    )
+                    {
+                        patchAttraction[pointI] = attraction;
+                        patchConstraints[pointI] = edgeConstr[featEdgeI][i];
+                    }
+                }
+                else
+                {
+                    WarningIn
+                    (
+                        "autoSnapDriver::featureAttractionUsingFeatureEdges(..)"
+                    )   << "Did not find pp point near " << featPt
+                        << endl;
+                }
+            }
+        }
+    }
+
+
+    // Different procs might have different patchAttraction,patchConstraints
+    // however these only contain geometric information, no topology
+    // so as long as we synchronise after overriding with feature points
+    // there is no problem, just possibly a small error.
+
+
+    // Find nearest mesh point to feature point
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // (overrides attraction to feature edge)
+    forAll(pointAttractor, featI)
+    {
+        const labelList& pointAttr = pointAttractor[featI];
+        const List<pointConstraint>& pointConstr = pointConstraints[featI];
+
+        forAll(pointAttr, featPointI)
+        {
+            if (pointAttr[featPointI] != -1)
+            {
+                const point& featPt = features[featI].points()
+                [
+                    featPointI
+                ];
+
+                // Find nearest pp point
+                pointIndexHit nearInfo = ppTree.findNearest(featPt, sqr(GREAT));
+
+                if (nearInfo.hit())
+                {
+                    label pointI = nearInfo.index();
+                    const point& pt = pp.localPoints()[pointI];
+                    const point attraction = featPt-pt;
+
+                    // - already attracted to feature edge : point always wins
+                    // - already attracted to feature point: nearest wins
+
+                    if (patchConstraints[pointI].first() <= 1)
+                    {
+                        patchAttraction[pointI] = attraction;
+                        patchConstraints[pointI] = pointConstr[featPointI];
+                    }
+                    else if (patchConstraints[pointI].first() == 2)
+                    {
+                        patchAttraction[pointI] = attraction;
+                        patchConstraints[pointI] = pointConstr[featPointI];
+                    }
+                    else if (patchConstraints[pointI].first() == 3)
+                    {
+                        // Only if nearer
+                        if
+                        (
+                            magSqr(attraction)
+                          < magSqr(patchAttraction[pointI])
+                        )
+                        {
+                            patchAttraction[pointI] = attraction;
+                            patchConstraints[pointI] = pointConstr[featPointI];
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+
+    // Dump
+    if (debug&meshRefinement::OBJINTERSECTIONS)
+    {
+        OFstream featureEdgeStr
+        (
+            meshRefiner_.mesh().time().path()
+          / "edgeAttractors_" + name(iter) + ".obj"
+        );
+        label featureEdgeVertI = 0;
+        Pout<< nl << "Dumping feature-edge attraction to "
+            << featureEdgeStr.name() << endl;
+
+        OFstream featurePointStr
+        (
+            meshRefiner_.mesh().time().path()
+          / "pointAttractors_" + name(iter) + ".obj"
+        );
+        label featurePointVertI = 0;
+        Pout<< nl << "Dumping feature-point attraction to "
+            << featurePointStr.name() << endl;
+
+        forAll(patchConstraints, pointI)
+        {
+            const point& pt = pp.localPoints()[pointI];
+
+            if (patchConstraints[pointI].first() == 2)
+            {
+                meshTools::writeOBJ(featureEdgeStr, pt);
+                featureEdgeVertI++;
+                meshTools::writeOBJ(featureEdgeStr, pt+patchAttraction[pointI]);
+                featureEdgeVertI++;
+                featureEdgeStr << "l " << featureEdgeVertI-1
+                    << ' ' << featureEdgeVertI << nl;
+            }
+            else if (patchConstraints[pointI].first() == 3)
+            {
+                meshTools::writeOBJ(featurePointStr, pt);
+                featurePointVertI++;
+                meshTools::writeOBJ
+                (
+                    featurePointStr,
+                    pt+patchAttraction[pointI]
+                );
+                featurePointVertI++;
+                featurePointStr << "l " << featurePointVertI-1
+                    << ' ' << featurePointVertI << nl;
+            }
+        }
+    }
+
+
+
+    // Snap edges to feature edges
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Walk existing edges and snap remaining ones (that are marked as
+    // feature edges in allPatchConstraints)
+
+    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() != 0)
+                    {
+                        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() == 0
+                         && allPatchConstraints[nbrPointI].first() != 0
+                        )
+                        {
+                            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;
+        }
+    }
+
+
+
+    // Avoid diagonal attraction
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Attract one of the non-diagonal points.
+
+    forAll(pp.localFaces(), faceI)
+    {
+        const face& f = pp.localFaces()[faceI];
+        // For now just detect any attraction. Improve this to look at
+        // actual attraction position and only if would form a problem add
+        // the non-diagonal point
+        if (f.size() == 4)
+        {
+            label nAttract = 0;
+            label firstAttract = -1;
+            forAll(f, fp)
+            {
+                label pointI = f[fp];
+                if (patchConstraints[pointI].first() == 2)
+                {
+                    nAttract++;
+                    if (firstAttract == -1)
+                    {
+                        firstAttract = fp;
+                    }
+                }
+            }
+            if (nAttract == 2)
+            {
+                label nextAttract = f.fcIndex(f.fcIndex(firstAttract));
+                label pointI = f[nextAttract];
+
+                if (patchConstraints[pointI].first() == 2)
+                {
+                    // Found two diagonal points that being attracted.
+                    // For now just attract my one to the average of those.
+                    const label i0 = f[firstAttract];
+                    const point pt0 = pp.localPoints()[i0]+patchAttraction[i0];
+                    const label i1 = f[nextAttract];
+                    const point pt1 = pp.localPoints()[i1]+patchAttraction[i1];
+                    const point mid = 0.5*(pt0+pt1);
+
+
+                    const scalar cosAngle = mag
+                    (
+                        patchConstraints[i0].second()
+                      & patchConstraints[i1].second()
+                    );
+
+                    //Pout<< "Found diagonal attraction at indices:"
+                    //    << firstAttract
+                    //    << " and " << nextAttract
+                    //    << " with cosAngle:" << cosAngle
+                    //    << " mid:" << mid << endl;
+
+                    if (cosAngle > featureCos)
+                    {
+                        // Add the nearest of the other two points as attractor
+                        label minFp = -1;
+                        scalar minDistSqr = GREAT;
+                        forAll(f, fp)
+                        {
+                            label pointI = f[fp];
+                            if (patchConstraints[pointI].first() == 0)
+                            {
+                                const point& pt = pp.localPoints()[pointI];
+                                scalar distSqr = magSqr(mid-pt);
+                                if (distSqr < minDistSqr)
+                                {
+                                    distSqr = minDistSqr;
+                                    minFp = fp;
+                                }
+                            }
+                        }
+                        if (minFp != -1)
+                        {
+                            label minPointI = f[minFp];
+                            patchAttraction[minPointI] =
+                                mid-pp.localPoints()[minPointI];
+                            patchConstraints[minPointI] =
+                                patchConstraints[f[firstAttract]];
+                        }
+                    }
+                    else
+                    {
+                        //Pout<< "Diagonal attractors at" << nl
+                        //    << "    pt0:" << pt0
+                        //    << "    constraint:"
+                        //    << patchConstraints[i0].second() << nl
+                        //    << "    pt1:" << pt1
+                        //    << "    constraint:"
+                        //    << patchConstraints[i1].second() << nl
+                        //    << "    make too large an angle:"
+                        //    <<  mag
+                        //        (
+                        //            patchConstraints[i0].second()
+                        //          & patchConstraints[i1].second()
+                        //        )
+                        //    << endl;
+                    }
+                }
+            }
+        }
+    }
+
+
+    if (debug&meshRefinement::OBJINTERSECTIONS)
+    {
+        dumpMove
+        (
+            meshRefiner_.mesh().time().path()
+          / "patchAttraction_" + name(iter) + ".obj",
+            pp.localPoints(),
+            pp.localPoints() + patchAttraction
+        );
+    }
+}
+
+
+// Correct for squeezing of face
+void Foam::autoSnapDriver::preventFaceSqueeze
+(
+    const label iter,
+    const scalar featureCos,
+
+    const indirectPrimitivePatch& pp,
+    const scalarField& snapDist,
+
+    vectorField& patchAttraction,
+    List<pointConstraint>& patchConstraints
+) const
+{
+    pointField points;
+    face singleF;
+    forAll(pp.localFaces(), faceI)
+    {
+        const face& f = pp.localFaces()[faceI];
+
+        if (f.size() != points.size())
+        {
+            points.setSize(f.size());
+            singleF.setSize(f.size());
+            for (label i = 0; i < f.size(); i++)
+            {
+                singleF[i] = i;
+            }
+        }
+        label nConstraints = 0;
+        forAll(f, fp)
+        {
+            label pointI = f[fp];
+            if (patchConstraints[pointI].first() != 0)
+            {
+                points[fp] = pp.localPoints()[pointI] + patchAttraction[pointI];
+                nConstraints++;
+            }
+            else
+            {
+                points[fp] = pp.localPoints()[pointI];
+            }
+        }
+
+        if (nConstraints == f.size())
+        {
+            scalar oldArea = f.mag(pp.localPoints());
+            scalar newArea = singleF.mag(points);
+            if (newArea < 0.1*oldArea)
+            {
+                // For now remove the point with largest distance
+                label maxFp = -1;
+                scalar maxS = -1;
+                forAll(f, fp)
+                {
+                    scalar s = magSqr(patchAttraction[f[fp]]);
+                    if (s > maxS)
+                    {
+                        maxS = s;
+                        maxFp = fp;
+                    }
+                }
+                if (maxFp != -1)
+                {
+                    label pointI = f[maxFp];
+                    // Lower attraction on pointI
+                    patchAttraction[pointI] *= 0.5;
+                }
+            }
+        }
+    }
+}
+
+
+Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
+(
+    const label iter,
+    const scalar featureCos,
+    const scalar featureAttract,
+    const scalarField& snapDist,
+    const vectorField& nearestDisp,
+    motionSmoother& meshMover
+) const
+{
+    Info<< "Calculating patchDisplacement as distance to nearest surface"
+        << " point ..." << endl;
+
+    const indirectPrimitivePatch& pp = meshMover.patch();
+    const pointField& localPoints = pp.localPoints();
+    const fvMesh& mesh = meshRefiner_.mesh();
+
+    // Displacement and orientation per pp face.
+    vectorField faceDisp(pp.size(), vector::zero);
+    vectorField faceSurfaceNormal(pp.size(), vector::zero);
+    vectorField faceRotation(pp.size(), vector::zero);
+
+    calcNearestFace
+    (
+        iter,
+        pp,
+        faceDisp,
+        faceSurfaceNormal,
+        faceRotation
+    );
+
+    //
+    // Interpolate face-wise data to points
+    //
+
+
+    // Displacement due to snapping of face centres
+    pointField patchDisp(localPoints.size(), vector::zero);
+    // Displacement due to rotation
+    pointField patchRotationDisp(localPoints.size(), vector::zero);
+
+    interpolateFaceToPoint
+    (
+        iter,
+        pp,
+        faceSurfaceNormal,
+        faceDisp,
+        faceRotation,
+
+        patchDisp,
+        patchRotationDisp
+    );
+
+    // Override patchDisp with exact nearest
+    {
+        pointField pointSurfaceNormal;
+        calcNearest(pp.localPoints(), patchDisp, pointSurfaceNormal);
+    }
+
+
+    // Collect (possibly remote) face-wise data on coupled points.
+    // - faceSurfaceNormal
+    // - faceDisp
+    // - faceRotation
+    // - faceCentres
+
+    // For now just get all surrounding face data. Expensive - should just
+    // store and sync data on coupled points only (see e.g PatchToolsNormals.C)
+
+    List<List<point> > pointFaceNormals(pp.nPoints());
+    List<List<point> > pointFaceDisp(pp.nPoints());
+    List<List<point> > pointFaceRotation(pp.nPoints());
+    List<List<point> > pointFaceCentres(pp.nPoints());
+
+    // Fill local data
+    forAll(pp.pointFaces(), pointI)
+    {
+        const labelList& pFaces = pp.pointFaces()[pointI];
+        List<point>& pNormals = pointFaceNormals[pointI];
+        pNormals.setSize(pFaces.size());
+        List<point>& pDisp = pointFaceDisp[pointI];
+        pDisp.setSize(pFaces.size());
+        List<point>& pRot = pointFaceRotation[pointI];
+        pRot.setSize(pFaces.size());
+        List<point>& pFc = pointFaceCentres[pointI];
+        pFc.setSize(pFaces.size());
+        forAll(pFaces, i)
+        {
+            pNormals[i] = faceSurfaceNormal[pFaces[i]];
+            pDisp[i] = faceDisp[pFaces[i]];
+            pRot[i] = faceRotation[pFaces[i]];
+            pFc[i] = pp.faceCentres()[pFaces[i]];
+        }
+    }
+
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        pointFaceNormals,
+        listPlusEqOp(),
+        List<point>(),
+        listTransform()
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        pointFaceDisp,
+        listPlusEqOp(),
+        List<point>(),
+        listTransform()
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        pointFaceRotation,
+        listPlusEqOp(),
+        List<point>(),
+        listTransform()
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        pointFaceCentres,
+        listPlusEqOp(),
+        List<point>(),
+        listTransform()
+    );
+
+
+
+//XXXXXXX
+//    // Indices of pp points that are coupled
+//    labelList coupledPointLabels(pp.nPoints());
+//    // Corresponding coupled point index
+//    labelList coupledPointAddr(pp.nPoints());
+//    {
+//        label coupledI = 0;
+//        forAll(pp.meshPoints(), pointI)
+//        {
+//            label meshPointI = pp.meshPoints()[pointI];
+//            Map<label>::const_iterator fnd = coupledPatchMP.find(meshPointI);
+//            if (fnd != coupledPatchMP.end())
+//            {
+//                coupledPointLabels[coupledI] = pointI;
+//                coupledPointAddr[coupledI] = fnd();
+//                coupledI++;
+//            }
+//        }
+//        coupledPointLabels.setSize(coupledI);
+//        coupledPointAddr.setSize(coupledI);
+//    }
+//
+//    List<List<point> > pointFaceNormals(map.constructSize());
+//    List<List<point> > pointFaceDisp(map.constructSize());
+//    List<List<point> > pointFaceRotation(map.constructSize());
+//    List<List<point> > pointFaceCentres(map.constructSize());
+//
+//    // Collect local data.
+//    forAll(coupledPointLabels, coupledI)
+//    {
+//        label pointI = coupledPointLabels[coupledI];
+//        label coupledPointI = coupledPointAddr[coupledI];
+//        const labelList& pFaces = pp.pointFaces()[pointI];
+//        List<point>& pNormals = pointFaceNormals[coupledPointI];
+//        pNormals.setSize(pFaces.size());
+//        List<point>& pDisp = pointFaceDisp[coupledPointI];
+//        pDisp.setSize(pFaces.size());
+//        List<point>& pRot = pointFaceRotation[coupledPointI];
+//        pRot.setSize(pFaces.size());
+//        List<point>& pFc = pointFaceCentres[coupledPointI];
+//        pFc.setSize(pFaces.size());
+//        forAll(pFaces, i)
+//        {
+//            pNormals[i] = faceSurfaceNormal[pFaces[i]];
+//            pDisp[i] = pointFaceDisp[pFaces[i]];
+//            pRot[i] = pointFaceRotation[pFaces[i]];
+//            pFc[i] = pp.faceCentres()[pFaces[i]];
+//        }
+//    }
+//
+//    // Pull remote data into local slots
+//XXXXXXX
+
+
+    // Nearest feature
+    vectorField patchAttraction(localPoints.size(), vector::zero);
+    // Constraints at feature
+    List<pointConstraint> patchConstraints(localPoints.size());
+    //featureAttractionUsingReconstruction
+    featureAttractionUsingFeatureEdges
+    (
+        iter,
+        featureCos,
+
+        pp,
+        snapDist,
+
+        pointFaceNormals,
+        pointFaceDisp,
+        pointFaceCentres,
+
+        patchAttraction,
+        patchConstraints
+    );
+
+    preventFaceSqueeze
+    (
+        iter,
+        featureCos,
+
+        pp,
+        snapDist,
+
+        patchAttraction,
+        patchConstraints
+    );
+
+
+    Info<< "Attraction:" << endl
+        << "     linear   : max:" << gMax(patchDisp)
+        << " avg:" << gAverage(patchDisp)
+        << endl
+        << "     rotation : max:" << gMax(patchRotationDisp)
+        << " avg:" << gAverage(patchRotationDisp)
+        << endl
+        << "     feature  : max:" << gMax(patchAttraction)
+        << " avg:" << gAverage(patchAttraction)
+        << endl;
+
+
+    // So now we have:
+    // - patchDisp          : point movement to go to nearest point on surface
+    //                       (either direct or through interpolation of
+    //                        face nearest)
+    // - patchRotationDisp  : point movement to align face normal with surface.
+    //
+    // - patchAttraction    : direct attraction to features
+    // - patchConstraints   : type of features
+
+    // Use any combination of patchDisp, patchRotationDisp and direct feature
+    // attraction.
+
+
+    // Bit of patchRotation
+    //patchDisp = 0.3*patchRotationDisp+0.7*patchDisp;
+    // Disable patchRotationDisp
+    patchRotationDisp = vector::zero;
+
+    // Mix with direct feature attraction
+    forAll(patchConstraints, pointI)
+    {
+        if (patchConstraints[pointI].first() != 0)
+        {
+            patchDisp[pointI] =
+                (1.0-featureAttract)*patchDisp[pointI]
+              + featureAttract*patchAttraction[pointI];
+        }
+    }
+
+    //dumpMove
+    //(
+    //    mesh.time().path()
+    //  / "linearPatchDisp_" + name(iter) + ".obj",
+    //    pp.localPoints(),
+    //    pp.localPoints() + patchDisp
+    //);
+    //dumpMove
+    //(
+    //    mesh.time().path()
+    //  / "rotationPatchDisp_" + name(iter) + ".obj",
+    //    pp.localPoints(),
+    //    pp.localPoints() + patchRotationDisp
+    //);
+
+
+
+    // Count
+    {
+        label nPlanar = 0;
+        label nEdge = 0;
+        label nPoint = 0;
+
+        forAll(patchConstraints, pointI)
+        {
+            if (patchConstraints[pointI].first() == 1)
+            {
+                nPlanar++;
+            }
+            else if (patchConstraints[pointI].first() == 2)
+            {
+                nEdge++;
+            }
+            else if (patchConstraints[pointI].first() == 3)
+            {
+                nPoint++;
+            }
+        }
+
+        label nTotPoints = returnReduce(pp.nPoints(), sumOp<label>());
+        reduce(nPlanar, sumOp<label>());
+        reduce(nEdge, sumOp<label>());
+        reduce(nPoint, sumOp<label>());
+        Info<< "Feature analysis : total points:"
+            << nTotPoints
+            << " attraction to :" << nl
+            << "    feature point   : " << nPoint << nl
+            << "    feature edge    : " << nEdge << nl
+            << "    nearest surface : " << nTotPoints-nPoint-nEdge
+            << " (rest)" << nl
+            << endl;
+    }
+
+
+
+    // Now we have the displacement per patch point to move onto the surface
+    // Split into tangential and normal direction.
+    // - start off with all non-constrained points following the constrained
+    //   ones since point normals not relevant.
+    // - finish with only tangential component smoothed.
+    // Note: tangential is most
+    // likely to come purely from face-centre snapping, not face rotation.
+    if (featureAttract < 1-0.001)
+    {
+        // 1. Smoothed all displacement
+        vectorField smoothedPatchDisp = patchDisp;
+        smoothAndConstrain
+        (
+            pp,
+            patchConstraints,
+            smoothedPatchDisp
+        );
+
+
+        // 2. Smoothed tangential component
+        vectorField tangPatchDisp = patchDisp;
+        tangPatchDisp -= (pp.pointNormals() & patchDisp) * pp.pointNormals();
+        smoothAndConstrain
+        (
+            pp,
+            patchConstraints,
+            tangPatchDisp
+        );
+
+        // Re-add normal component
+        tangPatchDisp += (pp.pointNormals() & patchDisp) * pp.pointNormals();
+
+        if (debug&meshRefinement::OBJINTERSECTIONS)
+        {
+            dumpMove
+            (
+                mesh.time().path()
+              / "tangPatchDispConstrained_" + name(iter) + ".obj",
+                pp.localPoints(),
+                pp.localPoints() + tangPatchDisp
+            );
+        }
+
+        patchDisp =
+             (1.0-featureAttract)*smoothedPatchDisp
+           + featureAttract*tangPatchDisp;
+    }
+
+
+
+    const scalar relax = featureAttract; //1.0;
+    patchDisp *= relax;
+
+
+    // Points on zones in one domain but only present as point on other
+    // will not do condition 2 on all. Sync explicitly.
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        patchDisp,
+        minMagSqrEqOp<point>(),         // combine op
+        vector(GREAT, GREAT, GREAT)     // null value (note: cannot use VGREAT)
+    );
+
+
+    // Check for displacement being outwards.
+    outwardsDisplacement(pp, patchDisp);
+
+    // Set initial distribution of displacement field (on patches) from
+    // patchDisp and make displacement consistent with b.c. on displacement
+    // pointVectorField.
+    meshMover.setDisplacement(patchDisp);
+
+    return patchDisp;
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C
index 4fb5a7778d..f9322efdb5 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,23 +27,14 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from dictionary
-Foam::snapParameters::snapParameters(const dictionary& dict, const label dummy)
-:
-    nSmoothPatch_(readLabel(dict.lookup("nSmoothPatch"))),
-    snapTol_(readScalar(dict.lookup("snapTol"))),
-    nSmoothDispl_(readLabel(dict.lookup("nSmoothDispl"))),
-    nSnap_(readLabel(dict.lookup("nSnap")))
-{}
-
-
 // Construct from dictionary
 Foam::snapParameters::snapParameters(const dictionary& dict)
 :
     nSmoothPatch_(readLabel(dict.lookup("nSmoothPatch"))),
     snapTol_(readScalar(dict.lookup("tolerance"))),
     nSmoothDispl_(readLabel(dict.lookup("nSolveIter"))),
-    nSnap_(readLabel(dict.lookup("nRelaxIter")))
+    nSnap_(readLabel(dict.lookup("nRelaxIter"))),
+    nFeatureSnap_(dict.lookupOrDefault("nFeatureSnapIter", -1))
 {}
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H
index 9be867e321..bd023565b1 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -61,6 +61,8 @@ class snapParameters
 
         const label nSnap_;
 
+        const label nFeatureSnap_;
+
 
     // Private Member Functions
 
@@ -75,10 +77,7 @@ public:
 
     // Constructors
 
-        //- Construct from dictionary - old syntax
-        snapParameters(const dictionary& dict, const label dummy);
-
-        //- Construct from dictionary - new syntax
+        //- Construct from dictionary
         snapParameters(const dictionary& dict);
 
 
@@ -114,6 +113,12 @@ public:
             {
                 return nSnap_;
             }
+
+            label nFeatureSnap() const
+            {
+                return nFeatureSnap_;
+            }
+
 };
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index 738a5cced3..e5065b91f0 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -31,6 +31,7 @@ License
 #include "Time.H"
 #include "refinementHistory.H"
 #include "refinementSurfaces.H"
+#include "refinementFeatures.H"
 #include "decompositionMethod.H"
 #include "regionSplit.H"
 #include "fvMeshDistribute.H"
@@ -497,6 +498,10 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
         mesh_.clearOut();
     }
 
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+    setInstance(mesh_.facesInstance());
+
     // Update local mesh data
     cellRemover.updateMesh(map);
 
@@ -857,6 +862,7 @@ Foam::meshRefinement::meshRefinement
     const scalar mergeDistance,
     const bool overwrite,
     const refinementSurfaces& surfaces,
+    const refinementFeatures& features,
     const shellSurfaces& shells
 )
 :
@@ -865,6 +871,7 @@ Foam::meshRefinement::meshRefinement
     overwrite_(overwrite),
     oldInstance_(mesh.pointsInstance()),
     surfaces_(surfaces),
+    features_(features),
     shells_(shells),
     meshCutter_
     (
@@ -1271,6 +1278,10 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
 
         // Update numbering of meshRefiner
         distribute(map);
+
+        // Set correct instance (for if overwrite)
+        mesh_.setInstance(timeName());
+        setInstance(mesh_.facesInstance());
     }
     return map;
 }
@@ -1449,6 +1460,8 @@ Foam::tmp<Foam::pointVectorField> Foam::meshRefinement::makeDisplacementField
         }
     }
 
+    // Note: time().timeName() instead of meshRefinement::timeName() since
+    // postprocessable field.
     tmp<pointVectorField> tfld
     (
         new pointVectorField
@@ -1456,7 +1469,7 @@ Foam::tmp<Foam::pointVectorField> Foam::meshRefinement::makeDisplacementField
             IOobject
             (
                 "pointDisplacement",
-                mesh.time().timeName(), //timeName(),
+                mesh.time().timeName(),
                 mesh,
                 IOobject::NO_READ,
                 IOobject::AUTO_WRITE
@@ -1930,10 +1943,6 @@ void Foam::meshRefinement::distribute(const mapDistributePolyMesh& map)
             pointMap.clear();
         }
     }
-
-    // If necessary reset the instance
-    mesh_.setInstance(timeName());
-    setInstance(mesh_.facesInstance());
 }
 
 
@@ -2059,10 +2068,6 @@ void Foam::meshRefinement::updateMesh
             data.transfer(newFaceData);
         }
     }
-
-    // If necessary reset the instance
-    mesh_.setInstance(timeName());
-    setInstance(mesh_.facesInstance());
 }
 
 
@@ -2118,12 +2123,32 @@ void Foam::meshRefinement::printMeshInfo(const bool debug, const string& msg)
             << "  points(local):" << mesh_.nPoints()
             << endl;
     }
-    else
+
     {
+        PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
+        label nMasterFaces = 0;
+        forAll(isMasterFace, i)
+        {
+            if (isMasterFace[i])
+            {
+                nMasterFaces++;
+            }
+        }
+
+        PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh_));
+        label nMasterPoints = 0;
+        forAll(isMasterPoint, i)
+        {
+            if (isMasterPoint[i])
+            {
+                nMasterPoints++;
+            }
+        }
+
         Info<< msg.c_str()
             << " : cells:" << pData.nTotalCells()
-            << "  faces:" << pData.nTotalFaces()
-            << "  points:" << pData.nTotalPoints()
+            << "  faces:" << returnReduce(nMasterFaces, sumOp<label>())
+            << "  points:" << returnReduce(nMasterPoints, sumOp<label>())
             << endl;
     }
 
@@ -2168,12 +2193,14 @@ Foam::word Foam::meshRefinement::timeName() const
 
 void Foam::meshRefinement::dumpRefinementLevel() const
 {
+    // Note: use time().timeName(), not meshRefinement::timeName()
+    // so as to dump the fields to 0, not to constant.
     volScalarField volRefLevel
     (
         IOobject
         (
             "cellLevel",
-            mesh_.time().timeName(),// Dump to current time, not to mesh inst
+            mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE,
@@ -2201,7 +2228,7 @@ void Foam::meshRefinement::dumpRefinementLevel() const
         IOobject
         (
             "pointLevel",
-            mesh_.time().timeName(),// Dump to current time, not to mesh inst
+            mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::NO_WRITE,
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 5a1d8ff821..0d7bf10bd2 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -62,13 +62,14 @@ class fvMesh;
 class mapDistributePolyMesh;
 class decompositionMethod;
 class refinementSurfaces;
+class refinementFeatures;
 class shellSurfaces;
 class removeCells;
-class featureEdgeMesh;
 class fvMeshDistribute;
 class searchableSurface;
 class regionSplit;
 class globalIndex;
+class removePoints;
 
 /*---------------------------------------------------------------------------*\
                            Class meshRefinement Declaration
@@ -118,6 +119,9 @@ private:
         //- All surface-intersection interaction
         const refinementSurfaces& surfaces_;
 
+        //- All feature-edge interaction
+        const refinementFeatures& features_;
+
         //- All shell-refinement interaction
         const shellSurfaces& shells_;
 
@@ -260,8 +264,6 @@ private:
             label markFeatureRefinement
             (
                 const point& keepPoint,
-                const PtrList<featureEdgeMesh>& featureMeshes,
-                const labelList& featureLevels,
                 const label nAllowRefine,
 
                 labelList& refineCell,
@@ -415,12 +417,6 @@ private:
                 const labelList& globalToPatch
             ) const;
 
-            ////- Initial test of marking faces using geometric information.
-            //labelList markFacesOnProblemCellsGeometric
-            //(
-            //    const dictionary& motionDict
-            //) const;
-
 
         // Baffle merging
 
@@ -501,6 +497,7 @@ public:
             const scalar mergeDistance,
             const bool overwrite,
             const refinementSurfaces&,
+            const refinementFeatures&,
             const shellSurfaces&
         );
 
@@ -542,6 +539,12 @@ public:
                 return surfaces_;
             }
 
+            //- reference to feature edge mesh
+            const refinementFeatures& features() const
+            {
+                return features_;
+            }
+
             //- reference to refinement shells (regions)
             const shellSurfaces& shells() const
             {
@@ -648,9 +651,6 @@ public:
                 const point& keepPoint,
                 const scalar curvature,
 
-                const PtrList<featureEdgeMesh>& featureMeshes,
-                const labelList& featureLevels,
-
                 const bool featureRefinement,
                 const bool internalRefinement,
                 const bool surfaceRefinement,
@@ -803,19 +803,45 @@ public:
                     const Map<label>& cellsToRestore
                 );
 
+            // Merging coplanar faces and edges
 
-            //- Merge faces on the same patch (usually from exposing refinement)
-            //  Returns global number of faces merged.
-            label mergePatchFaces
-            (
-                const scalar minCos,
-                const scalar concaveCos,
-                const labelList& patchIDs
-            );
+                label mergePatchFacesUndo
+                (
+                    const scalar minCos,
+                    const scalar concaveCos,
+                    const labelList& patchIDs,
+                    const dictionary& motionDict
+                );
 
-            //- Remove points not used by any face or points used
-            // by only two faces where the edges are in line
-            autoPtr<mapPolyMesh> mergeEdges(const scalar minCos);
+                autoPtr<mapPolyMesh> doRemovePoints
+                (
+                    removePoints& pointRemover,
+                    const boolList& pointCanBeDeleted
+                );
+
+                autoPtr<mapPolyMesh> doRestorePoints
+                (
+                    removePoints& pointRemover,
+                    const labelList& facesToRestore
+                );
+
+                labelList collectFaces
+                (
+                    const labelList& candidateFaces,
+                    const labelHashSet& set
+                ) const;
+
+                // Pick up faces of cells of faces in set.
+                labelList growFaceCellFace
+                (
+                    const labelHashSet& set
+                ) const;
+
+                label mergeEdgesUndo
+                (
+                    const scalar minCos,
+                    const dictionary& motionDict
+                );
 
 
         // Debug/IO
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index f3ab1e65c9..e318f37ed2 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -524,6 +524,10 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
         mesh_.clearOut();
     }
 
+
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+
     //- Redo the intersections on the newly create baffle faces. Note that
     //  this changes also the cell centre positions.
     faceSet baffledFacesSet(mesh_, "baffledFacesSet", 2*nBaffles);
@@ -644,7 +648,7 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::getDuplicateFaces
         }
         Pout<< "Writing duplicate faces (baffles) to faceSet "
             << duplicateFaceSet.name() << nl << endl;
-        duplicateFaceSet.instance() = mesh_.time().timeName();
+        duplicateFaceSet.instance() = timeName();
         duplicateFaceSet.write();
     }
 
@@ -724,9 +728,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createZoneBaffles
             if (debug)
             {
                 const_cast<Time&>(mesh_.time())++;
-                Pout<< "Writing baffled mesh to time "
-                    << mesh_.time().timeName() << endl;
-                mesh_.write();
+                Pout<< "Writing zone-baffled mesh to time " << timeName()
+                    << endl;
+                write(debug, mesh_.time().path()/"baffles");
             }
         }
         Info<< "Created " << nZoneFaces << " baffles in = "
@@ -746,14 +750,6 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
     const List<labelPair>& couples
 ) const
 {
-    // Construct addressing engine for all duplicate faces (only one
-    // for each pair)
-
-    // Duplicate faces in mesh labels (first face of each pair only)
-    // (reused later on to mark off filtered couples. see below)
-    labelList duplicateFaces(couples.size());
-
-
     // All duplicate faces on edge of the patch are to be merged.
     // So we count for all edges of duplicate faces how many duplicate
     // faces use them.
@@ -824,7 +820,8 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
 
 
     // Baffles which are not next to other boundaries and baffles will have
-    // value 2*1000000+2*1
+    // nBafflesPerEdge value 2*1000000+2*1 (from 2 boundary faces which are
+    // both baffle faces)
 
     List<labelPair> filteredCouples(couples.size());
     label filterI = 0;
@@ -967,6 +964,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeBaffles
         mesh_.clearOut();
     }
 
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+
     // Update intersections. Recalculate intersections on merged faces since
     // this seems to give problems? Note: should not be nessecary since
     // baffles preserve intersections from when they were created.
@@ -1789,7 +1789,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
                     problemTopo.insert(faceI);
                 }
             }
-            problemTopo.instance() = mesh_.time().timeName();
+            problemTopo.instance() = timeName();
             Pout<< "Dumping " << problemTopo.size()
                 << " problem faces to " << problemTopo.objectPath() << endl;
             problemTopo.write();
@@ -2246,6 +2246,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints()
         mesh_.clearOut();
     }
 
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+
     // Update intersections. Is mapping only (no faces created, positions stay
     // same) so no need to recalculate intersections.
     updateMesh(map, labelList(0));
@@ -2782,6 +2785,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
         mesh_.clearOut();
     }
 
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+
     // Print some stats (note: zones are synchronised)
     if (mesh_.cellZones().size() > 0)
     {
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementMerge.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementMerge.C
index e5a089c678..39a6833b55 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementMerge.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementMerge.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,42 +27,258 @@ License
 #include "combineFaces.H"
 #include "polyTopoChange.H"
 #include "removePoints.H"
+#include "faceSet.H"
+#include "Time.H"
+#include "motionSmoother.H"
+#include "syncTools.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-// Merge faces that are in-line.
-Foam::label Foam::meshRefinement::mergePatchFaces
+//// Merge faces that are in-line.
+//Foam::label Foam::meshRefinement::mergePatchFaces
+//(
+//    const scalar minCos,
+//    const scalar concaveCos,
+//    const labelList& patchIDs
+//)
+//{
+//    // Patch face merging engine
+//    combineFaces faceCombiner(mesh_);
+//
+//    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+//
+//    // Pick up all candidate cells on boundary
+//    labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
+//
+//    forAll(patchIDs, i)
+//    {
+//        label patchI = patchIDs[i];
+//
+//        const polyPatch& patch = patches[patchI];
+//
+//        if (!patch.coupled())
+//        {
+//            forAll(patch, i)
+//            {
+//                boundaryCells.insert(mesh_.faceOwner()[patch.start()+i]);
+//            }
+//        }
+//    }
+//
+//    // Get all sets of faces that can be merged
+//    labelListList mergeSets
+//    (
+//        faceCombiner.getMergeSets
+//        (
+//            minCos,
+//            concaveCos,
+//            boundaryCells
+//        )
+//    );
+//
+//    label nFaceSets = returnReduce(mergeSets.size(), sumOp<label>());
+//
+//    Info<< "mergePatchFaces : Merging " << nFaceSets
+//        << " sets of faces." << endl;
+//
+//    if (nFaceSets > 0)
+//    {
+//        // Topology changes container
+//        polyTopoChange meshMod(mesh_);
+//
+//        // Merge all faces of a set into the first face of the set. Remove
+//        // unused points.
+//        faceCombiner.setRefinement(mergeSets, meshMod);
+//
+//        // Change the mesh (no inflation)
+//        autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
+//
+//        // Update fields
+//        mesh_.updateMesh(map);
+//
+//        // Move mesh (since morphing does not do this)
+//        if (map().hasMotionPoints())
+//        {
+//            mesh_.movePoints(map().preMotionPoints());
+//        }
+//        else
+//        {
+//            // Delete mesh volumes. No other way to do this?
+//            mesh_.clearOut();
+//        }
+//
+//
+//        // Reset the instance for if in overwrite mode
+//        mesh_.setInstance(timeName());
+//
+//        faceCombiner.updateMesh(map);
+//
+//        // Get the kept faces that need to be recalculated.
+//        // Merging two boundary faces might shift the cell centre
+//        // (unless the faces are absolutely planar)
+//        labelHashSet retestFaces(6*mergeSets.size());
+//
+//        forAll(mergeSets, setI)
+//        {
+//            label oldMasterI = mergeSets[setI][0];
+//
+//            label faceI = map().reverseFaceMap()[oldMasterI];
+//
+//            // faceI is always uncoupled boundary face
+//            const cell& cFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
+//
+//            forAll(cFaces, i)
+//            {
+//                retestFaces.insert(cFaces[i]);
+//            }
+//        }
+//        updateMesh(map, retestFaces.toc());
+//    }
+//
+//
+//    return nFaceSets;
+//}
+//
+//
+//// Remove points not used by any face or points used by only two faces where
+//// the edges are in line
+//Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeEdges
+//(
+//    const scalar minCos
+//)
+//{
+//    // Point removal analysis engine
+//    removePoints pointRemover(mesh_);
+//
+//    // Count usage of points
+//    boolList pointCanBeDeleted;
+//    label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
+//
+//    Info<< "Removing " << nRemove
+//        << " straight edge points." << endl;
+//
+//    autoPtr<mapPolyMesh> map;
+//
+//    if (nRemove > 0)
+//    {
+//        // Save my local faces that will change. These changed faces might
+//        // cause a shift in the cell centre which needs to be retested.
+//        // Have to do this before changing mesh since point will be removed.
+//        labelHashSet retestOldFaces(nRemove / Pstream::nProcs());
+//
+//        {
+//            const faceList& faces = mesh_.faces();
+//
+//            forAll(faces, faceI)
+//            {
+//                const face& f = faces[faceI];
+//
+//                forAll(f, fp)
+//                {
+//                    if (pointCanBeDeleted[f[fp]])
+//                    {
+//                        retestOldFaces.insert(faceI);
+//                        break;
+//                    }
+//                }
+//            }
+//        }
+//
+//        // Topology changes container
+//        polyTopoChange meshMod(mesh_);
+//
+//        pointRemover.setRefinement(pointCanBeDeleted, meshMod);
+//
+//        // Change the mesh (no inflation)
+//        map = meshMod.changeMesh(mesh_, false, true);
+//
+//        // Update fields
+//        mesh_.updateMesh(map);
+//
+//        // Move mesh (since morphing does not do this)
+//        if (map().hasMotionPoints())
+//        {
+//            mesh_.movePoints(map().preMotionPoints());
+//        }
+//        else
+//        {
+//            // Delete mesh volumes. No other way to do this?
+//            mesh_.clearOut();
+//        }
+//
+//        // Reset the instance for if in overwrite mode
+//        mesh_.setInstance(timeName());
+//
+//        pointRemover.updateMesh(map);
+//
+//        // Get the kept faces that need to be recalculated.
+//        labelHashSet retestFaces(6*retestOldFaces.size());
+//
+//        const cellList& cells = mesh_.cells();
+//
+//        forAllConstIter(labelHashSet, retestOldFaces, iter)
+//        {
+//            label faceI = map().reverseFaceMap()[iter.key()];
+//
+//            const cell& ownFaces = cells[mesh_.faceOwner()[faceI]];
+//
+//            forAll(ownFaces, i)
+//            {
+//                retestFaces.insert(ownFaces[i]);
+//            }
+//
+//            if (mesh_.isInternalFace(faceI))
+//            {
+//                const cell& neiFaces = cells[mesh_.faceNeighbour()[faceI]];
+//
+//                forAll(neiFaces, i)
+//                {
+//                    retestFaces.insert(neiFaces[i]);
+//                }
+//            }
+//        }
+//        updateMesh(map, retestFaces.toc());
+//    }
+//
+//    return map;
+//}
+
+
+Foam::label Foam::meshRefinement::mergePatchFacesUndo
 (
     const scalar minCos,
     const scalar concaveCos,
-    const labelList& patchIDs
+    const labelList& patchIDs,
+    const dictionary& motionDict
 )
 {
     // Patch face merging engine
-    combineFaces faceCombiner(mesh_);
-
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    combineFaces faceCombiner(mesh_, true);
 
     // Pick up all candidate cells on boundary
     labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
 
-    forAll(patchIDs, i)
     {
-        label patchI = patchIDs[i];
+        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-        const polyPatch& patch = patches[patchI];
-
-        if (!patch.coupled())
+        forAll(patchIDs, i)
         {
-            forAll(patch, i)
+            label patchI = patchIDs[i];
+
+            const polyPatch& patch = patches[patchI];
+
+            if (!patch.coupled())
             {
-                boundaryCells.insert(mesh_.faceOwner()[patch.start()+i]);
+                forAll(patch, i)
+                {
+                    boundaryCells.insert(mesh_.faceOwner()[patch.start()+i]);
+                }
             }
         }
     }
 
     // Get all sets of faces that can be merged
-    labelListList mergeSets
+    labelListList allFaceSets
     (
         faceCombiner.getMergeSets
         (
@@ -72,19 +288,44 @@ Foam::label Foam::meshRefinement::mergePatchFaces
         )
     );
 
-    label nFaceSets = returnReduce(mergeSets.size(), sumOp<label>());
+    label nFaceSets = returnReduce(allFaceSets.size(), sumOp<label>());
 
-    Info<< "mergePatchFaces : Merging " << nFaceSets
-        << " sets of faces." << endl;
+    Info<< "Merging " << nFaceSets << " sets of faces." << nl << endl;
 
     if (nFaceSets > 0)
     {
+        if (debug)
+        {
+            faceSet allSets(mesh_, "allFaceSets", allFaceSets.size());
+            forAll(allFaceSets, setI)
+            {
+                forAll(allFaceSets[setI], i)
+                {
+                    allSets.insert(allFaceSets[setI][i]);
+                }
+            }
+            Pout<< "Writing all faces to be merged to set "
+                << allSets.objectPath() << endl;
+            allSets.instance() = timeName();
+            allSets.write();
+
+            const_cast<Time&>(mesh_.time())++;
+        }
+
+
         // Topology changes container
         polyTopoChange meshMod(mesh_);
 
-        // Merge all faces of a set into the first face of the set. Remove
-        // unused points.
-        faceCombiner.setRefinement(mergeSets, meshMod);
+        // Merge all faces of a set into the first face of the set.
+        faceCombiner.setRefinement(allFaceSets, meshMod);
+
+        // Experimental: store data for all the points that have been deleted
+        storeData
+        (
+            faceCombiner.savedPointLabels(),    // points to store
+            labelList(0),                       // faces to store
+            labelList(0)                        // cells to store
+        );
 
         // Change the mesh (no inflation)
         autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
@@ -99,138 +340,538 @@ Foam::label Foam::meshRefinement::mergePatchFaces
         }
         else
         {
-            // Delete mesh volumes. No other way to do this?
+            // Delete mesh volumes.
             mesh_.clearOut();
         }
 
+        // Reset the instance for if in overwrite mode
+        mesh_.setInstance(timeName());
+
         faceCombiner.updateMesh(map);
 
-        // Get the kept faces that need to be recalculated.
-        // Merging two boundary faces might shift the cell centre
-        // (unless the faces are absolutely planar)
-        labelHashSet retestFaces(6*mergeSets.size());
+        updateMesh(map, labelList(0));
 
-        forAll(mergeSets, setI)
+        if (debug)
         {
-            label oldMasterI = mergeSets[setI][0];
-
-            label faceI = map().reverseFaceMap()[oldMasterI];
-
-            // faceI is always uncoupled boundary face
-            const cell& cFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
-
-            forAll(cFaces, i)
-            {
-                retestFaces.insert(cFaces[i]);
-            }
+            Pout<< "Writing initial merged-faces mesh to time "
+                << timeName() << nl << endl;
+            write();
         }
-        updateMesh(map, retestFaces.toc());
-    }
 
+        for (label iteration = 0; iteration < 100; iteration++)
+        {
+            Info<< nl
+                << "Undo iteration " << iteration << nl
+                << "----------------" << endl;
+
+
+            // Check mesh for errors
+            // ~~~~~~~~~~~~~~~~~~~~~
+
+            faceSet errorFaces
+            (
+                mesh_,
+                "errorFaces",
+                mesh_.nFaces()-mesh_.nInternalFaces()
+            );
+            bool hasErrors = motionSmoother::checkMesh
+            (
+                false,  // report
+                mesh_,
+                motionDict,
+                errorFaces
+            );
+
+            //if (checkEdgeConnectivity)
+            //{
+            //    Info<< "Checking edge-face connectivity (duplicate faces"
+            //        << " or non-consecutive shared vertices)" << endl;
+            //
+            //    label nOldSize = errorFaces.size();
+            //
+            //    hasErrors =
+            //        mesh_.checkFaceFaces
+            //        (
+            //            false,
+            //            &errorFaces
+            //        )
+            //     || hasErrors;
+            //
+            //    Info<< "Detected additional "
+            //        <<  returnReduce
+            //            (
+            //                errorFaces.size() - nOldSize,
+            //                sumOp<label>()
+            //            )
+            //        << " faces with illegal face-face connectivity"
+            //        << endl;
+            //}
+
+            if (!hasErrors)
+            {
+                break;
+            }
+
+
+            if (debug)
+            {
+                errorFaces.instance() = timeName();
+                Pout<< "Writing all faces in error to faceSet "
+                    << errorFaces.objectPath() << nl << endl;
+                errorFaces.write();
+            }
+
+
+            // Check any master cells for using any of the error faces
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            DynamicList<label> mastersToRestore(allFaceSets.size());
+
+            forAll(allFaceSets, setI)
+            {
+                label masterFaceI = faceCombiner.masterFace()[setI];
+
+                if (masterFaceI != -1)
+                {
+                    label masterCellII = mesh_.faceOwner()[masterFaceI];
+
+                    const cell& cFaces = mesh_.cells()[masterCellII];
+
+                    forAll(cFaces, i)
+                    {
+                        if (errorFaces.found(cFaces[i]))
+                        {
+                            mastersToRestore.append(masterFaceI);
+                            break;
+                        }
+                    }
+                }
+            }
+            mastersToRestore.shrink();
+
+            label nRestore = returnReduce
+            (
+                mastersToRestore.size(),
+                sumOp<label>()
+            );
+
+            Info<< "Masters that need to be restored:"
+                << nRestore << endl;
+
+            if (debug)
+            {
+                faceSet restoreSet(mesh_, "mastersToRestore", mastersToRestore);
+                restoreSet.instance() = timeName();
+                Pout<< "Writing all " << mastersToRestore.size()
+                    << " masterfaces to be restored to set "
+                    << restoreSet.objectPath() << endl;
+                restoreSet.write();
+            }
+
+
+            if (nRestore == 0)
+            {
+                break;
+            }
+
+
+            // Undo
+            // ~~~~
+
+            if (debug)
+            {
+                const_cast<Time&>(mesh_.time())++;
+            }
+
+            // Topology changes container
+            polyTopoChange meshMod(mesh_);
+
+            // Merge all faces of a set into the first face of the set.
+            // Experimental:mark all points/faces/cells that have been restored.
+            Map<label> restoredPoints(0);
+            Map<label> restoredFaces(0);
+            Map<label> restoredCells(0);
+
+            faceCombiner.setUnrefinement
+            (
+                mastersToRestore,
+                meshMod,
+                restoredPoints,
+                restoredFaces,
+                restoredCells
+            );
+
+            // Change the mesh (no inflation)
+            autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
+
+            // Update fields
+            mesh_.updateMesh(map);
+
+            // Move mesh (since morphing does not do this)
+            if (map().hasMotionPoints())
+            {
+                mesh_.movePoints(map().preMotionPoints());
+            }
+            else
+            {
+                // Delete mesh volumes.
+                mesh_.clearOut();
+            }
+
+
+            // Reset the instance for if in overwrite mode
+            mesh_.setInstance(timeName());
+
+            faceCombiner.updateMesh(map);
+
+            // Renumber restore maps
+            inplaceMapKey(map().reversePointMap(), restoredPoints);
+            inplaceMapKey(map().reverseFaceMap(), restoredFaces);
+            inplaceMapKey(map().reverseCellMap(), restoredCells);
+
+            // Experimental:restore all points/face/cells in maps
+            updateMesh
+            (
+                map,
+                labelList(0),       // changedFaces
+                restoredPoints,
+                restoredFaces,
+                restoredCells
+            );
+
+            if (debug)
+            {
+                Pout<< "Writing merged-faces mesh to time "
+                    << timeName() << nl << endl;
+                write();
+            }
+
+            Info<< endl;
+        }
+    }
+    else
+    {
+        Info<< "No faces merged ..." << endl;
+    }
 
     return nFaceSets;
 }
 
 
-// Remove points not used by any face or points used by only two faces where
-// the edges are in line
-Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeEdges
+// Remove points. pointCanBeDeleted is parallel synchronised.
+Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemovePoints
 (
-    const scalar minCos
+    removePoints& pointRemover,
+    const boolList& pointCanBeDeleted
 )
 {
-    // Point removal analysis engine
-    removePoints pointRemover(mesh_);
+    // Topology changes container
+    polyTopoChange meshMod(mesh_);
 
-    // Count usage of points
-    boolList pointCanBeDeleted;
-    label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
+    pointRemover.setRefinement(pointCanBeDeleted, meshMod);
 
-    Info<< "Removing " << nRemove
-        << " straight edge points." << endl;
+    // Change the mesh (no inflation)
+    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
 
-    autoPtr<mapPolyMesh> map;
+    // Update fields
+    mesh_.updateMesh(map);
 
-    if (nRemove > 0)
+    // Move mesh (since morphing does not do this)
+    if (map().hasMotionPoints())
     {
-        // Save my local faces that will change. These changed faces might
-        // cause a shift in the cell centre which needs to be retested.
-        // Have to do this before changing mesh since point will be removed.
-        labelHashSet retestOldFaces(nRemove / Pstream::nProcs());
-
-        {
-            const faceList& faces = mesh_.faces();
-
-            forAll(faces, faceI)
-            {
-                const face& f = faces[faceI];
-
-                forAll(f, fp)
-                {
-                    if (pointCanBeDeleted[f[fp]])
-                    {
-                        retestOldFaces.insert(faceI);
-                        break;
-                    }
-                }
-            }
-        }
-
-        // Topology changes container
-        polyTopoChange meshMod(mesh_);
-
-        pointRemover.setRefinement(pointCanBeDeleted, meshMod);
-
-        // Change the mesh (no inflation)
-        map = meshMod.changeMesh(mesh_, false, true);
-
-        // Update fields
-        mesh_.updateMesh(map);
-
-        // Move mesh (since morphing does not do this)
-        if (map().hasMotionPoints())
-        {
-            mesh_.movePoints(map().preMotionPoints());
-        }
-        else
-        {
-            // Delete mesh volumes. No other way to do this?
-            mesh_.clearOut();
-        }
-
-        pointRemover.updateMesh(map);
-
-        // Get the kept faces that need to be recalculated.
-        labelHashSet retestFaces(6*retestOldFaces.size());
-
-        const cellList& cells = mesh_.cells();
-
-        forAllConstIter(labelHashSet, retestOldFaces, iter)
-        {
-            label faceI = map().reverseFaceMap()[iter.key()];
-
-            const cell& ownFaces = cells[mesh_.faceOwner()[faceI]];
-
-            forAll(ownFaces, i)
-            {
-                retestFaces.insert(ownFaces[i]);
-            }
-
-            if (mesh_.isInternalFace(faceI))
-            {
-                const cell& neiFaces = cells[mesh_.faceNeighbour()[faceI]];
-
-                forAll(neiFaces, i)
-                {
-                    retestFaces.insert(neiFaces[i]);
-                }
-            }
-        }
-        updateMesh(map, retestFaces.toc());
+        mesh_.movePoints(map().preMotionPoints());
     }
+    else
+    {
+        // Delete mesh volumes.
+        mesh_.clearOut();
+    }
+
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+
+    pointRemover.updateMesh(map);
+    updateMesh(map, labelList(0));
 
     return map;
 }
 
 
+// Restore faces (which contain removed points)
+Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRestorePoints
+(
+    removePoints& pointRemover,
+    const labelList& facesToRestore
+)
+{
+    // Topology changes container
+    polyTopoChange meshMod(mesh_);
+
+    // Determine sets of points and faces to restore
+    labelList localFaces, localPoints;
+    pointRemover.getUnrefimentSet
+    (
+        facesToRestore,
+        localFaces,
+        localPoints
+    );
+
+    // Undo the changes on the faces that are in error.
+    pointRemover.setUnrefinement
+    (
+        localFaces,
+        localPoints,
+        meshMod
+    );
+
+    // Change the mesh (no inflation)
+    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
+
+    // Update fields
+    mesh_.updateMesh(map);
+
+    // Move mesh (since morphing does not do this)
+    if (map().hasMotionPoints())
+    {
+        mesh_.movePoints(map().preMotionPoints());
+    }
+    else
+    {
+        // Delete mesh volumes.
+        mesh_.clearOut();
+    }
+
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+
+    pointRemover.updateMesh(map);
+    updateMesh(map, labelList(0));
+
+    return map;
+}
+
+
+// Collect all faces that are both in candidateFaces and in the set.
+// If coupled face also collects the coupled face.
+Foam::labelList Foam::meshRefinement::collectFaces
+(
+    const labelList& candidateFaces,
+    const labelHashSet& set
+) const
+{
+    // Has face been selected?
+    boolList selected(mesh_.nFaces(), false);
+
+    forAll(candidateFaces, i)
+    {
+        label faceI = candidateFaces[i];
+
+        if (set.found(faceI))
+        {
+            selected[faceI] = true;
+        }
+    }
+    syncTools::syncFaceList
+    (
+        mesh_,
+        selected,
+        orEqOp<bool>()      // combine operator
+    );
+
+    labelList selectedFaces(findIndices(selected, true));
+
+    return selectedFaces;
+}
+
+
+// Pick up faces of cells of faces in set.
+Foam::labelList Foam::meshRefinement::growFaceCellFace
+(
+    const labelHashSet& set
+) const
+{
+    boolList selected(mesh_.nFaces(), false);
+
+    forAllConstIter(faceSet, set, iter)
+    {
+        label faceI = iter.key();
+
+        label own = mesh_.faceOwner()[faceI];
+
+        const cell& ownFaces = mesh_.cells()[own];
+        forAll(ownFaces, ownFaceI)
+        {
+            selected[ownFaces[ownFaceI]] = true;
+        }
+
+        if (mesh_.isInternalFace(faceI))
+        {
+            label nbr = mesh_.faceNeighbour()[faceI];
+
+            const cell& nbrFaces = mesh_.cells()[nbr];
+            forAll(nbrFaces, nbrFaceI)
+            {
+                selected[nbrFaces[nbrFaceI]] = true;
+            }
+        }
+    }
+    syncTools::syncFaceList
+    (
+        mesh_,
+        selected,
+        orEqOp<bool>()      // combine operator
+    );
+    return findIndices(selected, true);
+}
+
+
+// Remove points not used by any face or points used by only two faces where
+// the edges are in line
+Foam::label Foam::meshRefinement::mergeEdgesUndo
+(
+    const scalar minCos,
+    const dictionary& motionDict
+)
+{
+    Info<< nl
+        << "Merging all points on surface that" << nl
+        << "- are used by only two boundary faces and" << nl
+        << "- make an angle with a cosine of more than " << minCos
+        << "." << nl << endl;
+
+    // Point removal analysis engine with undo
+    removePoints pointRemover(mesh_, true);
+
+    // Count usage of points
+    boolList pointCanBeDeleted;
+    label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
+
+    if (nRemove > 0)
+    {
+        Info<< "Removing " << nRemove
+            << " straight edge points ..." << nl << endl;
+
+        // Remove points
+        // ~~~~~~~~~~~~~
+
+        doRemovePoints(pointRemover, pointCanBeDeleted);
+
+
+        for (label iteration = 0; iteration < 100; iteration++)
+        {
+            Info<< nl
+                << "Undo iteration " << iteration << nl
+                << "----------------" << endl;
+
+
+            // Check mesh for errors
+            // ~~~~~~~~~~~~~~~~~~~~~
+
+            faceSet errorFaces
+            (
+                mesh_,
+                "errorFaces",
+                mesh_.nFaces()-mesh_.nInternalFaces()
+            );
+            bool hasErrors = motionSmoother::checkMesh
+            (
+                false,  // report
+                mesh_,
+                motionDict,
+                errorFaces
+            );
+            //if (checkEdgeConnectivity)
+            //{
+            //    Info<< "Checking edge-face connectivity (duplicate faces"
+            //        << " or non-consecutive shared vertices)" << endl;
+            //
+            //    label nOldSize = errorFaces.size();
+            //
+            //    hasErrors =
+            //        mesh_.checkFaceFaces
+            //        (
+            //            false,
+            //            &errorFaces
+            //        )
+            //     || hasErrors;
+            //
+            //    Info<< "Detected additional "
+            //        << returnReduce(errorFaces.size()-nOldSize,sumOp<label>())
+            //        << " faces with illegal face-face connectivity"
+            //        << endl;
+            //}
+
+            if (!hasErrors)
+            {
+                break;
+            }
+
+            if (debug)
+            {
+                errorFaces.instance() = timeName();
+                Pout<< "**Writing all faces in error to faceSet "
+                    << errorFaces.objectPath() << nl << endl;
+                errorFaces.write();
+            }
+
+            labelList masterErrorFaces
+            (
+                collectFaces
+                (
+                    pointRemover.savedFaceLabels(),
+                    errorFaces
+                )
+            );
+
+            label n = returnReduce(masterErrorFaces.size(), sumOp<label>());
+
+            Info<< "Detected " << n
+                << " error faces on boundaries that have been merged."
+                << " These will be restored to their original faces." << nl
+                << endl;
+
+            if (n == 0)
+            {
+                if (hasErrors)
+                {
+                    Info<< "Detected "
+                        << returnReduce(errorFaces.size(), sumOp<label>())
+                        << " error faces in mesh."
+                        << " Restoring neighbours of faces in error." << nl
+                        << endl;
+
+                    labelList expandedErrorFaces
+                    (
+                        growFaceCellFace
+                        (
+                            errorFaces
+                        )
+                    );
+
+                    doRestorePoints(pointRemover, expandedErrorFaces);
+                }
+
+                break;
+            }
+
+            doRestorePoints(pointRemover, masterErrorFaces);
+        }
+
+        if (debug)
+        {
+            const_cast<Time&>(mesh_.time())++;
+            Pout<< "Writing merged-edges mesh to time "
+                << timeName() << nl << endl;
+            write();
+        }
+    }
+    else
+    {
+        Info<< "No straight edges simplified and no points removed ..." << endl;
+    }
+
+    return nRemove;
+}
+
+
 // ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
index 1ad85acab6..ab1dda4c1c 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -202,7 +202,7 @@ Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
     if (debug)
     {
         faceSet fSet(mesh_, "edgeConnectedFaces", candidateFaces);
-        fSet.instance() = mesh_.time().timeName();
+        fSet.instance() = timeName();
         Pout<< "Writing " << fSet.size()
             << " with problematic topology to faceSet "
             << fSet.objectPath() << endl;
@@ -265,7 +265,7 @@ Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
 
     if (debug)
     {
-        perpFaces.instance() = mesh_.time().timeName();
+        perpFaces.instance() = timeName();
         Pout<< "Writing " << perpFaces.size()
             << " faces that are perpendicular to the surface to set "
             << perpFaces.objectPath() << endl;
@@ -427,6 +427,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 
     // Count of faces marked for baffling
     label nBaffleFaces = 0;
+    PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
 
     // Count of faces not baffled since would not cause a collapse
     label nPrevented = 0;
@@ -483,7 +484,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
         if (debug)
         {
             cellSet problemCellSet(mesh_, "problemCells", problemCells.toc());
-            problemCellSet.instance() = mesh_.time().timeName();
+            problemCellSet.instance() = timeName();
             Pout<< "Writing " << problemCellSet.size()
                 << " cells that are edge connected to coarser cell to set "
                 << problemCellSet.objectPath() << endl;
@@ -918,7 +919,10 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
                         else
                         {
                             facePatch[faceI] = getBafflePatch(facePatch, faceI);
-                            nBaffleFaces++;
+                            if (isMasterFace[faceI])
+                            {
+                                nBaffleFaces++;
+                            }
 
                             // Do NOT update boundary data since this would grow
                             // blocked faces across gaps.
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index a7cc79ee7a..20cbd076e1 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -28,6 +28,7 @@ License
 #include "syncTools.H"
 #include "Time.H"
 #include "refinementSurfaces.H"
+#include "refinementFeatures.H"
 #include "shellSurfaces.H"
 #include "faceSet.H"
 #include "decompositionMethod.H"
@@ -53,6 +54,11 @@ Foam::labelList Foam::meshRefinement::getChangedFaces
     const polyMesh& mesh = map.mesh();
 
     labelList changedFaces;
+
+    // For reporting: number of masterFaces changed
+    label nMasterChanged = 0;
+    PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
+
     {
         // Mark any face on a cell which has been added or changed
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -182,36 +188,26 @@ Foam::labelList Foam::meshRefinement::getChangedFaces
         // Now we have in changedFace marked all affected faces. Pack.
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        label nChanged = 0;
+        changedFaces = findIndices(changedFace, true);
+
+        // Count changed master faces.
+        nMasterChanged = 0;
 
         forAll(changedFace, faceI)
         {
-            if (changedFace[faceI])
+            if (changedFace[faceI] && isMasterFace[faceI])
             {
-                nChanged++;
+                nMasterChanged++;
             }
         }
 
-        changedFaces.setSize(nChanged);
-        nChanged = 0;
-
-        forAll(changedFace, faceI)
-        {
-            if (changedFace[faceI])
-            {
-                changedFaces[nChanged++] = faceI;
-            }
-        }
     }
 
-    label nChangedFaces = changedFaces.size();
-    reduce(nChangedFaces, sumOp<label>());
-
     if (debug)
     {
         Pout<< "getChangedFaces : Detected "
             << " local:" << changedFaces.size()
-            << " global:" << nChangedFaces
+            << " global:" << returnReduce(nMasterChanged, sumOp<label>())
             << " changed faces out of " << mesh.globalData().nTotalFaces()
             << endl;
 
@@ -252,8 +248,6 @@ bool Foam::meshRefinement::markForRefine
 Foam::label Foam::meshRefinement::markFeatureRefinement
 (
     const point& keepPoint,
-    const PtrList<featureEdgeMesh>& featureMeshes,
-    const labelList& featureLevels,
     const label nAllowRefine,
 
     labelList& refineCell,
@@ -287,9 +281,11 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
 
     if (cellI != -1)
     {
-        forAll(featureMeshes, featI)
+        forAll(features_, featI)
         {
-            const featureEdgeMesh& featureMesh = featureMeshes[featI];
+            const featureEdgeMesh& featureMesh = features_[featI];
+            const label featureLevel = features_.levels()[featI];
+
             const labelListList& pointEdges = featureMesh.pointEdges();
 
             forAll(pointEdges, pointI)
@@ -315,7 +311,7 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
                             tetFaceI,
                             tetPtI,
                             featureMesh.points()[pointI],   // endpos
-                            featureLevels[featI],           // level
+                            featureLevel,                   // level
                             featI,                          // featureMesh
                             pointI                          // end point
                         )
@@ -339,11 +335,11 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
     maxFeatureLevel = -1;
 
     // Whether edge has been visited.
-    List<PackedBoolList> featureEdgeVisited(featureMeshes.size());
+    List<PackedBoolList> featureEdgeVisited(features_.size());
 
-    forAll(featureMeshes, featI)
+    forAll(features_, featI)
     {
-        featureEdgeVisited[featI].setSize(featureMeshes[featI].edges().size());
+        featureEdgeVisited[featI].setSize(features_[featI].edges().size());
         featureEdgeVisited[featI] = 0u;
     }
 
@@ -359,7 +355,7 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
             label featI = tp.i();
             label pointI = tp.j();
 
-            const featureEdgeMesh& featureMesh = featureMeshes[featI];
+            const featureEdgeMesh& featureMesh = features_[featI];
             const labelList& pEdges = featureMesh.pointEdges()[pointI];
 
             // Particle now at pointI. Check connected edges to see which one
@@ -1068,9 +1064,6 @@ Foam::labelList Foam::meshRefinement::refineCandidates
     const point& keepPoint,
     const scalar curvature,
 
-    const PtrList<featureEdgeMesh>& featureMeshes,
-    const labelList& featureLevels,
-
     const bool featureRefinement,
     const bool internalRefinement,
     const bool surfaceRefinement,
@@ -1134,8 +1127,6 @@ Foam::labelList Foam::meshRefinement::refineCandidates
             label nFeatures = markFeatureRefinement
             (
                 keepPoint,
-                featureMeshes,
-                featureLevels,
                 nAllowRefine,
 
                 refineCell,
@@ -1251,6 +1242,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::refine
         mesh_.clearOut();
     }
 
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+
     // Update intersection info
     updateMesh(map, getChangedFaces(map, cellsToRefine));
 
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
new file mode 100644
index 0000000000..ec97e51fe3
--- /dev/null
+++ b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
@@ -0,0 +1,233 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "refinementFeatures.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::refinementFeatures::refinementFeatures
+(
+    const objectRegistry& io,
+    const PtrList<dictionary>& featDicts
+)
+:
+    PtrList<featureEdgeMesh>(featDicts.size()),
+    levels_(featDicts.size()),
+    edgeTrees_(featDicts.size()),
+    pointTrees_(featDicts.size())
+{
+    // Read features
+
+    forAll(featDicts, i)
+    {
+        const dictionary& dict = featDicts[i];
+
+        fileName featFileName(dict.lookup("file"));
+
+        set
+        (
+            i,
+            new featureEdgeMesh
+            (
+                IOobject
+                (
+                    featFileName,                       // name
+                    io.time().constant(),               // instance
+                    "triSurface",                       // local
+                    io.time(),                          // registry
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE,
+                    false
+                )
+            )
+        );
+
+        const featureEdgeMesh& eMesh = operator[](i);
+
+        //eMesh.mergePoints(meshRefiner_.mergeDistance());
+        levels_[i] = readLabel(dict.lookup("level"));
+
+        Info<< "Refinement level " << levels_[i]
+            << " for all cells crossed by feature " << featFileName
+            << " (" << eMesh.points().size() << " points, "
+            << eMesh.edges().size() << " edges)." << endl;
+    }
+
+
+    // Search engines
+
+    forAll(*this, i)
+    {
+        const featureEdgeMesh& eMesh = operator[](i);
+        const pointField& points = eMesh.points();
+        const edgeList& edges = eMesh.edges();
+
+        // Calculate bb of all points
+        const treeBoundBox bb(points);
+
+        edgeTrees_.set
+        (
+            i,
+            new indexedOctree<treeDataEdge>
+            (
+                treeDataEdge
+                (
+                    false,                  // do not cache bb
+                    edges,
+                    points,
+                    identity(edges.size())
+                ),
+                bb,     // overall search domain
+                8,      // maxLevel
+                10,     // leafsize
+                3.0     // duplicity
+            )
+        );
+
+
+        // Detect feature points from edges.
+        const labelListList& pointEdges = eMesh.pointEdges();
+        DynamicList<label> featurePoints;
+        forAll(pointEdges, pointI)
+        {
+            if (pointEdges[pointI].size() > 2)
+            {
+                featurePoints.append(pointI);
+            }
+        }
+
+        Info<< "Detected " << featurePoints.size()
+            << " featurePoints out of " << points.size() << endl;
+
+        pointTrees_.set
+        (
+            i,
+            new indexedOctree<treeDataPoint>
+            (
+                treeDataPoint(points, featurePoints),
+                bb,     // overall search domain
+                8,      // maxLevel
+                10,     // leafsize
+                3.0     // duplicity
+            )
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::refinementFeatures::findNearestEdge
+(
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+    labelList& nearFeature,
+    List<pointIndexHit>& nearInfo
+) const
+{
+    nearFeature.setSize(samples.size());
+    nearFeature = -1;
+    nearInfo.setSize(samples.size());
+
+    forAll(edgeTrees_, featI)
+    {
+        const indexedOctree<treeDataEdge>& tree = edgeTrees_[featI];
+        forAll(samples, sampleI)
+        {
+            const point& sample = samples[sampleI];
+
+            scalar distSqr;
+            if (nearInfo[sampleI].hit())
+            {
+                distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
+            }
+            else
+            {
+                distSqr = nearestDistSqr[sampleI];
+            }
+
+            pointIndexHit info = tree.findNearest(sample, distSqr);
+
+            if (info.hit())
+            {
+                nearInfo[sampleI] = info;
+                nearFeature[sampleI] = featI;
+            }
+        }
+    }
+}
+
+
+void Foam::refinementFeatures::findNearestPoint
+(
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+    labelList& nearFeature,
+    labelList& nearIndex
+) const
+{
+    nearFeature.setSize(samples.size());
+    nearFeature = -1;
+    nearIndex.setSize(samples.size());
+    nearIndex = -1;
+
+    forAll(pointTrees_, featI)
+    {
+        const indexedOctree<treeDataPoint>& tree = pointTrees_[featI];
+        forAll(samples, sampleI)
+        {
+            const point& sample = samples[sampleI];
+
+            scalar distSqr;
+            if (nearFeature[sampleI] != -1)
+            {
+                label nearFeatI = nearFeature[sampleI];
+                const indexedOctree<treeDataPoint>& nearTree =
+                    pointTrees_[nearFeatI];
+                label featPointI =
+                    nearTree.shapes().pointLabels()[nearIndex[sampleI]];
+                const point& featPt =
+                    operator[](nearFeatI).points()[featPointI];
+                distSqr = magSqr(featPt-sample);
+            }
+            else
+            {
+                distSqr = nearestDistSqr[sampleI];
+            }
+
+            pointIndexHit info = tree.findNearest(sample, distSqr);
+
+            if (info.hit())
+            {
+                nearFeature[sampleI] = featI;
+                nearIndex[sampleI] = info.index();
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H
new file mode 100644
index 0000000000..d13cad8d0f
--- /dev/null
+++ b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H
@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::refinementFeatures
+
+Description
+    Encapsulates queries for features.
+
+SourceFiles
+    refinementFeatures.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef refinementFeatures_H
+#define refinementFeatures_H
+
+#include "featureEdgeMesh.H"
+#include "indexedOctree.H"
+#include "treeDataEdge.H"
+#include "treeDataPoint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class refinementFeatures Declaration
+\*---------------------------------------------------------------------------*/
+
+class refinementFeatures
+:
+    public PtrList<featureEdgeMesh>
+{
+private:
+
+    // Private data
+
+        //- Refinement levels
+        labelList levels_;
+
+        //- Edge
+        PtrList<indexedOctree<treeDataEdge> > edgeTrees_;
+
+        //- Features points
+        PtrList<indexedOctree<treeDataPoint> > pointTrees_;
+
+
+    // Private Member Functions
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        refinementFeatures
+        (
+            const objectRegistry& io,
+            const PtrList<dictionary>& featDicts
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            const labelList& levels() const
+            {
+                return levels_;
+            }
+
+            const PtrList<indexedOctree<treeDataEdge> >& edgeTrees() const
+            {
+                return edgeTrees_;
+            }
+
+            const PtrList<indexedOctree<treeDataPoint> >& pointTrees() const
+            {
+                return pointTrees_;
+            }
+
+
+        // Query
+
+            //- Find nearest point on nearest feature edge
+            void findNearestEdge
+            (
+                const pointField& samples,
+                const scalarField& nearestDistSqr,
+                labelList& nearFeature,
+                List<pointIndexHit>& nearInfo
+            ) const;
+
+            //- Find nearest feature point. Is an index into feature points
+            //  which itself is an index into the edgeMesh points.
+            //  So the point index is
+            //      pointTrees()[nearFeature].shapes().pointLabels()[nearIndex]
+            //  Wip.
+            void findNearestPoint
+            (
+                const pointField& samples,
+                const scalarField& nearestDistSqr,
+                labelList& nearFeature,
+                labelList& nearIndex
+            ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //