diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index 1021e16432..d6002dc40f 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) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -663,18 +663,11 @@ void Foam::autoRefineDriver::mergePatchFaces
     const dictionary& motionDict
 )
 {
-    const fvMesh& mesh = meshRefiner_.mesh();
-
     Info<< nl
         << "Merge refined boundary faces" << nl
         << "----------------------------" << nl
         << endl;
 
-    if (debug)
-    {
-        const_cast<Time&>(mesh.time())++;
-    }
-
     meshRefiner_.mergePatchFacesUndo
     (
         Foam::cos(degToRad(45.0)),
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index 948930467a..7f612246f4 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -1465,6 +1465,7 @@ void Foam::autoSnapDriver::doSnap
             {
                 disp = calcNearestSurfaceFeature
                 (
+                    snapParams,
                     iter,
                     featureCos,
                     scalar(iter+1)/nFeatIter,
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index 35d0822d2a..9431f9b751 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -111,17 +111,43 @@ class autoSnapDriver
 
             // Feature line snapping
 
+                //- Is point on two feature edges that make a largish angle?
+                bool isFeaturePoint
+                (
+                    const scalar featureCos,
+                    const indirectPrimitivePatch& pp,
+                    const PackedBoolList& isFeatureEdge,
+                    const label pointI
+                ) const;
+
                 void smoothAndConstrain
                 (
                     const indirectPrimitivePatch& pp,
                     const List<pointConstraint>& constraints,
                     vectorField& disp
                 ) const;
-                //void calcNearest
+                void smoothAndConstrain2
+                (
+                    const bool applyConstraints,
+                    const indirectPrimitivePatch& pp,
+                    const List<pointConstraint>& constraints,
+                    vectorField& disp
+                ) const;
+                void calcNearest
+                (
+                    const label iter,
+                    const indirectPrimitivePatch& pp,
+                    vectorField& pointDisp,
+                    vectorField& pointSurfaceNormal,
+                    vectorField& pointRotation
+                ) const;
+                //void calcNearestFace
                 //(
-                //    const pointField& points,
-                //    vectorField& disp,
-                //    vectorField& surfaceNormal
+                //    const label iter,
+                //    const indirectPrimitivePatch& pp,
+                //    vectorField& faceDisp,
+                //    vectorField& faceSurfaceNormal,
+                //    vectorField& faceRotation
                 //) const;
                 void calcNearestFace
                 (
@@ -129,19 +155,18 @@ class autoSnapDriver
                     const indirectPrimitivePatch& pp,
                     vectorField& faceDisp,
                     vectorField& faceSurfaceNormal,
+                    labelList& faceSurfaceRegion,
                     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 List<List<point> >& pointFaceDisp,
+                    const List<List<point> >& pointFaceRotation,
+                    const List<List<point> >& pointFaceCentres,
+                    vectorField& patchDisp
+                    //vectorField& patchRotationDisp
                 ) const;
                 void correctAttraction
                 (
@@ -152,6 +177,14 @@ class autoSnapDriver
                     const point& pt,
                     vector& edgeOffset  // offset from pt to point on edge
                 ) const;
+
+                //- Return hit if on multiple points
+                pointIndexHit findMultiPatchPoint
+                (
+                    const point& pt,
+                    const labelList& patchIDs,
+                    const List<point>& faceCentres
+                ) const;
                 void binFeatureFace
                 (
                     const label iter,
@@ -186,6 +219,25 @@ class autoSnapDriver
                     DynamicList<label>& surfaceCount
                 ) const;
 
+
+                void featureAttractionUsingReconstruction
+                (
+                    const label iter,
+                    const scalar featureCos,
+
+                    const indirectPrimitivePatch& pp,
+                    const scalarField& snapDist,
+                    const label pointI,
+
+                    const List<List<point> >& pointFaceSurfNormals,
+                    const List<List<point> >& pointFaceDisp,
+                    const List<List<point> >& pointFaceCentres,
+                    const labelListList& pointFacePatchID,
+
+                    vector& patchAttraction,
+                    pointConstraint& patchConstraint
+                ) const;
+
                 void featureAttractionUsingReconstruction
                 (
                     const label iter,
@@ -196,31 +248,12 @@ class autoSnapDriver
                     const List<List<point> >& pointFaceNormals,
                     const List<List<point> >& pointFaceDisp,
                     const List<List<point> >& pointFaceCentres,
+                    const labelListList& pointFacePatchID,
 
                     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,
@@ -232,6 +265,7 @@ class autoSnapDriver
                     const List<List<point> >& pointFaceNormals,
                     const List<List<point> >& pointFaceDisp,
                     const List<List<point> >& pointFaceCentres,
+                    const labelListList& pointFacePatchID,
 
                     List<labelList>& pointAttractor,
                     List<List<pointConstraint> >& pointConstraints,
@@ -242,6 +276,47 @@ class autoSnapDriver
                     List<pointConstraint>& patchConstraints
                 ) const;
 
+                //- Find point on nearest feature edge (within searchDist).
+                //  Return point and feature
+                //  and store feature-edge to mesh-point and vice versa
+                pointIndexHit findNearFeatureEdge
+                (
+                    const indirectPrimitivePatch& pp,
+                    const scalarField& snapDist,
+                    const label pointI,
+                    const point& estimatedPt,
+
+                    label& featI,
+                    List<List<DynamicList<point> > >&,
+                    List<List<DynamicList<pointConstraint> > >&,
+                    vectorField&,
+                    List<pointConstraint>&
+                ) const;
+
+                //- Find nearest feature point (within searchDist).
+                //  Return feature point
+                //  and store feature-point to mesh-point and vice versa.
+                //  If another mesh point already referring to this feature
+                //  point and further away, reset that one to a near feature
+                //  edge (using findNearFeatureEdge above)
+                labelPair findNearFeaturePoint
+                (
+                    const indirectPrimitivePatch& pp,
+                    const scalarField& snapDist,
+                    const label pointI,
+                    const point& estimatedPt,
+
+                    // 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;
+
                 void featureAttractionUsingFeatureEdges
                 (
                     const label iter,
@@ -270,6 +345,7 @@ class autoSnapDriver
 
                 vectorField calcNearestSurfaceFeature
                 (
+                    const snapParameters& snapParams,
                     const label iter,
                     const scalar featureCos,
                     const scalar featureAttract,
@@ -342,7 +418,6 @@ public:
                 motionSmoother& meshMover
             ) const;
 
-
             //- Smooth the displacement field to the internal.
             void smoothDisplacement
             (
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index 950849e15e..e526f8faf9 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -37,6 +37,7 @@ License
 #include "featureEdgeMesh.H"
 #include "treeDataPoint.H"
 #include "indexedOctree.H"
+#include "snapParameters.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -95,6 +96,67 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+bool Foam::autoSnapDriver::isFeaturePoint
+(
+    const scalar featureCos,
+    const indirectPrimitivePatch& pp,
+    const PackedBoolList& isFeatureEdge,
+    const label pointI
+) const
+{
+    const pointField& points = pp.localPoints();
+    const edgeList& edges = pp.edges();
+    const labelList& pEdges = pp.pointEdges()[pointI];
+
+    label nFeatEdges = 0;
+
+    forAll(pEdges, i)
+    {
+        if (isFeatureEdge[pEdges[i]])
+        {
+            nFeatEdges++;
+
+            for (label j = i+1; j < pEdges.size(); j++)
+            {
+                if (isFeatureEdge[pEdges[j]])
+                {
+                    const edge& eI = edges[pEdges[i]];
+                    const edge& eJ = edges[pEdges[j]];
+
+                    const point& p = points[pointI];
+                    const point& pI = points[eI.otherVertex(pointI)];
+                    const point& pJ = points[eJ.otherVertex(pointI)];
+
+                    vector vI = p-pI;
+                    scalar vIMag = mag(vI);
+
+                    vector vJ = pJ-p;
+                    scalar vJMag = mag(vJ);
+
+                    if
+                    (
+                        vIMag > SMALL
+                     && vJMag > SMALL
+                     && ((vI/vIMag & vJ/vJMag) < featureCos)
+                    )
+                    {
+                        return true;
+                    }
+                }
+            }
+        }
+    }
+
+    if (nFeatEdges == 1)
+    {
+        // End of feature-edge string
+        return true;
+    }
+
+    return false;
+}
+
+
 void Foam::autoSnapDriver::smoothAndConstrain
 (
     const indirectPrimitivePatch& pp,
@@ -102,12 +164,100 @@ void Foam::autoSnapDriver::smoothAndConstrain
     vectorField& disp
 ) const
 {
+    const fvMesh& mesh = meshRefiner_.mesh();
+
     for (label avgIter = 0; avgIter < 20; avgIter++)
     {
         // Calculate average displacement of neighbours
-        vectorField dispAvg(pp.nPoints(), vector::zero);
+        // - unconstrained (i.e. surface) points use average of all
+        //   neighbouring points
+        // - from testing it has been observed that it is not beneficial
+        //   to have edge constrained points use average of all edge or point
+        //   constrained neighbours since they're already attracted to
+        //   the nearest point on the feature.
+        //   Having them attract to point-constrained neighbours does not
+        //   make sense either since there is usually just one of them so
+        //   it severely distorts it.
+        // - same for feature points. They are already attracted to the
+        //   nearest feature point.
+
+        vectorField dispSum(pp.nPoints(), vector::zero);
+        labelList dispCount(pp.nPoints(), 0);
 
         const labelListList& pointEdges = pp.pointEdges();
+        const edgeList& edges = pp.edges();
+
+        forAll(pointEdges, pointI)
+        {
+            const labelList& pEdges = pointEdges[pointI];
+
+            label nConstraints = constraints[pointI].first();
+
+            if (nConstraints <= 1)
+            {
+                forAll(pEdges, i)
+                {
+                    label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
+                    if (constraints[nbrPointI].first() >= nConstraints)
+                    {
+                        dispSum[pointI] += disp[nbrPointI];
+                        dispCount[pointI]++;
+                    }
+                }
+            }
+        }
+
+        syncTools::syncPointList
+        (
+            mesh,
+            pp.meshPoints(),
+            dispSum,
+            plusEqOp<point>(),
+            vector::zero,
+            mapDistribute::transform()
+        );
+        syncTools::syncPointList
+        (
+            mesh,
+            pp.meshPoints(),
+            dispCount,
+            plusEqOp<label>(),
+            0,
+            mapDistribute::transform()
+        );
+
+        // Constraints
+        forAll(constraints, pointI)
+        {
+            if (dispCount[pointI] > 0)
+            {
+                // Mix my displacement with neighbours' displacement
+                disp[pointI] =
+                    0.5
+                   *(disp[pointI] + dispSum[pointI]/dispCount[pointI]);
+            }
+        }
+    }
+}
+//XXXXXX
+//TODO: make proper parallel so coupled edges don't have double influence
+void Foam::autoSnapDriver::smoothAndConstrain2
+(
+    const bool applyConstraints,
+    const indirectPrimitivePatch& pp,
+    const List<pointConstraint>& constraints,
+    vectorField& disp
+) const
+{
+    const fvMesh& mesh = meshRefiner_.mesh();
+
+    for (label avgIter = 0; avgIter < 20; avgIter++)
+    {
+        vectorField dispSum(pp.nPoints(), vector::zero);
+        labelList dispCount(pp.nPoints(), 0);
+
+        const labelListList& pointEdges = pp.pointEdges();
+        const edgeList& edges = pp.edges();
 
         forAll(pointEdges, pointI)
         {
@@ -115,24 +265,54 @@ void Foam::autoSnapDriver::smoothAndConstrain
 
             forAll(pEdges, i)
             {
-                const edge& e = pp.edges()[pEdges[i]];
-                label nbrPointI = e.otherVertex(pointI);
-                dispAvg[pointI] += disp[nbrPointI];
+                label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
+                dispSum[pointI] += disp[nbrPointI];
+                dispCount[pointI]++;
             }
-            dispAvg[pointI] /= pEdges.size();
         }
 
+        syncTools::syncPointList
+        (
+            mesh,
+            pp.meshPoints(),
+            dispSum,
+            plusEqOp<point>(),
+            vector::zero,
+            mapDistribute::transform()
+        );
+        syncTools::syncPointList
+        (
+            mesh,
+            pp.meshPoints(),
+            dispCount,
+            plusEqOp<label>(),
+            0,
+            mapDistribute::transform()
+        );
+
         // Constraints
         forAll(constraints, pointI)
         {
-            if (constraints[pointI].first() == 0)
+            if (dispCount[pointI] > 0)// && constraints[pointI].first() <= 1)
             {
                 // Mix my displacement with neighbours' displacement
-                disp[pointI] = 0.5*disp[pointI] + 0.5*dispAvg[pointI];
+                disp[pointI] =
+                    0.5
+                   *(disp[pointI] + dispSum[pointI]/dispCount[pointI]);
+
+                if (applyConstraints)
+                {
+                    disp[pointI] = transform
+                    (
+                        constraints[pointI].constraintTransformation(),
+                        disp[pointI]
+                    );
+                }
             }
         }
     }
 }
+//XXXXXX
 
 
 void Foam::autoSnapDriver::calcNearestFace
@@ -141,6 +321,7 @@ void Foam::autoSnapDriver::calcNearestFace
     const indirectPrimitivePatch& pp,
     vectorField& faceDisp,
     vectorField& faceSurfaceNormal,
+    labelList& faceSurfaceRegion,
     vectorField& faceRotation
 ) const
 {
@@ -152,6 +333,8 @@ void Foam::autoSnapDriver::calcNearestFace
     faceDisp = vector::zero;
     faceSurfaceNormal.setSize(pp.size());
     faceSurfaceNormal = vector::zero;
+    faceSurfaceRegion.setSize(pp.size());
+    faceSurfaceRegion = -1;
 
     // Divide surfaces into zoned and unzoned
     labelList zonedSurfaces = surfaces.getNamedSurfaces();
@@ -236,6 +419,11 @@ void Foam::autoSnapDriver::calcNearestFace
                 label faceI = ppFaces[hitI];
                 faceDisp[faceI] = hitInfo[hitI].hitPoint() - fc[hitI];
                 faceSurfaceNormal[faceI] = hitNormal[hitI];
+                faceSurfaceRegion[faceI] = surfaces.globalRegion
+                (
+                    hitSurface[hitI],
+                    hitRegion[hitI]
+                );
             }
         }
     }
@@ -255,8 +443,8 @@ void Foam::autoSnapDriver::calcNearestFace
             meshFaces.append(pp.addressing()[i]);
         }
     }
-    //Pout<< "Found " << ppFaces.size() << " unzoned faces out of " << pp.size()
-    //    << endl;
+    //Pout<< "Found " << ppFaces.size() << " unzoned faces out of "
+    //   << pp.size() << endl;
 
     pointField fc
     (
@@ -289,6 +477,7 @@ void Foam::autoSnapDriver::calcNearestFace
             label faceI = ppFaces[hitI];
             faceDisp[faceI] = hitInfo[hitI].hitPoint() - fc[hitI];
             faceSurfaceNormal[faceI] = hitNormal[hitI];
+            faceSurfaceRegion[faceI] = hitRegion[hitI];
         }
     }
 
@@ -305,7 +494,7 @@ void Foam::autoSnapDriver::calcNearestFace
         // Note: extend to >180 degrees checking
         faceRotation[faceI] =
             pp.faceNormals()[faceI]
-          ^ -faceSurfaceNormal[faceI];
+          ^ faceSurfaceNormal[faceI];
     }
 
     if (debug&meshRefinement::OBJINTERSECTIONS)
@@ -369,6 +558,79 @@ void Foam::autoSnapDriver::correctAttraction
 }
 
 
+Foam::pointIndexHit Foam::autoSnapDriver::findMultiPatchPoint
+(
+    const point& pt,
+    const labelList& patchIDs,
+    const List<point>& faceCentres
+) const
+{
+    // Determine if multiple patchIDs
+    if (patchIDs.size())
+    {
+        label patch0 = patchIDs[0];
+
+        for (label i = 1; i < patchIDs.size(); i++)
+        {
+            if (patchIDs[i] != patch0)
+            {
+                return pointIndexHit(true, pt, labelMax);
+            }
+        }
+    }
+    return pointIndexHit(false, vector::zero, labelMax);
+}
+////XXXXXXXX
+//void Foam::autoSnapDriver::attractMultiPatchPoint
+//(
+//    const label iter,
+//    const scalar featureCos,
+//
+//    const indirectPrimitivePatch& pp,
+//    const scalarField& snapDist,
+//    const label pointI,
+//
+//    const List<List<point> >& pointFaceSurfNormals,
+//    const labelListList& pointFaceSurfaceRegion,
+//    const List<List<point> >& pointFaceDisp,
+//    const List<List<point> >& pointFaceCentres,
+//    const labelListList& pointFacePatchID,
+//
+//    vector& patchAttraction,
+//    pointConstraint& patchConstraint
+//) const
+//{
+//    // Collect
+//
+//        );
+//
+//        if
+//        (
+//            (constraint.first() > patchConstraints[pointI].first())
+//         || (magSqr(attraction) < magSqr(patchAttraction[pointI]))
+//        )
+//        {
+//            patchAttraction[pointI] = attraction;
+//            patchConstraints[pointI] = constraint;
+//
+//            // Check the number of directions
+//            if (patchConstraints[pointI].first() == 1)
+//            {
+//                // Flat surface. Check for different patchIDs
+//                pointIndexHit multiPatchPt
+//                (
+//                    findMultiPatchPoint
+//                    (
+//                        pt,
+//                        pointFacePatchID[pointI],
+//                        pointFaceCentres[pointI]
+//                    )
+//                );
+//                if (multiPatchPt.hit())
+//                {
+//                    // Behave like when having two surface normals so
+////XXXXXXXX
+
 void Foam::autoSnapDriver::binFeatureFace
 (
     const label iter,
@@ -452,10 +714,6 @@ void Foam::autoSnapDriver::binFeatureFace
                     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);
@@ -464,15 +722,6 @@ void Foam::autoSnapDriver::binFeatureFace
             }
         }
     }
-    //else
-    //{
-    //    Pout<< "binFeatureFace : for point:" << pp.localPoints()[pointI]
-    //        << " found face-nearest point:"
-    //        << fc + faceDisp
-    //        << " at distance:" << mag(faceDisp)
-    //        << " further away than:" << snapDist
-    //        << endl;
-    //}
 }
 
 
@@ -523,6 +772,131 @@ void Foam::autoSnapDriver::binFeatureFaces
 }
 
 
+void Foam::autoSnapDriver::featureAttractionUsingReconstruction
+(
+    const label iter,
+    const scalar featureCos,
+
+    const indirectPrimitivePatch& pp,
+    const scalarField& snapDist,
+    const label pointI,
+
+    const List<List<point> >& pointFaceSurfNormals,
+    const List<List<point> >& pointFaceDisp,
+    const List<List<point> >& pointFaceCentres,
+    const labelListList& pointFacePatchID,
+
+    vector& patchAttraction,
+    pointConstraint& patchConstraint
+) const
+{
+    patchAttraction = vector::zero;
+    patchConstraint = pointConstraint();
+
+    // Collect all different directions
+    DynamicList<point> surfacePoints(4);
+    DynamicList<vector> surfaceNormals(4);
+    DynamicList<label> surfaceCount(4);
+
+    binFeatureFaces
+    (
+        iter,
+        featureCos,
+
+        pp,
+        snapDist,
+        pointI,
+
+        pointFaceSurfNormals,
+        pointFaceDisp,
+        pointFaceCentres,
+
+        surfacePoints,
+        surfaceNormals,
+        surfaceCount
+    );
+
+    const point& pt = pp.localPoints()[pointI];
+
+    // Check the number of directions
+    if (surfaceNormals.size() == 1)
+    {
+        // Normal distance to plane
+        vector d =
+            ((surfacePoints[0]-pt) & surfaceNormals[0])
+           *surfaceNormals[0];
+
+        // Trim to snap distance
+        if (magSqr(d) > sqr(snapDist[pointI]))
+        {
+            d *= Foam::sqrt(sqr(snapDist[pointI])/magSqr(d));
+        }
+
+        patchAttraction = d;
+
+        // Store constraints
+        patchConstraint.applyConstraint(surfaceNormals[0]);
+    }
+    else 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 = d;
+
+        // Store constraints
+        patchConstraint.applyConstraint(surfaceNormals[0]);
+        patchConstraint.applyConstraint(surfaceNormals[1]);
+    }
+    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));
+        }
+
+        patchAttraction = d;
+
+        // Store constraints
+        patchConstraint.applyConstraint(surfaceNormals[0]);
+        patchConstraint.applyConstraint(surfaceNormals[1]);
+        patchConstraint.applyConstraint(surfaceNormals[2]);
+    }
+}
+
+
 // Special version that calculates attraction in one go
 void Foam::autoSnapDriver::featureAttractionUsingReconstruction
 (
@@ -535,54 +909,49 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
     const List<List<point> >& pointFaceSurfNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
+    const labelListList& pointFacePatchID,
 
     vectorField& patchAttraction,
     List<pointConstraint>& patchConstraints
 ) const
 {
-    autoPtr<OFstream> featureEdgeStr;
-    label featureEdgeVertI = 0;
-    autoPtr<OFstream> featurePointStr;
-    label featurePointVertI = 0;
+    autoPtr<OFstream> feStr;
+    label feVertI = 0;
+    autoPtr<OFstream> fpStr;
+    label fpVertI = 0;
 
     if (debug&meshRefinement::OBJINTERSECTIONS)
     {
-        featureEdgeStr.reset
+        feStr.reset
         (
             new OFstream
             (
                 meshRefiner_.mesh().time().path()
-              / "featureEdge_" + name(iter) + ".obj"
+              / "implicitFeatureEdge_" + name(iter) + ".obj"
             )
         );
-        Pout<< nl << "Dumping feature-edge direction to "
-            << featureEdgeStr().name() << endl;
+        Pout<< "Dumping implicit feature-edge direction to "
+            << feStr().name() << endl;
 
-        featurePointStr.reset
+        fpStr.reset
         (
             new OFstream
             (
                 meshRefiner_.mesh().time().path()
-              / "featurePoint_" + name(iter) + ".obj"
+              / "implicitFeaturePoint_" + name(iter) + ".obj"
             )
         );
-        Pout<< nl << "Dumping feature-point direction to "
-            << featurePointStr().name() << endl;
+        Pout<< "Dumping implicit feature-point direction to "
+            << fpStr().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;
+        vector attraction = vector::zero;
+        pointConstraint constraint;
 
-        binFeatureFaces
+        featureAttractionUsingReconstruction
         (
             iter,
             featureCos,
@@ -594,137 +963,104 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
             pointFaceSurfNormals,
             pointFaceDisp,
             pointFaceCentres,
+            pointFacePatchID,
 
-            surfacePoints,
-            surfaceNormals,
-            surfaceCount
+            attraction,
+            constraint
         );
 
-        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
+        if
         (
-            meshRefiner_.mesh().time().path()
-          / "patchAttraction_" + name(iter) + ".obj",
-            pp.localPoints(),
-            pp.localPoints() + patchAttraction
-        );
+            (constraint.first() > patchConstraints[pointI].first())
+         || (magSqr(attraction) < magSqr(patchAttraction[pointI]))
+        )
+        {
+            patchAttraction[pointI] = attraction;
+            patchConstraints[pointI] = constraint;
+
+            const point& pt = pp.localPoints()[pointI];
+
+            if (patchConstraints[pointI].first() == 2 && feStr.valid())
+            {
+                meshTools::writeOBJ(feStr(), pt);
+                feVertI++;
+                meshTools::writeOBJ(feStr(), pt+patchAttraction[pointI]);
+                feVertI++;
+                feStr() << "l " << feVertI-1 << ' ' << feVertI << nl;
+            }
+            else if (patchConstraints[pointI].first() == 3 && fpStr.valid())
+            {
+                meshTools::writeOBJ(fpStr(), pt);
+                fpVertI++;
+                meshTools::writeOBJ(fpStr(), pt+patchAttraction[pointI]);
+                fpVertI++;
+                fpStr() << "l " << fpVertI-1 << ' ' << fpVertI << nl;
+            }
+        }
     }
 }
 
 
-// Finds nearest feature (within snapDist) for all points of pp.
-void Foam::autoSnapDriver::determineAllFeatures
+Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
 (
-    const label iter,
-    const scalar featureCos,
-
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
+    const label pointI,
+    const point& estimatedPt,
 
-    const List<List<point> >& pointFaceSurfNormals,
-    const List<List<point> >& pointFaceDisp,
-    const List<List<point> >& pointFaceCentres,
+    label& featI,
+    List<List<DynamicList<point> > >& edgeAttractors,
+    List<List<DynamicList<pointConstraint> > >& edgeConstraints,
+    vectorField& patchAttraction,
+    List<pointConstraint>& patchConstraints
+) const
+{
+    const refinementFeatures& features = meshRefiner_.features();
 
+    labelList nearEdgeFeat;
+    List<pointIndexHit> nearEdgeInfo;
+    features.findNearestEdge
+    (
+        pointField(1, estimatedPt),
+        scalarField(1, sqr(snapDist[pointI])),
+        nearEdgeFeat,
+        nearEdgeInfo
+    );
+
+    const pointIndexHit& nearInfo = nearEdgeInfo[0];
+    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;
+        const edge e = features[featI].edges()[nearInfo.index()];
+        vector eVec = e.vec(features[featI].points());
+        eVec /= mag(eVec)+VSMALL;
+        c.first() = 2;
+        c.second() = eVec;
+        edgeConstraints[featI][nearInfo.index()].append(c);
+
+        // Store for later use
+        patchAttraction[pointI] =
+            nearInfo.hitPoint()-pp.localPoints()[pointI];
+        patchConstraints[pointI] = c;
+    }
+    return nearInfo;
+}
+Foam::labelPair Foam::autoSnapDriver::findNearFeaturePoint
+(
+    const indirectPrimitivePatch& pp,
+    const scalarField& snapDist,
+    const label pointI,
+    const point& estimatedPt,
+
+    // Feature-point to pp point
     List<labelList>& pointAttractor,
     List<List<pointConstraint> >& pointConstraints,
     // Feature-edge to pp point
@@ -735,134 +1071,94 @@ void Foam::autoSnapDriver::determineAllFeatures
     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 nearFeat;
+    labelList nearIndex;
+    features.findNearestPoint
+    (
+        pointField(1, estimatedPt),
+        scalarField(1, sqr(snapDist[pointI])),
+        nearFeat,
+        nearIndex
+    );
 
-    labelList nearPointFeat;
-    labelList nearPointIndex;
-    {
-        scalarField snapDistSqr(sqr(snapDist));
-        features.findNearestEdge
-        (
-            pp.localPoints(),
-            snapDistSqr,
-            nearEdgeFeat,
-            nearEdgeInfo
-        );
+    label featI = nearFeat[0];
+    label featPointI = -1;
 
-        // Look at near feature points
-        features.findNearestPoint
-        (
-            pp.localPoints(),
-            snapDistSqr,
-            nearPointFeat,
-            nearPointIndex
-        );
-    }
-
-    forAll(pp.localPoints(), pointI)
+    if (featI != -1)
     {
         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())
+        const treeDataPoint& shapes =
+            features.pointTrees()[featI].shapes();
+        featPointI = shapes.pointLabels()[nearIndex[0]];
+        const point& featPt = shapes.points()[featPointI];
+        scalar distSqr = magSqr(featPt-pt);
+
+        // Check if already attracted
+        label oldPointI = pointAttractor[featI][featPointI];
+
+        if (oldPointI != -1)
         {
-            // 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())
+            // Check distance
+            if (distSqr >= magSqr(featPt-pp.localPoints()[oldPointI]))
             {
-                meshTools::writeOBJ(featureEdgeStr(), pt);
-                featureEdgeVertI++;
-                meshTools::writeOBJ(featureEdgeStr(), nearEdge.hitPoint());
-                featureEdgeVertI++;
-                featureEdgeStr()
-                    << "l " << featureEdgeVertI-1 << ' '
-                    << featureEdgeVertI << nl;
+                // oldPointI nearest. Keep.
+                featI = -1;
+                featPointI = -1;
+            }
+            else
+            {
+                // Current pointI nearer.
+                pointAttractor[featI][featPointI] = pointI;
+                pointConstraints[featI][featPointI].first() = 3;
+                pointConstraints[featI][featPointI].second() = vector::zero;
+
+                // Store for later use
+                patchAttraction[pointI] = featPt-pt;
+                patchConstraints[pointI] =
+                    pointConstraints[featI][featPointI];
+
+                // Reset oldPointI to nearest on feature edge
+                patchAttraction[oldPointI] = vector::zero;
+                patchConstraints[oldPointI] = pointConstraint();
+
+                label edgeFeatI;
+                const pointIndexHit nearInfo = findNearFeatureEdge
+                (
+                    pp,
+                    snapDist,
+                    oldPointI,
+                    pp.localPoints()[oldPointI],
+
+                    edgeFeatI,
+                    edgeAttractors,
+                    edgeConstraints,
+                    patchAttraction,
+                    patchConstraints
+                );
             }
         }
-
-
-        // Mark point on the nearest feature point.
-        if (nearPointFeat[pointI] != -1)
+        else
         {
-            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];
-
+            // Current pointI nearer.
             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;
-            }
+            // Store for later use
+            patchAttraction[pointI] = featPt-pt;
+            patchConstraints[pointI] = pointConstraints[featI][featPointI];
         }
     }
+
+    return labelPair(featI, featPointI);
 }
 
 
+// Determines for every pp point - that is on multiple faces that form
+// a feature - the nearest feature edge/point.
 void Foam::autoSnapDriver::determineFeatures
 (
     const label iter,
@@ -871,12 +1167,10 @@ void Foam::autoSnapDriver::determineFeatures
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
 
-    //const vectorField& faceSurfaceNormal,
-    //const vectorField& faceDisp,
-    //const vectorField& faceRotation,
     const List<List<point> >& pointFaceSurfNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
+    const labelListList& pointFacePatchID,
 
     // Feature-point to pp point
     List<labelList>& pointAttractor,
@@ -906,7 +1200,7 @@ void Foam::autoSnapDriver::determineFeatures
               / "featureEdge_" + name(iter) + ".obj"
             )
         );
-        Pout<< nl << "Dumping feature-edge sampling to "
+        Pout<< "Dumping feature-edge sampling to "
             << featureEdgeStr().name() << endl;
 
         missedEdgeStr.reset
@@ -917,7 +1211,7 @@ void Foam::autoSnapDriver::determineFeatures
               / "missedFeatureEdge_" + name(iter) + ".obj"
             )
         );
-        Pout<< nl << "Dumping feature-edges that are too far away to "
+        Pout<< "Dumping feature-edges that are too far away to "
             << missedEdgeStr().name() << endl;
 
         featurePointStr.reset
@@ -928,7 +1222,7 @@ void Foam::autoSnapDriver::determineFeatures
               / "featurePoint_" + name(iter) + ".obj"
             )
         );
-        Pout<< nl << "Dumping feature-point sampling to "
+        Pout<< "Dumping feature-point sampling to "
             << featurePointStr().name() << endl;
     }
 
@@ -936,90 +1230,137 @@ void Foam::autoSnapDriver::determineFeatures
 
     forAll(pp.localPoints(), pointI)
     {
-        if (patchConstraints[pointI].first() == 0)
+        const point& pt = pp.localPoints()[pointI];
+
+        vector attraction = vector::zero;
+        pointConstraint constraint;
+
+        featureAttractionUsingReconstruction
+        (
+            iter,
+            featureCos,
+
+            pp,
+            snapDist,
+            pointI,
+
+            pointFaceSurfNormals,
+            pointFaceDisp,
+            pointFaceCentres,
+            pointFacePatchID,
+
+            attraction,
+            constraint
+        );
+
+        if
+        (
+            (constraint.first() > patchConstraints[pointI].first())
+         || (magSqr(attraction) < magSqr(patchAttraction[pointI]))
+        )
         {
-            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,
-                pointFaceSurfNormals,
-                pointFaceDisp,
-                pointFaceCentres,
-
-                surfacePoints,
-                surfaceNormals,
-                surfaceCount
-            );
+            patchAttraction[pointI] = attraction;
+            patchConstraints[pointI] = constraint;
 
             // Check the number of directions
-            if (surfaceNormals.size() == 2)
+            if (patchConstraints[pointI].first() == 1)
             {
-                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]))
+                // Flat surface. Check for different patchIDs
+                pointIndexHit multiPatchPt
+                (
+                    findMultiPatchPoint
+                    (
+                        pt,
+                        pointFacePatchID[pointI],
+                        pointFaceCentres[pointI]
+                    )
+                );
+                if (multiPatchPt.hit())
                 {
-                    d *= Foam::sqrt(sqr(snapDist[pointI])/magSqr(d));
+                    // Behave like when having two surface normals so
+                    // attract to nearest feature edge (with a guess for
+                    // the multipatch point as starting point)
+                    label featI = -1;
+                    pointIndexHit nearInfo = findNearFeatureEdge
+                    (
+                        pp,
+                        snapDist,
+                        pointI,
+                        multiPatchPt.hitPoint(),        //estimatedPt
+
+                        featI,
+                        edgeAttractors,
+                        edgeConstraints,
+
+                        patchAttraction,
+                        patchConstraints
+                    );
+
+                    if (nearInfo.hit())
+                    {
+                        // 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 (patchConstraints[pointI].first() == 2)
+            {
                 // 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);
+                const point estimatedPt(pt + patchAttraction[pointI]);
 
-                labelList nearEdgeFeat;
-                List<pointIndexHit> nearEdgeInfo;
-                features.findNearestEdge
+                // Determine nearest point on feature edge. Store constraint
+                // (calculated from feature edge, alternative would be to
+                //  use constraint calculated from both surfaceNormals)
+                label featI = -1;
+                pointIndexHit nearInfo = findNearFeatureEdge
                 (
-                    pointField(1, estimatedPt),
-                    scalarField(1, sqr(snapDist[pointI])),
-                    nearEdgeFeat,
-                    nearEdgeInfo
-                );
+                    pp,
+                    snapDist,
+                    pointI,
+                    estimatedPt,
 
-                const pointIndexHit& nearInfo = nearEdgeInfo[0];
-                label featI = nearEdgeFeat[0];
+                    featI,
+                    edgeAttractors,
+                    edgeConstraints,
+
+                    patchAttraction,
+                    patchConstraints
+                );
 
                 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())
                     {
@@ -1054,59 +1395,39 @@ void Foam::autoSnapDriver::determineFeatures
                     }
                 }
             }
-            else if (surfaceNormals.size() == 3)
+            else if (patchConstraints[pointI].first() == 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);
+                const point estimatedPt(pt + patchAttraction[pointI]);
 
-                labelList nearPointFeat;
-                labelList nearPointIndex;
-                features.findNearestPoint
+                labelPair nearInfo = findNearFeaturePoint
                 (
-                    pointField(1, estimatedPt),
-                    scalarField(1, sqr(snapDist[pointI])),
-                    nearPointFeat,
-                    nearPointIndex
+                    pp,
+                    snapDist,
+                    pointI,
+                    estimatedPt,
+
+                    // Feature-point to pp point
+                    pointAttractor,
+                    pointConstraints,
+                    // Feature-edge to pp point
+                    edgeAttractors,
+                    edgeConstraints,
+                    // pp point to nearest feature
+                    patchAttraction,
+                    patchConstraints
                 );
 
-
-                if (nearPointIndex[0] != -1)
+                if (nearInfo.first() != -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())
                     {
+                        const treeDataPoint& shapes =
+                            features.pointTrees()[nearInfo.first()].shapes();
+                        const point& featPt =
+                            shapes.points()[nearInfo.second()];
+
                         meshTools::writeOBJ(featurePointStr(), pt);
                         featurePointVertI++;
                         meshTools::writeOBJ(featurePointStr(), featPt);
@@ -1130,9 +1451,6 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
 
-    //const vectorField& faceSurfaceNormal,
-    //const vectorField& faceDisp,
-    //const vectorField& faceRotation,
     const List<List<point> >& pointFaceSurfNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
@@ -1145,12 +1463,15 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     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());
+    List<List<DynamicList<pointConstraint> > > edgeConstraints
+    (
+        features.size()
+    );
     forAll(features, featI)
     {
         label nFeatEdges = features[featI].edges().size();
@@ -1175,7 +1496,6 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     List<pointConstraint> allPatchConstraints(pp.nPoints());
 
     determineFeatures
-    //determineAllFeatures
     (
         iter,
         featureCos,
@@ -1186,6 +1506,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
         pointFaceSurfNormals,
         pointFaceDisp,
         pointFaceCentres,
+        pointFacePatchID,
 
         // Feature-point to pp point
         pointAttractor,
@@ -1199,6 +1520,263 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     );
 
 
+
+    // Baffle handling
+    // ~~~~~~~~~~~~~~~
+    // Override pointAttractor, edgeAttractor, allPatchAttration etc. to
+    // implement 'baffle' handling.
+    // Baffle: the mesh pp point originates from a loose standing
+    // baffle.
+    // Sampling the surface with the surrounding face-centres only picks up
+    // a single triangle normal so above determineFeatures will not have
+    // detected anything. So explicitly pick up feature edges on the pp
+    // (after duplicating points & smoothing so will already have been
+    // expanded) and match these to the features.
+    {
+        const fvMesh& mesh = meshRefiner_.mesh();
+
+        // Calculate edge-faces
+        List<List<point> > edgeFaceNormals(pp.nEdges());
+
+        // Fill local data
+        forAll(pp.edgeFaces(), edgeI)
+        {
+            const labelList& eFaces = pp.edgeFaces()[edgeI];
+            List<point>& eFc = edgeFaceNormals[edgeI];
+            eFc.setSize(eFaces.size());
+            forAll(eFaces, i)
+            {
+                label faceI = eFaces[i];
+                eFc[i] = pp.faceNormals()[faceI];
+            }
+        }
+
+        {
+            // Precalculate mesh edges for pp.edges.
+            const labelList meshEdges
+            (
+                pp.meshEdges(mesh.edges(), mesh.pointEdges())
+            );
+            syncTools::syncEdgeList
+            (
+                mesh,
+                meshEdges,
+                edgeFaceNormals,
+                listPlusEqOp<point>(),
+                List<point>(),
+                listTransform()
+            );
+        }
+
+        // Detect baffle edges. Assume initial mesh will have 0,90 or 180
+        // (baffle) degree angles so smoothing should make 0,90
+        // to be less than 90.
+        const scalar baffleFeatureCos = Foam::cos(degToRad(91));
+
+
+        autoPtr<OFstream> baffleEdgeStr;
+        label baffleEdgeVertI = 0;
+        if (debug&meshRefinement::OBJINTERSECTIONS)
+        {
+            baffleEdgeStr.reset
+            (
+                new OFstream
+                (
+                    meshRefiner_.mesh().time().path()
+                  / "baffleEdge_" + name(iter) + ".obj"
+                )
+            );
+            Info<< nl << "Dumping baffle-edges to "
+                << baffleEdgeStr().name() << endl;
+        }
+
+
+        // Is edge on baffle
+        PackedBoolList isBaffleEdge(pp.nEdges());
+        // Is point on
+        //  0 : baffle-edge (0)
+        //  1 : baffle-feature-point (1)
+        // -1 : rest
+        labelList pointStatus(pp.nPoints(), -1);
+
+        forAll(edgeFaceNormals, edgeI)
+        {
+            const List<point>& efn = edgeFaceNormals[edgeI];
+
+            if (efn.size() == 2 && (efn[0]&efn[1]) < baffleFeatureCos)
+            {
+                isBaffleEdge[edgeI] = true;
+                const edge& e = pp.edges()[edgeI];
+                pointStatus[e[0]] = 0;
+                pointStatus[e[1]] = 0;
+
+                if (baffleEdgeStr.valid())
+                {
+                    const point& p0 = pp.localPoints()[e[0]];
+                    const point& p1 = pp.localPoints()[e[1]];
+                    meshTools::writeOBJ(baffleEdgeStr(), p0);
+                    baffleEdgeVertI++;
+                    meshTools::writeOBJ(baffleEdgeStr(), p1);
+                    baffleEdgeVertI++;
+                    baffleEdgeStr() << "l " << baffleEdgeVertI-1
+                        << ' ' << baffleEdgeVertI << nl;
+                }
+            }
+        }
+
+        forAll(pp.pointEdges(), pointI)
+        {
+            if
+            (
+                isFeaturePoint
+                (
+                    featureCos,
+                    pp,
+                    isBaffleEdge,
+                    pointI
+                )
+            )
+            {
+                //Pout<< "Detected feature point:" << pp.localPoints()[pointI]
+                //    << endl;
+                //-TEMPORARILY DISABLED:
+                //pointStatus[pointI] = 1;
+            }
+        }
+
+        forAll(pointStatus, pointI)
+        {
+            const point& pt = pp.localPoints()[pointI];
+
+            if (pointStatus[pointI] == 0)   // baffle edge
+            {
+                label featI;
+                const pointIndexHit nearInfo = findNearFeatureEdge
+                (
+                    pp,
+                    snapDist,
+                    pointI,
+                    pt,
+
+                    featI,
+                    edgeAttractors,
+                    edgeConstraints,
+                    allPatchAttraction,
+                    allPatchConstraints
+                );
+
+                if (!nearInfo.hit())
+                {
+                    //Pout<< "*** Failed to find close edge to point " << pt
+                    //    << endl;
+                }
+            }
+            else if (pointStatus[pointI] == 1)   // baffle point
+            {
+                labelList nearFeat;
+                labelList nearIndex;
+                features.findNearestPoint
+                (
+                    pointField(1, pt),
+                    scalarField(1, sqr(snapDist[pointI])),
+                    nearFeat,
+                    nearIndex
+                );
+
+                label featI = nearFeat[0];
+
+                if (featI != -1)
+                {
+                    const treeDataPoint& shapes =
+                        features.pointTrees()[featI].shapes();
+                    label featPointI = shapes.pointLabels()[nearIndex[0]];
+                    const point& featPt = shapes.points()[featPointI];
+                    scalar distSqr = magSqr(featPt-pt);
+
+                    // Check if already attracted
+                    label oldPointI = pointAttractor[featI][featPointI];
+
+                    if
+                    (
+                        oldPointI == -1
+                     || (
+                            distSqr
+                          < magSqr(featPt-pp.localPoints()[oldPointI])
+                        )
+                    )
+                    {
+                        pointAttractor[featI][featPointI] = pointI;
+                        pointConstraints[featI][featPointI].first() = 3;
+                        pointConstraints[featI][featPointI].second() =
+                            vector::zero;
+
+                        // Store for later use
+                        allPatchAttraction[pointI] = featPt-pt;
+                        allPatchConstraints[pointI] =
+                            pointConstraints[featI][featPointI];
+
+                        if (oldPointI != -1)
+                        {
+                            // The current point is closer so wins. Reset
+                            // the old point to attract to nearest edge
+                            // instead.
+                            label edgeFeatI;
+                            const pointIndexHit nearInfo = findNearFeatureEdge
+                            (
+                                pp,
+                                snapDist,
+                                oldPointI,
+                                pp.localPoints()[oldPointI],
+
+                                edgeFeatI,
+                                edgeAttractors,
+                                edgeConstraints,
+                                allPatchAttraction,
+                                allPatchConstraints
+                            );
+                        }
+                    }
+                    else
+                    {
+                        // Make it fall through to check below
+                        featI = -1;
+                    }
+                }
+
+                // Not found a feature point or another point is already
+                // closer to that feature
+                if (featI == -1)
+                {
+                    //Pout<< "*** Falling back to finding nearest feature edge"
+                    //    << " for baffle-feature-point " << pt
+                    //    << endl;
+
+                    label featI;
+                    const pointIndexHit nearInfo = findNearFeatureEdge
+                    (
+                        pp,
+                        snapDist,
+                        pointI,
+                        pt,                     // starting point
+
+                        featI,
+                        edgeAttractors,
+                        edgeConstraints,
+                        allPatchAttraction,
+                        allPatchConstraints
+                    );
+                }
+            }
+        }
+    }
+
+
+    //
+    // Reverse lookup
+    // ~~~~~~~~~~~~~~
+    //
+
+
     // Find nearest mesh point to feature edge
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Reverse lookup : go through all edgeAttractors and find the
@@ -1244,9 +1822,12 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
             forAll(attr, i)
             {
                 // Find nearest pp point
-
                 const point& featPt = attr[i];
-                pointIndexHit nearInfo = ppTree.findNearest(featPt, sqr(GREAT));
+                pointIndexHit nearInfo = ppTree.findNearest
+                (
+                    featPt,
+                    sqr(GREAT)
+                );
 
                 if (nearInfo.hit())
                 {
@@ -1257,7 +1838,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                     // override it only if nearer.
                     if
                     (
-                        patchConstraints[pointI].first() == 0
+                        patchConstraints[pointI].first() <= 1
                      || magSqr(attraction) < magSqr(patchAttraction[pointI])
                     )
                     {
@@ -1269,7 +1850,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                 {
                     WarningIn
                     (
-                        "autoSnapDriver::featureAttractionUsingFeatureEdges(..)"
+                        "autoSnapDriver::featureAttractionUsingFeatureEdges"
+                        "(..)"
                     )   << "Did not find pp point near " << featPt
                         << endl;
                 }
@@ -1302,7 +1884,11 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                 ];
 
                 // Find nearest pp point
-                pointIndexHit nearInfo = ppTree.findNearest(featPt, sqr(GREAT));
+                pointIndexHit nearInfo = ppTree.findNearest
+                (
+                    featPt,
+                    sqr(GREAT)
+                );
 
                 if (nearInfo.hit())
                 {
@@ -1333,7 +1919,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                         )
                         {
                             patchAttraction[pointI] = attraction;
-                            patchConstraints[pointI] = pointConstr[featPointI];
+                            patchConstraints[pointI] =
+                                pointConstr[featPointI];
                         }
                     }
                 }
@@ -1344,7 +1931,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
 
 
     //MEJ: any faces that have multi-patch points only keep the multi-patch
-    //     points.
+    //     points. The other points on the face will be dragged along
+    //     (hopefully)
     {
         autoPtr<OFstream> multiPatchStr;
         if (debug&meshRefinement::OBJINTERSECTIONS)
@@ -1357,7 +1945,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                   / "multiPatch_" + name(iter) + ".obj"
                 )
             );
-            Pout<< nl << "Dumping removed constraints due to same-face"
+            Pout<< "Dumping removed constraints due to same-face"
                 << " multi-patch points to "
                 << multiPatchStr().name() << endl;
         }
@@ -1368,16 +1956,13 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
 
         forAll(pointFacePatchID, pointI)
         {
-            const labelList& patches = pointFacePatchID[pointI];
-            label patch0 = patches[0];
-            for (label i = 1; i < patches.size(); i++)
-            {
-                if (patch0 != patches[i])
-                {
-                    isMultiPatchPoint[pointI] = 1u;
-                    break;
-                }
-            }
+            pointIndexHit multiPatchPt = findMultiPatchPoint
+            (
+                pp.localPoints()[pointI],
+                pointFacePatchID[pointI],
+                pointFaceCentres[pointI]
+            );
+            isMultiPatchPoint[pointI] = multiPatchPt.hit();
         }
 
         // 2. Make sure multi-patch points are also attracted
@@ -1387,8 +1972,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
             {
                 if
                 (
-                    patchConstraints[pointI].first() == 0
-                 && allPatchConstraints[pointI].first() > 0
+                    patchConstraints[pointI].first() <= 1
+                 && allPatchConstraints[pointI].first() > 1
                 )
                 {
                     patchAttraction[pointI] = allPatchAttraction[pointI];
@@ -1422,7 +2007,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                 if
                 (
                     isMultiPatchPoint[pointI]
-                 && patchConstraints[pointI].first() != 0
+                 && patchConstraints[pointI].first() > 1
                 )
                 {
                     nMultiPatchPoints++;
@@ -1437,7 +2022,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                     if
                     (
                        !isMultiPatchPoint[pointI]
-                     && patchConstraints[pointI].first() != 0
+                     && patchConstraints[pointI].first() > 1
                     )
                     {
                         //Pout<< "Knocking out constraint"
@@ -1476,7 +2061,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
           / "edgeAttractors_" + name(iter) + ".obj"
         );
         label featureEdgeVertI = 0;
-        Pout<< nl << "Dumping feature-edge attraction to "
+        Pout<< "Dumping feature-edge attraction to "
             << featureEdgeStr.name() << endl;
 
         OFstream featurePointStr
@@ -1485,7 +2070,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
           / "pointAttractors_" + name(iter) + ".obj"
         );
         label featurePointVertI = 0;
-        Pout<< nl << "Dumping feature-point attraction to "
+        Pout<< "Dumping feature-point attraction to "
             << featurePointStr.name() << endl;
 
         forAll(patchConstraints, pointI)
@@ -1496,7 +2081,11 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
             {
                 meshTools::writeOBJ(featureEdgeStr, pt);
                 featureEdgeVertI++;
-                meshTools::writeOBJ(featureEdgeStr, pt+patchAttraction[pointI]);
+                meshTools::writeOBJ
+                (
+                    featureEdgeStr,
+                    pt+patchAttraction[pointI]
+                );
                 featureEdgeVertI++;
                 featureEdgeStr << "l " << featureEdgeVertI-1
                     << ' ' << featureEdgeVertI << nl;
@@ -1547,10 +2136,11 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                     const edge& e = pp.edges()[pEdges[pEdgeI]];
                     label nbrPointI = e.otherVertex(pointI);
 
-                    if (patchConstraints[nbrPointI].first() != 0)
+                    if (patchConstraints[nbrPointI].first() > 1)
                     {
                         const point& nbrPt = pp.localPoints()[nbrPointI];
-                        const point featPt = nbrPt + patchAttraction[nbrPointI];
+                        const point featPt =
+                            nbrPt + patchAttraction[nbrPointI];
                         const scalar cosAngle = (featVec & (featPt-pt));
 
                         if (cosAngle > 0)
@@ -1582,8 +2172,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
 
                         if
                         (
-                            patchConstraints[nbrPointI].first() == 0
-                         && allPatchConstraints[nbrPointI].first() != 0
+                            patchConstraints[nbrPointI].first() <= 1
+                         && allPatchConstraints[nbrPointI].first() > 1
                         )
                         {
                             const vector& nbrFeatVec =
@@ -1602,7 +2192,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                                 const point featPt =
                                     pp.localPoints()[nbrPointI]
                                   + allPatchAttraction[nbrPointI];
-                                const scalar cosAngle = (featVec & (featPt-pt));
+                                const scalar cosAngle =
+                                    (featVec & (featPt-pt));
 
                                 if (cosAngle > 0)
                                 {
@@ -1628,7 +2219,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                     {
                         // Use reconstructed-feature attraction. Use only
                         // part of it since not sure...
-                        //const point& bestPt = pp.localPoints()[bestPosPointI];
+                        //const point& bestPt =
+                        //    pp.localPoints()[bestPosPointI];
                         //Pout<< "**Overriding point " << bestPt
                         //    << " on reconstructed feature edge at "
                         //    << allPatchAttraction[bestPosPointI]+bestPt
@@ -1644,7 +2236,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                     {
                         // Use reconstructed-feature attraction. Use only
                         // part of it since not sure...
-                        //const point& bestPt = pp.localPoints()[bestNegPointI];
+                        //const point& bestPt =
+                        //    pp.localPoints()[bestNegPointI];
                         //Pout<< "**Overriding point " << bestPt
                         //    << " on reconstructed feature edge at "
                         //    << allPatchAttraction[bestNegPointI]+bestPt
@@ -1708,9 +2301,11 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                     // 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 point pt0 =
+                        pp.localPoints()[i0]+patchAttraction[i0];
                     const label i1 = f[nextAttract];
-                    const point pt1 = pp.localPoints()[i1]+patchAttraction[i1];
+                    const point pt1 =
+                        pp.localPoints()[i1]+patchAttraction[i1];
                     const point mid = 0.5*(pt0+pt1);
 
 
@@ -1728,13 +2323,14 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
 
                     if (cosAngle > featureCos)
                     {
-                        // Add the nearest of the other two points as attractor
+                        // 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)
+                            if (patchConstraints[pointI].first() <= 1)
                             {
                                 const point& pt = pp.localPoints()[pointI];
                                 scalar distSqr = magSqr(mid-pt);
@@ -1822,14 +2418,16 @@ void Foam::autoSnapDriver::preventFaceSqueeze
         forAll(f, fp)
         {
             label pointI = f[fp];
-            if (patchConstraints[pointI].first() != 0)
+            const point& pt = pp.localPoints()[pointI];
+
+            if (patchConstraints[pointI].first() > 1)
             {
-                points[fp] = pp.localPoints()[pointI] + patchAttraction[pointI];
+                points[fp] = pt + patchAttraction[pointI];
                 nConstraints++;
             }
             else
             {
-                points[fp] = pp.localPoints()[pointI];
+                points[fp] = pt;
             }
         }
 
@@ -1865,6 +2463,7 @@ void Foam::autoSnapDriver::preventFaceSqueeze
 
 Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 (
+    const snapParameters& snapParams,
     const label iter,
     const scalar featureCos,
     const scalar featureAttract,
@@ -1873,8 +2472,14 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     motionSmoother& meshMover
 ) const
 {
-    Info<< "Calculating patchDisplacement as distance to nearest surface"
-        << " point ..." << endl;
+    const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
+    const Switch explicitFeatureAttraction = snapParams.explicitFeatureSnap();
+
+    Info<< "Overriding displacement on features :" << nl
+        << "   implicit features : " << implicitFeatureAttraction << nl
+        << "   explicit features : " << explicitFeatureAttraction << nl
+        << endl;
+
 
     const indirectPrimitivePatch& pp = meshMover.patch();
     const pointField& localPoints = pp.localPoints();
@@ -1887,6 +2492,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     vectorField faceDisp(pp.size(), vector::zero);
     // normal of surface at point on surface
     vectorField faceSurfaceNormal(pp.size(), vector::zero);
+    labelList faceSurfaceRegion(pp.size(), -1);
     vectorField faceRotation(pp.size(), vector::zero);
 
     calcNearestFace
@@ -1895,11 +2501,24 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         pp,
         faceDisp,
         faceSurfaceNormal,
+        faceSurfaceRegion,
         faceRotation
     );
 
-    // Start off with nearest point on surface.
-    vectorField patchDisp = nearestDisp;
+
+    //// Displacement and orientation per pp point
+    //// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    //vectorField pointDisp(pp.nPoints(), vector::zero);
+    //vectorField pointSurfaceNormal(pp.nPoints(), vector::zero);
+    //vectorField pointRotation(pp.nPoints(), vector::zero);
+    //calcNearest
+    //(
+    //    iter,
+    //    pp,
+    //    pointDisp,
+    //    pointSurfaceNormal,
+    //    pointRotation
+    //);
 
 
     // Collect (possibly remote) per point data of all surrounding faces
@@ -1910,10 +2529,12 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     // - faceCentres&faceNormal
 
     // For now just get all surrounding face data. Expensive - should just
-    // store and sync data on coupled points only (see e.g PatchToolsNormals.C)
+    // store and sync data on coupled points only
+    // (see e.g PatchToolsNormals.C)
 
     List<List<point> > pointFaceSurfNormals(pp.nPoints());
     List<List<point> > pointFaceDisp(pp.nPoints());
+    //List<List<point> > pointFaceRotation(pp.nPoints());
     List<List<point> > pointFaceCentres(pp.nPoints());
     List<labelList>    pointFacePatchID(pp.nPoints());
 
@@ -1925,6 +2546,8 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         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());
         labelList& pFid = pointFacePatchID[pointI];
@@ -1935,6 +2558,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
             label faceI = pFaces[i];
             pNormals[i] = faceSurfaceNormal[faceI];
             pDisp[i] = faceDisp[faceI];
+            //pRot[i] = faceRotation[faceI];
             pFc[i] = pp.faceCentres()[faceI];
             label meshFaceI = pp.addressing()[faceI];
             pFid[i] = mesh.boundaryMesh().whichPatch(meshFaceI);
@@ -1959,6 +2583,15 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         List<point>(),
         listTransform()
     );
+    //syncTools::syncPointList
+    //(
+    //    mesh,
+    //    pp.meshPoints(),
+    //    pointFaceRotation,
+    //    listPlusEqOp<point>(),
+    //    List<point>(),
+    //    listTransform()
+    //);
     syncTools::syncPointList
     (
         mesh,
@@ -1978,6 +2611,9 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     );
 
 
+    // Start off with nearest point on surface
+    vectorField patchDisp = nearestDisp;
+
 
     // Main calculation
     // ~~~~~~~~~~~~~~~~
@@ -1989,23 +2625,56 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     vectorField patchAttraction(localPoints.size(), vector::zero);
     // Constraints at feature
     List<pointConstraint> patchConstraints(localPoints.size());
-    //featureAttractionUsingReconstruction
-    featureAttractionUsingFeatureEdges
-    (
-        iter,
-        featureCos,
 
-        pp,
-        snapDist,
 
-        pointFaceSurfNormals,
-        pointFaceDisp,
-        pointFaceCentres,
-        pointFacePatchID,
+    if (implicitFeatureAttraction)
+    {
+        // Sample faces around each point and see if nearest surface normal
+        // differs. Reconstruct a feature edge/point if possible and snap to
+        // it.
+        featureAttractionUsingReconstruction
+        (
+            iter,
+            featureCos,
 
-        patchAttraction,
-        patchConstraints
-    );
+            pp,
+            snapDist,
+
+            pointFaceSurfNormals,
+            pointFaceDisp,
+            pointFaceCentres,
+            pointFacePatchID,
+
+            patchAttraction,
+            patchConstraints
+        );
+    }
+
+    if (explicitFeatureAttraction)
+    {
+        // Sample faces around each point and see if nearest surface normal
+        // differs. For those find the nearest real feature edge/point and
+        // store the correspondence. Then loop over feature edge/point
+        // and attract those nearest mesh point. (the first phase just is
+        // a subsetting of candidate points, the second makes sure that only
+        // one mesh point gets attracted per feature)
+        featureAttractionUsingFeatureEdges
+        (
+            iter,
+            featureCos,
+
+            pp,
+            snapDist,
+
+            pointFaceSurfNormals,
+            pointFaceDisp,
+            pointFaceCentres,
+            pointFacePatchID,
+
+            patchAttraction,
+            patchConstraints
+        );
+    }
 
     preventFaceSqueeze
     (
@@ -2036,14 +2705,13 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     // - patchAttraction    : direct attraction to features
     // - patchConstraints   : type of features
 
-    // Use any combination of patchDisp and direct feature
-    // attraction.
+    // Use any combination of patchDisp and direct feature attraction.
 
 
     // Mix with direct feature attraction
     forAll(patchConstraints, pointI)
     {
-        if (patchConstraints[pointI].first() != 0)
+        if (patchConstraints[pointI].first() > 1)
         {
             patchDisp[pointI] =
                 (1.0-featureAttract)*patchDisp[pointI]
@@ -2051,14 +2719,6 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         }
     }
 
-    //dumpMove
-    //(
-    //    mesh.time().path()
-    //  / "linearPatchDisp_" + name(iter) + ".obj",
-    //    pp.localPoints(),
-    //    pp.localPoints() + patchDisp
-    //);
-
 
 
     // Count
@@ -2092,13 +2752,13 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
             << " attraction to :" << nl
             << "    feature point   : " << nPoint << nl
             << "    feature edge    : " << nEdge << nl
-            << "    nearest surface : " << nTotPoints-nPoint-nEdge
-            << " (rest)" << nl
+            << "    nearest surface : " << nPlanar << nl
+            << "    rest            : " << nTotPoints-nPoint-nEdge-nPlanar
+            << 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
@@ -2106,6 +2766,10 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     // - finish with only tangential component smoothed.
     // Note: tangential is most
     // likely to come purely from face-centre snapping, not face rotation.
+    // Note: could use the constraints here (constraintTransformation())
+    //       but this is not necessarily accurate and we're smoothing to
+    //       get out of problems.
+
     if (featureAttract < 1-0.001)
     {
         // 1. Smoothed all displacement
@@ -2148,8 +2812,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     }
 
 
-
-    const scalar relax = featureAttract; //1.0;
+    const scalar relax = featureAttract;
     patchDisp *= relax;
 
 
@@ -2161,7 +2824,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         pp.meshPoints(),
         patchDisp,
         minMagSqrEqOp<point>(),         // combine op
-        vector(GREAT, GREAT, GREAT)     // null value (note: cannot use VGREAT)
+        vector(GREAT, GREAT, GREAT)     // null value (note: cant use VGREAT)
     );
 
     return patchDisp;
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C
index 54231f8bd1..b9d6f3b0c6 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) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,7 +34,9 @@ Foam::snapParameters::snapParameters(const dictionary& dict)
     snapTol_(readScalar(dict.lookup("tolerance"))),
     nSmoothDispl_(readLabel(dict.lookup("nSolveIter"))),
     nSnap_(readLabel(dict.lookup("nRelaxIter"))),
-    nFeatureSnap_(dict.lookupOrDefault("nFeatureSnapIter", -1))
+    nFeatureSnap_(dict.lookupOrDefault("nFeatureSnapIter", -1)),
+    explicitFeatureSnap_(dict.lookupOrDefault("explicitFeatureSnap", true)),
+    implicitFeatureSnap_(dict.lookupOrDefault("implicitFeatureSnap", false))
 {}
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H
index 3895c58aad..8103cd4082 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) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,6 +37,7 @@ SourceFiles
 
 #include "dictionary.H"
 #include "scalar.H"
+#include "Switch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -63,6 +64,10 @@ class snapParameters
 
         const label nFeatureSnap_;
 
+        const Switch explicitFeatureSnap_;
+
+        const Switch implicitFeatureSnap_;
+
 
     // Private Member Functions
 
@@ -119,6 +124,16 @@ public:
                 return nFeatureSnap_;
             }
 
+            Switch explicitFeatureSnap() const
+            {
+                return explicitFeatureSnap_;
+            }
+
+            Switch implicitFeatureSnap() const
+            {
+                return implicitFeatureSnap_;
+            }
+
 };
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 4deb6f3eed..a3b521eb12 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -326,6 +326,16 @@ private:
                 const labelList& globalToPatch
             ) const;
 
+            //- Geometric test on see whether face needs to be baffled:
+            //  is face boundary face and perpendicular to surface normal?
+            bool validBaffleTopology
+            (
+                const label faceI,
+                const vector& n1,
+                const vector& n2,
+                const vector& testDir
+            ) const;
+
             //- Determine patches for baffles
             void getBafflePatches
             (
@@ -414,7 +424,7 @@ private:
 
             //- Extract those baffles (duplicate) faces that are on the edge
             //  of a baffle region. These are candidates for merging.
-            List<labelPair> filterDuplicateFaces(const List<labelPair>&) const;
+            List<labelPair> freeStandingBaffles(const List<labelPair>&) const;
 
 
         // Zone handling
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 977f75dbd4..291ae3cb08 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -43,6 +43,7 @@ License
 #include "OFstream.H"
 #include "regionSplit.H"
 #include "removeCells.H"
+#include "unitConversion.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -216,6 +217,45 @@ Foam::label Foam::meshRefinement::getBafflePatch
 }
 
 
+// Check if we are a boundary face and normal of surface does
+// not align with test vector. In this case there'd probably be
+// a freestanding 'baffle' so we might as well not create it.
+// Note that since it is not a proper baffle we cannot detect it
+// afterwards so this code cannot be merged with the
+// filterDuplicateFaces code.
+bool Foam::meshRefinement::validBaffleTopology
+(
+    const label faceI,
+    const vector& n1,
+    const vector& n2,
+    const vector& testDir
+) const
+{
+
+    label patchI = mesh_.boundaryMesh().whichPatch(faceI);
+    if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
+    {
+        return true;
+    }
+    else if (mag(n1&n2) > cos(degToRad(30)))
+    {
+        // Both normals aligned. Check that test vector perpendicularish to
+        // surface normal
+        scalar magTestDir = mag(testDir);
+        if (magTestDir > VSMALL)
+        {
+            if (mag(n1&(testDir/magTestDir)) < cos(degToRad(45)))
+            {
+                //Pout<< "** disabling baffling face "
+                //    << mesh_.faceCentres()[faceI] << endl;
+                return false;
+            }
+        }
+    }
+    return true;
+}
+
+
 // Determine patches for baffles on all intersected unnamed faces
 void Foam::meshRefinement::getBafflePatches
 (
@@ -745,7 +785,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createZoneBaffles
 // Done by counting the number of baffles faces per mesh edge. If edge
 // has 2 boundary faces and both are baffle faces it is the edge of a baffle
 // region.
-Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
+Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
 (
     const List<labelPair>& couples
 ) const
@@ -852,12 +892,112 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
     }
     filteredCouples.setSize(filterI);
 
-    //Info<< "filterDuplicateFaces : from "
+    //Info<< "freeStandingBaffles : from "
     //    << returnReduce(couples.size(), sumOp<label>())
     //    << " down to "
     //    << returnReduce(filteredCouples.size(), sumOp<label>())
     //    << " baffles." << nl << endl;
 
+
+
+//XXXXXX
+//    {
+//        // Collect segments
+//        // ~~~~~~~~~~~~~~~~
+//
+//        pointField start(filteredCouples.size());
+//        pointField end(filteredCouples.size());
+//
+//        const pointField& cellCentres = mesh_.cellCentres();
+//
+//        forAll(filteredCouples, i)
+//        {
+//            const labelPair& couple = couples[i];
+//            start[i] = cellCentres[mesh_.faceOwner()[couple.first()]];
+//            end[i] = cellCentres[mesh_.faceOwner()[couple.second()]];
+//        }
+//
+//        // Extend segments a bit
+//        {
+//            const vectorField smallVec(Foam::sqrt(SMALL)*(end-start));
+//            start -= smallVec;
+//            end += smallVec;
+//        }
+//
+//
+//        // Do test for intersections
+//        // ~~~~~~~~~~~~~~~~~~~~~~~~~
+//        labelList surface1;
+//        List<pointIndexHit> hit1;
+//        labelList region1;
+//        vectorField normal1;
+//
+//        labelList surface2;
+//        List<pointIndexHit> hit2;
+//        labelList region2;
+//        vectorField normal2;
+//
+//        surfaces_.findNearestIntersection
+//        (
+//            surfacesToBaffle,
+//            start,
+//            end,
+//
+//            surface1,
+//            hit1,
+//            region1,
+//            normal1,
+//
+//            surface2,
+//            hit2,
+//            region2,
+//            normal2
+//        );
+//
+//        forAll(testFaces, i)
+//        {
+//            if (hit1[i].hit() && hit2[i].hit())
+//            {
+//                bool createBaffle = true;
+//
+//                label faceI = couples[i].first();
+//                label patchI = mesh_.boundaryMesh().whichPatch(faceI);
+//                if (patchI != -1 && !mesh_.boundaryMesh()[patchI].coupled())
+//                {
+//                    // Check if we are a boundary face and normal of surface
+//                    // does
+//                    // not align with test vector. In this case there'd
+//                    // probably be
+//                    // a freestanding 'baffle' so we might as well not
+//                    // create it.
+//                    // Note that since it is not a proper baffle we cannot
+//                    // detect it
+//                    // afterwards so this code cannot be merged with the
+//                    // filterDuplicateFaces code.
+//                    if (mag(normal1[i]&normal2[i]) > cos(degToRad(30)))
+//                    {
+//                        // Both normals aligned
+//                        vector n = end[i]-start[i];
+//                        scalar magN = mag(n);
+//                        if (magN > VSMALL)
+//                        {
+//                            n /= magN;
+//
+//                            if (mag(normal1[i]&n) < cos(degToRad(45)))
+//                            {
+//                                Pout<< "** disabling baffling face "
+//                                    << mesh_.faceCentres()[faceI] << endl;
+//                                createBaffle = false;
+//                            }
+//                        }
+//                    }
+//                }
+//
+//
+//        }
+//XXXXXX
+
+
     return filteredCouples;
 }
 
@@ -1717,11 +1857,6 @@ void Foam::meshRefinement::baffleAndSplitMesh
         neiPatch
     );
 
-    if (debug)
-    {
-        runTime++;
-    }
-
     createBaffles(ownPatch, neiPatch);
 
     if (debug)
@@ -1874,15 +2009,11 @@ void Foam::meshRefinement::baffleAndSplitMesh
             << "---------------------------" << nl
             << endl;
 
-        if (debug)
-        {
-            runTime++;
-        }
 
         // List of pairs of freestanding baffle faces.
         List<labelPair> couples
         (
-            filterDuplicateFaces    // filter out freestanding baffles
+            freeStandingBaffles    // filter out freestanding baffles
             (
                 getDuplicateFaces   // get all baffles
                 (
@@ -2484,11 +2615,13 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
         // the information already in surfaceIndex_.
 
         labelList surface1;
+        List<pointIndexHit> hit1;
+        vectorField normal1;
         labelList surface2;
+        List<pointIndexHit> hit2;
+        vectorField normal2;
         {
-            List<pointIndexHit> hit1;
             labelList region1;
-            List<pointIndexHit> hit2;
             labelList region2;
             surfaces_.findNearestIntersection
             (
@@ -2511,9 +2644,36 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
 
             if (surface1[i] != -1)
             {
-                // If both hit should probably choose nearest. For later.
-                namedSurfaceIndex[faceI] = surface1[i];
-                nSurfFaces[surface1[i]]++;
+                //- Not allowed not to create baffle - is vital for regioning.
+                //  Have logic instead at erosion!
+                //bool createBaffle = validBaffleTopology
+                //(
+                //    faceI,
+                //    normal1[i],
+                //    normal2[i],
+                //    end[i]-start[i]
+                //);
+                //
+
+
+                // If both hit should probably choose 'nearest'
+                if
+                (
+                    surface2[i] != -1
+                 && (
+                        magSqr(hit2[i].hitPoint())
+                      < magSqr(hit1[i].hitPoint())
+                    )
+                )
+                {
+                    namedSurfaceIndex[faceI] = surface2[i];
+                    nSurfFaces[surface2[i]]++;
+                }
+                else
+                {
+                    namedSurfaceIndex[faceI] = surface1[i];
+                    nSurfFaces[surface1[i]]++;
+                }
             }
             else if (surface2[i] != -1)
             {