diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index 07ff697df6..adb7c38590 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -910,6 +910,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
     Info<< "Detecting near surfaces ..." << endl;
 
     const pointField& localPoints = pp.localPoints();
+    const labelList& meshPoints = pp.meshPoints();
     const refinementSurfaces& surfaces = meshRefiner_.surfaces();
     const fvMesh& mesh = meshRefiner_.mesh();
 
@@ -1220,6 +1221,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
     }
 
 
+    const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
     label nOverride = 0;
 
     // 1. All points to non-interface surfaces
@@ -1332,7 +1334,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
                 }
             }
 
-            if (override)
+            if (override && isMasterPoint[meshPoints[pointI]])
             {
                 nOverride++;
             }
@@ -1399,8 +1401,6 @@ void Foam::autoSnapDriver::detectNearSurfaces
             );
 
 
-            label nOverride = 0;
-
             forAll(hit1, i)
             {
                 label pointI = zonePointIndices[i];
@@ -1472,7 +1472,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
                     }
                 }
 
-                if (override)
+                if (override && isMasterPoint[meshPoints[pointI]])
                 {
                     nOverride++;
                 }
@@ -1690,7 +1690,13 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
             scalarField magDisp(mag(patchDisp));
 
             Info<< "Wanted displacement : average:"
-                << gSum(magDisp)/returnReduce(patchDisp.size(), sumOp<label>())
+                <<  meshRefinement::gAverage
+                    (
+                        mesh,
+                        syncTools::getMasterPoints(mesh),
+                        pp.meshPoints(),
+                        magDisp
+                    )
                 << " min:" << gMin(magDisp)
                 << " max:" << gMax(magDisp) << endl;
         }
@@ -2548,25 +2554,30 @@ void Foam::autoSnapDriver::doSnap
                 adaptPatchIDs
             )
         );
-        indirectPrimitivePatch& pp = ppPtr();
+
 
         // Distance to attract to nearest feature on surface
-        const scalarField snapDist(calcSnapDistance(mesh, snapParams, pp));
+        const scalarField snapDist(calcSnapDistance(mesh, snapParams, ppPtr()));
 
 
         // Construct iterative mesh mover.
         Info<< "Constructing mesh displacer ..." << endl;
         Info<< "Using mesh parameters " << motionDict << nl << endl;
 
-        const pointMesh& pMesh = pointMesh::New(mesh);
-
-        motionSmoother meshMover
+        autoPtr<motionSmoother> meshMoverPtr
         (
-            mesh,
-            pp,
-            adaptPatchIDs,
-            meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs),
-            motionDict
+            new motionSmoother
+            (
+                mesh,
+                ppPtr(),
+                adaptPatchIDs,
+                meshRefinement::makeDisplacementField
+                (
+                    pointMesh::New(mesh),
+                    adaptPatchIDs
+                ),
+                motionDict
+            )
         );
 
 
@@ -2595,16 +2606,95 @@ void Foam::autoSnapDriver::doSnap
             snapParams,
             nInitErrors,
             baffles,
-            meshMover
+            meshMoverPtr()
         );
 
 
+
+        //- Only if in feature attraction mode:
+        // Nearest feature
+        vectorField patchAttraction;
+        // Constraints at feature
+        List<pointConstraint> patchConstraints;
+
+
         for (label iter = 0; iter < nFeatIter; iter++)
         {
             Info<< nl
                 << "Morph iteration " << iter << nl
                 << "-----------------" << endl;
 
+
+            //if (iter > 0 && iter == nFeatIter/2)
+            //{
+            //    Info<< "Splitting diagonal attractions" << endl;
+            //    const labelList& bFaces = ppPtr().addressing();
+            //
+            //    DynamicList<label> splitFaces(bFaces.size());
+            //    DynamicList<labelPair> splits(bFaces.size());
+            //
+            //    forAll(bFaces, faceI)
+            //    {
+            //        const labelPair split
+            //        (
+            //            findDiagonalAttraction
+            //            (
+            //                ppPtr(),
+            //                patchAttraction,
+            //                patchConstraints,
+            //                faceI
+            //            )
+            //        );
+            //
+            //        if (split != labelPair(-1, -1))
+            //        {
+            //            splitFaces.append(bFaces[faceI]);
+            //            splits.append(split);
+            //        }
+            //    }
+            //
+            //    autoPtr<mapPolyMesh> mapPtr = meshRefiner_.splitFaces
+            //    (
+            //        splitFaces,
+            //        splits
+            //    );
+            //
+            //    const labelList& faceMap = mapPtr().faceMap();
+            //    meshRefinement::updateList(faceMap, -1, duplicateFace);
+            //    const labelList& reverseFaceMap = mapPtr().reverseFaceMap();
+            //    forAll(baffles, i)
+            //    {
+            //        labelPair& baffle = baffles[i];
+            //        baffle.first() = reverseFaceMap[baffle.first()];
+            //        baffle.second() = reverseFaceMap[baffle.second()];
+            //    }
+            //
+            //    meshMoverPtr.clear();
+            //    ppPtr.clear();
+            //
+            //    ppPtr = meshRefinement::makePatch(mesh, adaptPatchIDs);
+            //    meshMoverPtr.reset
+            //    (
+            //        new motionSmoother
+            //        (
+            //            mesh,
+            //            ppPtr(),
+            //            adaptPatchIDs,
+            //            meshRefinement::makeDisplacementField
+            //            (
+            //                pointMesh::New(mesh),
+            //                adaptPatchIDs
+            //            ),
+            //            motionDict
+            //        )
+            //    );
+            //}
+
+
+            indirectPrimitivePatch& pp = ppPtr();
+            motionSmoother& meshMover = meshMoverPtr();
+
+
             // Calculate displacement at every patch point. Insert into
             // meshMover.
             // Calculate displacement at every patch point
@@ -2652,7 +2742,9 @@ void Foam::autoSnapDriver::doSnap
                     scalar(iter+1)/nFeatIter,
                     snapDist,
                     disp,
-                    meshMover
+                    meshMover,
+                    patchAttraction,
+                    patchConstraints
                 );
             }
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index 77327c0689..b8282bedfc 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -125,7 +125,9 @@ class autoSnapDriver
 
                 void smoothAndConstrain
                 (
+                    const PackedBoolList& isMasterEdge,
                     const indirectPrimitivePatch& pp,
+                    const labelList& meshEdges,
                     const List<pointConstraint>& constraints,
                     vectorField& disp
                 ) const;
@@ -196,6 +198,16 @@ class autoSnapDriver
                     List<pointConstraint>& patchConstraints
                 ) const;
 
+                //- Detect any diagonal attraction. Returns indices in face
+                //  or (-1, -1) if none
+                labelPair findDiagonalAttraction
+                (
+                    const indirectPrimitivePatch& pp,
+                    const vectorField& patchAttraction,
+                    const List<pointConstraint>& patchConstraints,
+                    const label faceI
+                ) const;
+
                 //- Return hit if on multiple points
                 pointIndexHit findMultiPatchPoint
                 (
@@ -371,7 +383,9 @@ class autoSnapDriver
                     const scalar featureAttract,
                     const scalarField& snapDist,
                     const vectorField& nearestDisp,
-                    motionSmoother& meshMover
+                    motionSmoother& meshMover,
+                    vectorField& patchAttraction,
+                    List<pointConstraint>& patchConstraints
                 ) const;
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index 887674ea7a..32133ad5f8 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -38,50 +38,18 @@ License
 #include "treeDataPoint.H"
 #include "indexedOctree.H"
 #include "snapParameters.H"
+#include "PatchTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    class listTransform
-    {
-    public:
-
-        void operator()
-        (
-            const vectorTensorTransform& vt,
-            const bool forward,
-            List<List<point> >& fld
-        ) const
-        {
-            const tensor T
-            (
-                forward
-              ? vt.R()
-              : vt.R().T()
-            );
-
-            forAll(fld, i)
-            {
-                List<point>& elems = fld[i];
-                forAll(elems, elemI)
-                {
-                    elems[elemI] = transform(T, elems[elemI]);
-                }
-            }
-        }
-    };
-
     template<class T>
     class listPlusEqOp
     {
     public:
 
-        void operator()
-        (
-            List<T>& x,
-            const List<T>& y
-        ) const
+        void operator()(List<T>& x, const List<T>& y) const
         {
             label sz = x.size();
             x.setSize(sz+y.size());
@@ -159,7 +127,9 @@ bool Foam::autoSnapDriver::isFeaturePoint
 
 void Foam::autoSnapDriver::smoothAndConstrain
 (
+    const PackedBoolList& isMasterEdge,
     const indirectPrimitivePatch& pp,
+    const labelList& meshEdges,
     const List<pointConstraint>& constraints,
     vectorField& disp
 ) const
@@ -197,11 +167,16 @@ void Foam::autoSnapDriver::smoothAndConstrain
             {
                 forAll(pEdges, i)
                 {
-                    label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
-                    if (constraints[nbrPointI].first() >= nConstraints)
+                    label edgeI = pEdges[i];
+
+                    if (isMasterEdge[meshEdges[edgeI]])
                     {
-                        dispSum[pointI] += disp[nbrPointI];
-                        dispCount[pointI]++;
+                        label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
+                        if (constraints[nbrPointI].first() >= nConstraints)
+                        {
+                            dispSum[pointI] += disp[nbrPointI];
+                            dispCount[pointI]++;
+                        }
                     }
                 }
             }
@@ -564,6 +539,9 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
 {
     const fvMesh& mesh = meshRefiner_.mesh();
 
+    const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
+
+
     // For now just get all surrounding face data. Expensive - should just
     // store and sync data on coupled points only
     // (see e.g PatchToolsNormals.C)
@@ -583,7 +561,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
         forAll(pFaces, i)
         {
             label faceI = pFaces[i];
-            if (faceSurfaceGlobalRegion[faceI] != -1)
+            if (isMasterFace[faceI] && faceSurfaceGlobalRegion[faceI] != -1)
             {
                 nFaces++;
             }
@@ -605,7 +583,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
             label faceI = pFaces[i];
             label globalRegionI = faceSurfaceGlobalRegion[faceI];
 
-            if (globalRegionI != -1)
+            if (isMasterFace[faceI] && globalRegionI != -1)
             {
                 pNormals[nFaces] = faceSurfaceNormal[faceI];
                 pDisp[nFaces] = faceDisp[faceI];
@@ -694,7 +672,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
         pointFaceSurfNormals,
         listPlusEqOp<point>(),
         List<point>(),
-        listTransform()
+        mapDistribute::transform()
     );
     syncTools::syncPointList
     (
@@ -703,7 +681,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
         pointFaceDisp,
         listPlusEqOp<point>(),
         List<point>(),
-        listTransform()
+        mapDistribute::transform()
     );
     syncTools::syncPointList
     (
@@ -712,7 +690,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
         pointFaceCentres,
         listPlusEqOp<point>(),
         List<point>(),
-        listTransform()
+        mapDistribute::transformPosition()
     );
     syncTools::syncPointList
     (
@@ -722,6 +700,25 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
         listPlusEqOp<label>(),
         List<label>()
     );
+
+
+    // Sort the data according to the face centres. This is only so we get
+    // consistent behaviour serial and parallel.
+    labelList visitOrder;
+    forAll(pointFaceDisp, pointI)
+    {
+        List<point>& pNormals = pointFaceSurfNormals[pointI];
+        List<point>& pDisp = pointFaceDisp[pointI];
+        List<point>& pFc = pointFaceCentres[pointI];
+        labelList& pFid = pointFacePatchID[pointI];
+
+        sortedOrder(mag(pFc)(), visitOrder);
+
+        pNormals = List<point>(pNormals, visitOrder);
+        pDisp = List<point>(pDisp, visitOrder);
+        pFc = List<point>(pFc, visitOrder);
+        pFid = UIndirectList<label>(pFid, visitOrder);
+    }
 }
 
 
@@ -1360,6 +1357,70 @@ void Foam::autoSnapDriver::stringFeatureEdges
 }
 
 
+// If only two attractions and across face return the face indices
+Foam::labelPair Foam::autoSnapDriver::findDiagonalAttraction
+(
+    const indirectPrimitivePatch& pp,
+    const vectorField& patchAttraction,
+    const List<pointConstraint>& patchConstraints,
+    const label faceI
+) const
+{
+    const face& f = pp.localFaces()[faceI];
+    // For now just detect any attraction. Improve this to look at
+    // actual attraction position and orientation
+
+    labelPair attractIndices(-1, -1);
+
+    if (f.size() >= 4)
+    {
+        forAll(f, fp)
+        {
+            label pointI = f[fp];
+            if (patchConstraints[pointI].first() >= 2)
+            {
+                // Attract to feature edge or point
+                if (attractIndices[0] == -1)
+                {
+                    // First attraction. Store
+                    attractIndices[0] = fp;
+                }
+                else if (attractIndices[1] == -1)
+                {
+                    // Second attraction. Check if not consecutive to first
+                    // attraction
+                    label fp0 = attractIndices[0];
+                    if (f.fcIndex(fp0) == fp || f.fcIndex(fp) == fp0)
+                    {
+                        // Consecutive. Skip.
+                        attractIndices = labelPair(-1, -1);
+                        break;
+                    }
+                    else
+                    {
+                        attractIndices[1] = fp;
+                    }
+                }
+                else
+                {
+                    // More than two attractions. Skip.
+                    attractIndices = labelPair(-1, -1);
+                    break;
+                }
+            }
+        }
+
+
+        if (attractIndices[1] == -1)
+        {
+            // Found only one attraction. Skip.
+            attractIndices = labelPair(-1, -1);
+        }
+    }
+    return attractIndices;
+}
+
+
 Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
 (
     const indirectPrimitivePatch& pp,
@@ -1929,7 +1990,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                 edgeFaceNormals,
                 listPlusEqOp<point>(),
                 List<point>(),
-                listTransform()
+                mapDistribute::transform()
             );
         }
 
@@ -2778,7 +2839,9 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     const scalar featureAttract,
     const scalarField& snapDist,
     const vectorField& nearestDisp,
-    motionSmoother& meshMover
+    motionSmoother& meshMover,
+    vectorField& patchAttraction,
+    List<pointConstraint>& patchConstraints
 ) const
 {
     const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
@@ -2851,7 +2914,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // - faceSurfaceNormal
     // - faceDisp
-    // - faceCentres&faceNormal
+    // - faceCentres
     List<List<point> > pointFaceSurfNormals(pp.nPoints());
     List<List<point> > pointFaceDisp(pp.nPoints());
     List<List<point> > pointFaceCentres(pp.nPoints());
@@ -2884,10 +2947,11 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     // here.
 
     // Nearest feature
-    vectorField patchAttraction(localPoints.size(), vector::zero);
+    patchAttraction.setSize(localPoints.size());
+    patchAttraction = vector::zero;
     // Constraints at feature
-    List<pointConstraint> patchConstraints(localPoints.size());
-
+    patchConstraints.setSize(localPoints.size());
+    patchConstraints = pointConstraint();
 
     if (implicitFeatureAttraction)
     {
@@ -2951,15 +3015,29 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         patchConstraints
     );
 
+    const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
+    {
+        vector avgPatchDisp = meshRefinement::gAverage
+        (
+            mesh,
+            isMasterPoint,
+            pp.meshPoints(),
+            patchDisp
+        );
+        vector avgPatchAttr = meshRefinement::gAverage
+        (
+            mesh,
+            isMasterPoint,
+            pp.meshPoints(),
+            patchAttraction
+        );
 
-    Info<< "Attraction:" << endl
-        << "     linear   : max:" << gMax(patchDisp)
-        << " avg:" << gAverage(patchDisp)
-        << endl
-        << "     feature  : max:" << gMax(patchAttraction)
-        << " avg:" << gAverage(patchAttraction)
-        << endl;
-
+        Info<< "Attraction:" << endl
+            << "     linear   : max:" << gMaxMagSqr(patchDisp)
+            << " avg:" << avgPatchDisp << endl
+            << "     feature  : max:" << gMaxMagSqr(patchAttraction)
+            << " avg:" << avgPatchAttr << endl;
+    }
 
     // So now we have:
     // - patchDisp          : point movement to go to nearest point on surface
@@ -2986,37 +3064,45 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 
     // Count
     {
+        const labelList& meshPoints = pp.meshPoints();
+
+        label nMasterPoints = 0;
         label nPlanar = 0;
         label nEdge = 0;
         label nPoint = 0;
 
         forAll(patchConstraints, pointI)
         {
-            if (patchConstraints[pointI].first() == 1)
+            if (isMasterPoint[meshPoints[pointI]])
             {
-                nPlanar++;
-            }
-            else if (patchConstraints[pointI].first() == 2)
-            {
-                nEdge++;
-            }
-            else if (patchConstraints[pointI].first() == 3)
-            {
-                nPoint++;
+                nMasterPoints++;
+
+                if (patchConstraints[pointI].first() == 1)
+                {
+                    nPlanar++;
+                }
+                else if (patchConstraints[pointI].first() == 2)
+                {
+                    nEdge++;
+                }
+                else if (patchConstraints[pointI].first() == 3)
+                {
+                    nPoint++;
+                }
             }
         }
 
-        label nTotPoints = returnReduce(pp.nPoints(), sumOp<label>());
+        reduce(nMasterPoints, sumOp<label>());
         reduce(nPlanar, sumOp<label>());
         reduce(nEdge, sumOp<label>());
         reduce(nPoint, sumOp<label>());
-        Info<< "Feature analysis : total points:"
-            << nTotPoints
+        Info<< "Feature analysis : total master points:"
+            << nMasterPoints
             << " attraction to :" << nl
             << "    feature point   : " << nPoint << nl
             << "    feature edge    : " << nEdge << nl
             << "    nearest surface : " << nPlanar << nl
-            << "    rest            : " << nTotPoints-nPoint-nEdge-nPlanar
+            << "    rest            : " << nMasterPoints-nPoint-nEdge-nPlanar
             << nl
             << endl;
     }
@@ -3035,11 +3121,29 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 
     if (featureAttract < 1-0.001)
     {
+        const PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
+
+        const vectorField pointNormals
+        (
+            PatchTools::pointNormals
+            (
+                mesh,
+                pp
+            )
+        );
+        const labelList meshEdges
+        (
+            pp.meshEdges(mesh.edges(), mesh.pointEdges())
+        );
+
+
         // 1. Smoothed all displacement
         vectorField smoothedPatchDisp = patchDisp;
         smoothAndConstrain
         (
+            isMasterEdge,
             pp,
+            meshEdges,
             patchConstraints,
             smoothedPatchDisp
         );
@@ -3047,16 +3151,18 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 
         // 2. Smoothed tangential component
         vectorField tangPatchDisp = patchDisp;
-        tangPatchDisp -= (pp.pointNormals() & patchDisp) * pp.pointNormals();
+        tangPatchDisp -= (pointNormals & patchDisp) * pointNormals;
         smoothAndConstrain
         (
+            isMasterEdge,
             pp,
+            meshEdges,
             patchConstraints,
             tangPatchDisp
         );
 
         // Re-add normal component
-        tangPatchDisp += (pp.pointNormals() & patchDisp) * pp.pointNormals();
+        tangPatchDisp += (pointNormals & patchDisp) * pointNormals;
 
         if (debug&meshRefinement::OBJINTERSECTIONS)
         {
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 8a5f327b41..873167bc4f 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -259,23 +259,14 @@ private:
                 label& nRefine
             ) const;
 
-            //- Mark cell if local curvature > curvature or
-            //  markDifferingRegions = true and intersections with different
-            //  regions.
-            bool checkCurvature
+            //- Helper: count number of normals1 that are in normals2
+            label countMatches
             (
-                const scalar curvature,
-                const label nAllowRefine,
-                const label surfaceLevel,
-                const vector& surfaceNormal,
-                const label cellI,
-                label& cellMaxLevel,
-                vector& cellMaxNormal,
-                labelList& refineCell,
-                label& nRefine
+                const List<point>& normals1,
+                const List<point>& normals2,
+                const scalar tol = 1e-6
             ) const;
 
-
             //- Mark cells for surface curvature based refinement. Marks if
             //  local curvature > curvature or if on different regions
             //  (markDifferingRegions)
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index 7444493e29..ce84ceec38 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -38,7 +38,74 @@ License
 #include "featureEdgeMesh.H"
 #include "Cloud.H"
 //#include "globalIndex.H"
-//#include "OBJstream.H"
+#include "OBJstream.H"
+#include "cellSet.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    //- To compare normals
+    static bool less(const vector& x, const vector& y)
+    {
+        for (direction i = 0; i < vector::nComponents; i++)
+        {
+            if (x[i] < y[i])
+            {
+                return true;
+            }
+            else if (x[i] > y[i])
+            {
+                return false;
+            }
+        }
+        // All components the same
+        return false;
+    }
+
+
+    //- To compare normals
+    class normalLess
+    {
+        const vectorList& values_;
+
+    public:
+
+        normalLess(const vectorList& values)
+        :
+            values_(values)
+        {}
+
+        bool operator()(const label a, const label b) const
+        {
+            return less(values_[a], values_[b]);
+        }
+    };
+
+
+    //- template specialization for pTraits<labelList> so we can have fields
+    template<>
+    class pTraits<labelList>
+    {
+
+    public:
+
+        //- Component type
+        typedef labelList cmptType;
+    };
+
+    //- template specialization for pTraits<labelList> so we can have fields
+    template<>
+    class pTraits<vectorList>
+    {
+
+    public:
+
+        //- Component type
+        typedef vectorList cmptType;
+    };
+}
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -954,64 +1021,33 @@ Foam::label Foam::meshRefinement::markSurfaceRefinement
 }
 
 
-// Checks if multiple intersections of a cell (by a surface with a higher
-// max than the cell level) and if so if the normals at these intersections
-// make a large angle.
-// Returns false if the nRefine limit has been reached, true otherwise.
-bool Foam::meshRefinement::checkCurvature
+// Count number of matches of first argument in second argument
+Foam::label Foam::meshRefinement::countMatches
 (
-    const scalar curvature,
-    const label nAllowRefine,
-
-    const label surfaceLevel,   // current intersection max level
-    const vector& surfaceNormal,// current intersection normal
-
-    const label cellI,
-
-    label& cellMaxLevel,        // cached max surface level for this cell
-    vector& cellMaxNormal,      // cached surface normal for this cell
-
-    labelList& refineCell,
-    label& nRefine
+    const List<point>& normals1,
+    const List<point>& normals2,
+    const scalar tol
 ) const
 {
-    const labelList& cellLevel = meshCutter_.cellLevel();
+    label nMatches = 0;
 
-    // Test if surface applicable
-    if (surfaceLevel > cellLevel[cellI])
+    forAll(normals1, i)
     {
-        if (cellMaxLevel == -1)
-        {
-            // First visit of cell. Store
-            cellMaxLevel = surfaceLevel;
-            cellMaxNormal = surfaceNormal;
-        }
-        else
-        {
-            // Second or more visit. Check curvature.
-            if ((cellMaxNormal & surfaceNormal) < curvature)
-            {
-                return markForRefine
-                (
-                    surfaceLevel,   // mark with any non-neg number.
-                    nAllowRefine,
-                    refineCell[cellI],
-                    nRefine
-                );
-            }
+        const vector& n1 = normals1[i];
 
-            // Set normal to that of highest surface. Not really necessary
-            // over here but we reuse cellMax info when doing coupled faces.
-            if (surfaceLevel > cellMaxLevel)
+        forAll(normals2, j)
+        {
+            const vector& n2 = normals2[j];
+
+            if (magSqr(n1-n2) < tol)
             {
-                cellMaxLevel = surfaceLevel;
-                cellMaxNormal = surfaceNormal;
+                nMatches++;
+                break;
             }
         }
     }
 
-    // Did not reach refinement limit.
-    return true;
+    return nMatches;
 }
 
 
@@ -1039,6 +1075,9 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
     // on a different surface gets refined (if its current level etc.)
 
 
+    const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
+
+
     // Collect candidate faces (i.e. intersecting any surface and
     // owner/neighbour not yet refined.
     labelList testFaces(getRefineCandidateFaces(refineCell));
@@ -1068,6 +1107,12 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
 
             start[i] = cellCentres[own];
             end[i] = neiCc[bFaceI];
+
+            if (!isMasterFace[faceI])
+            {
+                Swap(start[i], end[i]);
+            }
+
             minLevel[i] = min(cellLevel[own], neiLevel[bFaceI]);
         }
     }
@@ -1081,10 +1126,9 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
 
 
     // Test for all intersections (with surfaces of higher max level than
-    // minLevel) and cache per cell the max surface level and the local normal
-    // on that surface.
-    labelList cellMaxLevel(mesh_.nCells(), -1);
-    vectorField cellMaxNormal(mesh_.nCells(), vector::zero);
+    // minLevel) and cache per cell the interesting inter
+    labelListList cellSurfLevels(mesh_.nCells());
+    List<vectorList> cellSurfNormals(mesh_.nCells());
 
     {
         // Per segment the normals of the surfaces
@@ -1104,12 +1148,29 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
             surfaceNormal,
             surfaceLevel
         );
+
+
+        // Sort the data according to surface location. This will guarantee
+        // that on coupled faces both sides visit the intersections in
+        // the same order so will decide the same
+        labelList visitOrder;
+        forAll(surfaceNormal, pointI)
+        {
+            vectorList& pNormals = surfaceNormal[pointI];
+            labelList& pLevel = surfaceLevel[pointI];
+
+            sortedOrder(pNormals, visitOrder, normalLess(pNormals));
+
+            pNormals = List<point>(pNormals, visitOrder);
+            pLevel = UIndirectList<label>(pLevel, visitOrder);
+        }
+
         // Clear out unnecessary data
         start.clear();
         end.clear();
         minLevel.clear();
 
-        // Extract per cell information on the surface with the highest max
+        // Convert face-wise data to cell.
         forAll(testFaces, i)
         {
             label faceI = testFaces[i];
@@ -1118,164 +1179,280 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
             const vectorList& fNormals = surfaceNormal[i];
             const labelList& fLevels = surfaceLevel[i];
 
-            forAll(fLevels, hitI)
+            forAll(fNormals, hitI)
             {
-                checkCurvature
-                (
-                    curvature,
-                    nAllowRefine,
-
-                    fLevels[hitI],
-                    fNormals[hitI],
-
-                    own,
-                    cellMaxLevel[own],
-                    cellMaxNormal[own],
-
-                    refineCell,
-                    nRefine
-                );
-            }
-
-            if (mesh_.isInternalFace(faceI))
-            {
-                label nei = mesh_.faceNeighbour()[faceI];
-
-                forAll(fLevels, hitI)
+                if (fLevels[hitI] > cellLevel[own])
                 {
-                    checkCurvature
-                    (
-                        curvature,
-                        nAllowRefine,
+                    cellSurfLevels[own].append(fLevels[hitI]);
+                    cellSurfNormals[own].append(fNormals[hitI]);
+                }
 
-                        fLevels[hitI],
-                        fNormals[hitI],
-
-                        nei,
-                        cellMaxLevel[nei],
-                        cellMaxNormal[nei],
-
-                        refineCell,
-                        nRefine
-                    );
+                if (mesh_.isInternalFace(faceI))
+                {
+                    label nei = mesh_.faceNeighbour()[faceI];
+                    if (fLevels[hitI] > cellLevel[nei])
+                    {
+                        cellSurfLevels[nei].append(fLevels[hitI]);
+                        cellSurfNormals[nei].append(fNormals[hitI]);
+                    }
                 }
             }
         }
     }
 
-    // 2. Find out a measure of surface curvature and region edges.
-    // Send over surface region and surface normal to neighbour cell.
 
-    labelList neiBndMaxLevel(mesh_.nFaces()-mesh_.nInternalFaces());
-    vectorField neiBndMaxNormal(mesh_.nFaces()-mesh_.nInternalFaces());
 
-    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
+    // Bit of statistics
     {
-        label bFaceI = faceI-mesh_.nInternalFaces();
-        label own = mesh_.faceOwner()[faceI];
-
-        neiBndMaxLevel[bFaceI] = cellMaxLevel[own];
-        neiBndMaxNormal[bFaceI] = cellMaxNormal[own];
+        label nSet = 0;
+        label nNormals = 9;
+        forAll(cellSurfNormals, cellI)
+        {
+            const vectorList& normals = cellSurfNormals[cellI];
+            if (normals.size())
+            {
+                nSet++;
+                nNormals += normals.size();
+            }
+        }
+        reduce(nSet, sumOp<label>());
+        reduce(nNormals, sumOp<label>());
+        Info<< "markSurfaceCurvatureRefinement :"
+            << " cells:" << mesh_.globalData().nTotalCells()
+            << " of which with normals:" << nSet
+            << " ; total normals stored:" << nNormals
+            << endl;
     }
-    syncTools::swapBoundaryFaceList(mesh_, neiBndMaxLevel);
-    syncTools::swapBoundaryFaceList(mesh_, neiBndMaxNormal);
 
-    // Loop over all faces. Could only be checkFaces.. except if they're coupled
+
+
+    bool reachedLimit = false;
+
+
+    // 1. Check normals per cell
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    for
+    (
+        label cellI = 0;
+        !reachedLimit && cellI < cellSurfNormals.size();
+        cellI++
+    )
+    {
+        const vectorList& normals = cellSurfNormals[cellI];
+        const labelList& levels = cellSurfLevels[cellI];
+
+        // n^2 comparison of all normals in a cell
+        for (label i = 0; !reachedLimit && i < normals.size(); i++)
+        {
+            for (label j = i+1; !reachedLimit && j < normals.size(); j++)
+            {
+                if ((normals[i] & normals[j]) < curvature)
+                {
+                    label maxLevel = max(levels[i], levels[j]);
+
+                    if (cellLevel[cellI] < maxLevel)
+                    {
+                        if
+                        (
+                            !markForRefine
+                            (
+                                maxLevel,
+                                nAllowRefine,
+                                refineCell[cellI],
+                                nRefine
+                            )
+                        )
+                        {
+                            if (debug)
+                            {
+                                Pout<< "Stopped refining since reaching my cell"
+                                    << " limit of " << mesh_.nCells()+7*nRefine
+                                    << endl;
+                            }
+                            reachedLimit = true;
+                            break;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+
+
+    // 2. Find out a measure of surface curvature
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Look at normals between neighbouring surfaces
+    // Loop over all faces. Could only be checkFaces, except if they're coupled
 
     // Internal faces
-    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    for
+    (
+        label faceI = 0;
+        !reachedLimit && faceI < mesh_.nInternalFaces();
+        faceI++
+    )
     {
         label own = mesh_.faceOwner()[faceI];
         label nei = mesh_.faceNeighbour()[faceI];
 
-        if (cellMaxLevel[own] != -1 && cellMaxLevel[nei] != -1)
-        {
-            // Have valid data on both sides. Check curvature.
-            if ((cellMaxNormal[own] & cellMaxNormal[nei]) < curvature)
-            {
-                // See which side to refine
-                if (cellLevel[own] < cellMaxLevel[own])
-                {
-                    if
-                    (
-                        !markForRefine
-                        (
-                            cellMaxLevel[own],
-                            nAllowRefine,
-                            refineCell[own],
-                            nRefine
-                        )
-                    )
-                    {
-                        if (debug)
-                        {
-                            Pout<< "Stopped refining since reaching my cell"
-                                << " limit of " << mesh_.nCells()+7*nRefine
-                                << endl;
-                        }
-                        break;
-                    }
-                }
+        const vectorList& ownNormals = cellSurfNormals[own];
+        const labelList& ownLevels = cellSurfLevels[own];
+        const vectorList& neiNormals = cellSurfNormals[nei];
+        const labelList& neiLevels = cellSurfLevels[nei];
 
-                if (cellLevel[nei] < cellMaxLevel[nei])
+
+        // Special case: owner normals set is a subset of the neighbour
+        // normals. Do not do curvature refinement since other cell's normals
+        // do not add any information. Avoids weird corner extensions of
+        // refinement regions.
+        bool ownIsSubset =
+            countMatches(ownNormals, neiNormals)
+         == ownNormals.size();
+
+        bool neiIsSubset =
+            countMatches(neiNormals, ownNormals)
+         == neiNormals.size();
+
+
+        if (!ownIsSubset && !neiIsSubset)
+        {
+            // n^2 comparison of between ownNormals and neiNormals
+            for (label i = 0; !reachedLimit &&  i < ownNormals.size(); i++)
+            {
+                for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
                 {
-                    if
-                    (
-                        !markForRefine
-                        (
-                            cellMaxLevel[nei],
-                            nAllowRefine,
-                            refineCell[nei],
-                            nRefine
-                        )
-                    )
+                    // Have valid data on both sides. Check curvature.
+                    if ((ownNormals[i] & neiNormals[j]) < curvature)
                     {
-                        if (debug)
+                        // See which side to refine.
+                        if (cellLevel[own] < ownLevels[i])
                         {
-                            Pout<< "Stopped refining since reaching my cell"
-                                << " limit of " << mesh_.nCells()+7*nRefine
-                                << endl;
+                            if
+                            (
+                                !markForRefine
+                                (
+                                    ownLevels[i],
+                                    nAllowRefine,
+                                    refineCell[own],
+                                    nRefine
+                                )
+                            )
+                            {
+                                if (debug)
+                                {
+                                    Pout<< "Stopped refining since reaching"
+                                        << " my cell limit of "
+                                        << mesh_.nCells()+7*nRefine << endl;
+                                }
+                                reachedLimit = true;
+                                break;
+                            }
+                        }
+                        if (cellLevel[nei] < neiLevels[j])
+                        {
+                            if
+                            (
+                                !markForRefine
+                                (
+                                    neiLevels[j],
+                                    nAllowRefine,
+                                    refineCell[nei],
+                                    nRefine
+                                )
+                            )
+                            {
+                                if (debug)
+                                {
+                                    Pout<< "Stopped refining since reaching"
+                                        << " my cell limit of "
+                                        << mesh_.nCells()+7*nRefine << endl;
+                                }
+                                reachedLimit = true;
+                                break;
+                            }
                         }
-                        break;
                     }
                 }
             }
         }
     }
+
+
+    // Send over surface normal to neighbour cell.
+    List<vectorList> neiSurfaceNormals;
+    syncTools::swapBoundaryCellList(mesh_, cellSurfNormals, neiSurfaceNormals);
+
     // Boundary faces
-    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
+    for
+    (
+        label faceI = mesh_.nInternalFaces();
+        !reachedLimit && faceI < mesh_.nFaces();
+        faceI++
+    )
     {
         label own = mesh_.faceOwner()[faceI];
         label bFaceI = faceI - mesh_.nInternalFaces();
 
-        if (cellLevel[own] < cellMaxLevel[own] && neiBndMaxLevel[bFaceI] != -1)
+        const vectorList& ownNormals = cellSurfNormals[own];
+        const labelList& ownLevels = cellSurfLevels[own];
+        const vectorList& neiNormals = neiSurfaceNormals[bFaceI];
+
+        // Special case: owner normals set is a subset of the neighbour
+        // normals. Do not do curvature refinement since other cell's normals
+        // do not add any information. Avoids weird corner extensions of
+        // refinement regions.
+        bool ownIsSubset =
+            countMatches(ownNormals, neiNormals)
+         == ownNormals.size();
+
+        bool neiIsSubset =
+            countMatches(neiNormals, ownNormals)
+         == neiNormals.size();
+
+
+        if (!ownIsSubset && !neiIsSubset)
         {
-            // Have valid data on both sides. Check curvature.
-            if ((cellMaxNormal[own] & neiBndMaxNormal[bFaceI]) < curvature)
+            // n^2 comparison of between ownNormals and neiNormals
+            for (label i = 0; !reachedLimit &&  i < ownNormals.size(); i++)
             {
-                if
-                (
-                    !markForRefine
-                    (
-                        cellMaxLevel[own],
-                        nAllowRefine,
-                        refineCell[own],
-                        nRefine
-                    )
-                )
+                for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
                 {
-                    if (debug)
+                    // Have valid data on both sides. Check curvature.
+                    if ((ownNormals[i] & neiNormals[j]) < curvature)
                     {
-                        Pout<< "Stopped refining since reaching my cell"
-                            << " limit of " << mesh_.nCells()+7*nRefine
-                            << endl;
+                        if (cellLevel[own] < ownLevels[i])
+                        {
+                            if
+                            (
+                                !markForRefine
+                                (
+                                    ownLevels[i],
+                                    nAllowRefine,
+                                    refineCell[own],
+                                    nRefine
+                                )
+                            )
+                            {
+                                if (debug)
+                                {
+                                    Pout<< "Stopped refining since reaching"
+                                        << " my cell limit of "
+                                        << mesh_.nCells()+7*nRefine
+                                        << endl;
+                                }
+                                reachedLimit = true;
+                                break;
+                            }
+                        }
                     }
-                    break;
                 }
             }
         }
     }
 
+
     if
     (
         returnReduce(nRefine, sumOp<label>())