From 1a0a6b51dd33520014195b7c332af5e9b01378d5 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 13 Sep 2013 16:18:04 +0100
Subject: [PATCH] ENH: snappyHexMesh: detect small distance snapping

---
 .../autoHexMeshDriver/autoSnapDriver.C        | 750 +++++++++++++++++-
 .../autoHexMeshDriver/autoSnapDriver.H        |  23 +-
 .../autoHexMeshDriver/autoSnapDriverFeature.C | 168 +++-
 .../snapParameters/snapParameters.C           |   9 +-
 .../snapParameters/snapParameters.H           |   9 +-
 .../meshRefinement/meshRefinement.H           |  30 +-
 .../meshRefinementProblemCells.C              |   6 +-
 .../meshRefinement/meshRefinementRefine.C     |  56 +-
 8 files changed, 987 insertions(+), 64 deletions(-)

diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index e2b2989889..81fb2a36ca 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -792,11 +792,636 @@ Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
 }
 
 
+Foam::tmp<Foam::pointField> Foam::autoSnapDriver::avgCellCentres
+(
+    const fvMesh& mesh,
+    const indirectPrimitivePatch& pp
+)
+{
+    const labelListList& pointFaces = pp.pointFaces();
+
+
+    tmp<pointField> tavgBoundary
+    (
+        new pointField(pointFaces.size(), vector::zero)
+    );
+    pointField& avgBoundary = tavgBoundary();
+    labelList nBoundary(pointFaces.size(), 0);
+
+    forAll(pointFaces, pointI)
+    {
+        const labelList& pFaces = pointFaces[pointI];
+
+        forAll(pFaces, pfI)
+        {
+            label faceI = pFaces[pfI];
+            label meshFaceI = pp.addressing()[faceI];
+
+            label own = mesh.faceOwner()[meshFaceI];
+            avgBoundary[pointI] += mesh.cellCentres()[own];
+            nBoundary[pointI]++;
+        }
+    }
+
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        avgBoundary,
+        plusEqOp<point>(),  // combine op
+        vector::zero        // null value
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        nBoundary,
+        plusEqOp<label>(),  // combine op
+        label(0)            // null value
+    );
+
+    forAll(avgBoundary, i)
+    {
+        avgBoundary[i] /= nBoundary[i];
+    }
+    return tavgBoundary;
+}
+
+
+//Foam::tmp<Foam::scalarField> Foam::autoSnapDriver::calcEdgeLen
+//(
+//    const indirectPrimitivePatch& pp
+//) const
+//{
+//    // Get local edge length based on refinement level
+//    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+//    // (Ripped from autoLayerDriver)
+//
+//    tmp<scalarField> tedgeLen(new scalarField(pp.nPoints()));
+//    scalarField& edgeLen = tedgeLen();
+//    {
+//        const fvMesh& mesh = meshRefiner_.mesh();
+//        const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
+//        const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
+//
+//        labelList maxPointLevel(pp.nPoints(), labelMin);
+//
+//        forAll(pp, i)
+//        {
+//            label ownLevel = cellLevel[mesh.faceOwner()[pp.addressing()[i]]];
+//            const face& f = pp.localFaces()[i];
+//            forAll(f, fp)
+//            {
+//                maxPointLevel[f[fp]] = max(maxPointLevel[f[fp]], ownLevel);
+//            }
+//        }
+//
+//        syncTools::syncPointList
+//        (
+//            mesh,
+//            pp.meshPoints(),
+//            maxPointLevel,
+//            maxEqOp<label>(),
+//            labelMin            // null value
+//        );
+//
+//
+//        forAll(maxPointLevel, pointI)
+//        {
+//            // Find undistorted edge size for this level.
+//            edgeLen[pointI] = edge0Len/(1<<maxPointLevel[pointI]);
+//        }
+//    }
+//    return tedgeLen;
+//}
+
+
+void Foam::autoSnapDriver::detectNearSurfaces
+(
+    const scalar planarCos,
+    const indirectPrimitivePatch& pp,
+    const pointField& nearestPoint,
+    const vectorField& nearestNormal,
+
+    vectorField& disp
+) const
+{
+    Info<< "Detecting near surfaces ..." << endl;
+
+    const pointField& localPoints = pp.localPoints();
+    const refinementSurfaces& surfaces = meshRefiner_.surfaces();
+    const fvMesh& mesh = meshRefiner_.mesh();
+
+    //// Get local edge length based on refinement level
+    //const scalarField edgeLen(calcEdgeLen(pp));
+    //
+    //// Generate rays for every surface point
+    //// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    //
+    //{
+    //    const scalar cos45 = Foam::cos(degToRad(45));
+    //    vector n(cos45, cos45, cos45);
+    //    n /= mag(n);
+    //
+    //    pointField start(14*pp.nPoints());
+    //    pointField end(start.size());
+    //
+    //    label rayI = 0;
+    //    forAll(localPoints, pointI)
+    //    {
+    //        const point& pt = localPoints[pointI];
+    //
+    //        // Along coordinate axes
+    //
+    //        {
+    //            start[rayI] = pt;
+    //            point& endPt = end[rayI++];
+    //            endPt = pt;
+    //            endPt.x() -= edgeLen[pointI];
+    //        }
+    //        {
+    //            start[rayI] = pt;
+    //            point& endPt = end[rayI++];
+    //            endPt = pt;
+    //            endPt.x() += edgeLen[pointI];
+    //        }
+    //        {
+    //            start[rayI] = pt;
+    //            point& endPt = end[rayI++];
+    //            endPt = pt;
+    //            endPt.y() -= edgeLen[pointI];
+    //        }
+    //        {
+    //            start[rayI] = pt;
+    //            point& endPt = end[rayI++];
+    //            endPt = pt;
+    //            endPt.y() += edgeLen[pointI];
+    //        }
+    //        {
+    //            start[rayI] = pt;
+    //            point& endPt = end[rayI++];
+    //            endPt = pt;
+    //            endPt.z() -= edgeLen[pointI];
+    //        }
+    //        {
+    //            start[rayI] = pt;
+    //            point& endPt = end[rayI++];
+    //            endPt = pt;
+    //            endPt.z() += edgeLen[pointI];
+    //        }
+    //
+    //        // At 45 degrees
+    //
+    //        const vector vec(edgeLen[pointI]*n);
+    //
+    //        {
+    //            start[rayI] = pt;
+    //            point& endPt = end[rayI++];
+    //            endPt = pt;
+    //            endPt.x() += vec.x();
+    //            endPt.y() += vec.y();
+    //            endPt.z() += vec.z();
+    //        }
+    //        {
+    //            start[rayI] = pt;
+    //            point& endPt = end[rayI++];
+    //            endPt = pt;
+    //            endPt.x() -= vec.x();
+    //            endPt.y() += vec.y();
+    //            endPt.z() += vec.z();
+    //        }
+    //        {
+    //            start[rayI] = pt;
+    //            point& endPt = end[rayI++];
+    //            endPt = pt;
+    //            endPt.x() += vec.x();
+    //            endPt.y() -= vec.y();
+    //            endPt.z() += vec.z();
+    //        }
+    //        {
+    //            start[rayI] = pt;
+    //            point& endPt = end[rayI++];
+    //            endPt = pt;
+    //            endPt.x() -= vec.x();
+    //            endPt.y() -= vec.y();
+    //            endPt.z() += vec.z();
+    //        }
+    //        {
+    //            start[rayI] = pt;
+    //            point& endPt = end[rayI++];
+    //            endPt = pt;
+    //            endPt.x() += vec.x();
+    //            endPt.y() += vec.y();
+    //            endPt.z() -= vec.z();
+    //        }
+    //        {
+    //            start[rayI] = pt;
+    //            point& endPt = end[rayI++];
+    //            endPt = pt;
+    //            endPt.x() -= vec.x();
+    //            endPt.y() += vec.y();
+    //            endPt.z() -= vec.z();
+    //        }
+    //        {
+    //            start[rayI] = pt;
+    //            point& endPt = end[rayI++];
+    //            endPt = pt;
+    //            endPt.x() += vec.x();
+    //            endPt.y() -= vec.y();
+    //            endPt.z() -= vec.z();
+    //        }
+    //        {
+    //            start[rayI] = pt;
+    //            point& endPt = end[rayI++];
+    //            endPt = pt;
+    //            endPt.x() -= vec.x();
+    //            endPt.y() -= vec.y();
+    //            endPt.z() -= vec.z();
+    //        }
+    //    }
+    //
+    //    labelList surface1;
+    //    List<pointIndexHit> hit1;
+    //    labelList region1;
+    //    vectorField normal1;
+    //
+    //    labelList surface2;
+    //    List<pointIndexHit> hit2;
+    //    labelList region2;
+    //    vectorField normal2;
+    //    surfaces.findNearestIntersection
+    //    (
+    //        unzonedSurfaces,    // surfacesToTest,
+    //        start,
+    //        end,
+    //
+    //        surface1,
+    //        hit1,
+    //        region1,
+    //        normal1,
+    //
+    //        surface2,
+    //        hit2,
+    //        region2,
+    //        normal2
+    //    );
+    //
+    //    // All intersections
+    //    {
+    //        OBJstream str
+    //        (
+    //            mesh.time().path()
+    //          / "surfaceHits_" + meshRefiner_.timeName() + ".obj"
+    //        );
+    //
+    //        Info<< "Dumping intersections with rays to " << str.name()
+    //            << endl;
+    //
+    //        forAll(hit1, i)
+    //        {
+    //            if (hit1[i].hit())
+    //            {
+    //                str.write(linePointRef(start[i], hit1[i].hitPoint()));
+    //            }
+    //            if (hit2[i].hit())
+    //            {
+    //                str.write(linePointRef(start[i], hit2[i].hitPoint()));
+    //            }
+    //        }
+    //    }
+    //
+    //    // Co-planar intersections
+    //    {
+    //        OBJstream str
+    //        (
+    //            mesh.time().path()
+    //          / "coplanarHits_" + meshRefiner_.timeName() + ".obj"
+    //        );
+    //
+    //        Info<< "Dumping intersections with co-planar surfaces to "
+    //            << str.name() << endl;
+    //
+    //        forAll(localPoints, pointI)
+    //        {
+    //            bool hasNormal = false;
+    //            point surfPointA;
+    //            vector surfNormalA;
+    //            point surfPointB;
+    //            vector surfNormalB;
+    //
+    //            bool isCoplanar = false;
+    //
+    //            label rayI = 14*pointI;
+    //            for (label i = 0; i < 14; i++)
+    //            {
+    //                if (hit1[rayI].hit())
+    //                {
+    //                    const point& pt = hit1[rayI].hitPoint();
+    //                    const vector& n = normal1[rayI];
+    //
+    //                    if (!hasNormal)
+    //                    {
+    //                        hasNormal = true;
+    //                        surfPointA = pt;
+    //                        surfNormalA = n;
+    //                    }
+    //                    else
+    //                    {
+    //                        if
+    //                        (
+    //                            meshRefiner_.isGap
+    //                            (
+    //                                planarCos,
+    //                                surfPointA,
+    //                                surfNormalA,
+    //                                pt,
+    //                                n
+    //                            )
+    //                        )
+    //                        {
+    //                            isCoplanar = true;
+    //                            surfPointB = pt;
+    //                            surfNormalB = n;
+    //                            break;
+    //                        }
+    //                    }
+    //                }
+    //                if (hit2[rayI].hit())
+    //                {
+    //                    const point& pt = hit2[rayI].hitPoint();
+    //                    const vector& n = normal2[rayI];
+    //
+    //                    if (!hasNormal)
+    //                    {
+    //                        hasNormal = true;
+    //                        surfPointA = pt;
+    //                        surfNormalA = n;
+    //                    }
+    //                    else
+    //                    {
+    //                        if
+    //                        (
+    //                            meshRefiner_.isGap
+    //                            (
+    //                                planarCos,
+    //                                surfPointA,
+    //                                surfNormalA,
+    //                                pt,
+    //                                n
+    //                            )
+    //                        )
+    //                        {
+    //                            isCoplanar = true;
+    //                            surfPointB = pt;
+    //                            surfNormalB = n;
+    //                            break;
+    //                        }
+    //                    }
+    //                }
+    //
+    //                rayI++;
+    //            }
+    //
+    //            if (isCoplanar)
+    //            {
+    //                str.write(linePointRef(surfPointA, surfPointB));
+    //            }
+    //        }
+    //    }
+    //}
+
+
+    const pointField avgCc(avgCellCentres(mesh, pp));
+
+    // Construct rays from localPoints to beyond cell centre
+    pointField end(pp.nPoints());
+    forAll(localPoints, pointI)
+    {
+        const point& pt = localPoints[pointI];
+        end[pointI] = pt + 2*(avgCc[pointI]-pt);
+    }
+
+
+    label nOverride = 0;
+
+    // 1. All points to non-interface surfaces
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    {
+        const labelList unzonedSurfaces =
+            surfaceZonesInfo::getUnnamedSurfaces
+            (
+                meshRefiner_.surfaces().surfZones()
+            );
+
+        // Do intersection test
+        labelList surface1;
+        List<pointIndexHit> hit1;
+        labelList region1;
+        vectorField normal1;
+
+        labelList surface2;
+        List<pointIndexHit> hit2;
+        labelList region2;
+        vectorField normal2;
+        surfaces.findNearestIntersection
+        (
+            unzonedSurfaces,
+            localPoints,
+            end,
+
+            surface1,
+            hit1,
+            region1,
+            normal1,
+
+            surface2,
+            hit2,
+            region2,
+            normal2
+        );
+
+
+        forAll(localPoints, pointI)
+        {
+            // Current location
+            const point& pt = localPoints[pointI];
+
+            bool override = false;
+
+            if (hit1[pointI].hit())
+            {
+                if
+                (
+                    meshRefiner_.isGap
+                    (
+                        planarCos,
+                        nearestPoint[pointI],
+                        nearestNormal[pointI],
+                        hit1[pointI].hitPoint(),
+                        normal1[pointI]
+                    )
+                )
+                {
+                    disp[pointI] = hit1[pointI].hitPoint()-pt;
+                    override = true;
+                }
+            }
+            if (hit2[pointI].hit())
+            {
+                if
+                (
+                    meshRefiner_.isGap
+                    (
+                        planarCos,
+                        nearestPoint[pointI],
+                        nearestNormal[pointI],
+                        hit2[pointI].hitPoint(),
+                        normal2[pointI]
+                    )
+                )
+                {
+                    disp[pointI] = hit2[pointI].hitPoint()-pt;
+                    override = true;
+                }
+            }
+
+            if (override)
+            {
+                nOverride++;
+            }
+        }
+    }
+
+
+    // 2. All points on zones to their respective surface
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    {
+        // Surfaces with zone information
+        const PtrList<surfaceZonesInfo>& surfZones = surfaces.surfZones();
+
+        const labelList zonedSurfaces = surfaceZonesInfo::getNamedSurfaces
+        (
+            surfZones
+        );
+
+        forAll(zonedSurfaces, i)
+        {
+            label zoneSurfI = zonedSurfaces[i];
+
+            const word& faceZoneName = surfZones[zoneSurfI].faceZoneName();
+
+            const labelList surfacesToTest(1, zoneSurfI);
+
+            // Get indices of points both on faceZone and on pp.
+            labelList zonePointIndices
+            (
+                getZoneSurfacePoints
+                (
+                    mesh,
+                    pp,
+                    faceZoneName
+                )
+            );
+
+            // Do intersection test
+            labelList surface1;
+            List<pointIndexHit> hit1;
+            labelList region1;
+            vectorField normal1;
+
+            labelList surface2;
+            List<pointIndexHit> hit2;
+            labelList region2;
+            vectorField normal2;
+            surfaces.findNearestIntersection
+            (
+                surfacesToTest,
+                pointField(localPoints, zonePointIndices),
+                pointField(end, zonePointIndices),
+
+                surface1,
+                hit1,
+                region1,
+                normal1,
+
+                surface2,
+                hit2,
+                region2,
+                normal2
+            );
+
+
+            label nOverride = 0;
+
+            forAll(hit1, i)
+            {
+                label pointI = zonePointIndices[i];
+
+                // Current location
+                const point& pt = localPoints[pointI];
+
+                bool override = false;
+
+                if (hit1[i].hit())
+                {
+                    if
+                    (
+                        meshRefiner_.isGap
+                        (
+                            planarCos,
+                            nearestPoint[pointI],
+                            nearestNormal[pointI],
+                            hit1[i].hitPoint(),
+                            normal1[i]
+                        )
+                    )
+                    {
+                        disp[pointI] = hit1[i].hitPoint()-pt;
+                        override = true;
+                    }
+                }
+                if (hit2[i].hit())
+                {
+                    if
+                    (
+                        meshRefiner_.isGap
+                        (
+                            planarCos,
+                            nearestPoint[pointI],
+                            nearestNormal[pointI],
+                            hit2[i].hitPoint(),
+                            normal2[i]
+                        )
+                    )
+                    {
+                        disp[pointI] = hit2[i].hitPoint()-pt;
+                        override = true;
+                    }
+                }
+
+                if (override)
+                {
+                    nOverride++;
+                }
+            }
+        }
+    }
+
+    Info<< "Overriding nearest with intersection of close gaps at "
+        << returnReduce(nOverride, sumOp<label>())
+        << " out of " << returnReduce(pp.nPoints(), sumOp<label>())
+        << " points." << endl;
+}
+
+
 Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
 (
     const meshRefinement& meshRefiner,
     const scalarField& snapDist,
-    const indirectPrimitivePatch& pp
+    const indirectPrimitivePatch& pp,
+    pointField& nearestPoint,
+    vectorField& nearestNormal
 )
 {
     Info<< "Calculating patchDisplacement as distance to nearest surface"
@@ -815,12 +1440,12 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
         labelList snapSurf(localPoints.size(), -1);
 
         // Divide surfaces into zoned and unzoned
-        labelList zonedSurfaces =
+        const labelList zonedSurfaces =
             surfaceZonesInfo::getNamedSurfaces
             (
                 meshRefiner.surfaces().surfZones()
             );
-        labelList unzonedSurfaces =
+        const labelList unzonedSurfaces =
             surfaceZonesInfo::getUnnamedSurfaces
             (
                 meshRefiner.surfaces().surfZones()
@@ -833,14 +1458,42 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
         {
             List<pointIndexHit> hitInfo;
             labelList hitSurface;
-            surfaces.findNearest
-            (
-                unzonedSurfaces,
-                localPoints,
-                sqr(snapDist),        // sqr of attract distance
-                hitSurface,
-                hitInfo
-            );
+
+            if (nearestNormal.size() == localPoints.size())
+            {
+                labelList hitRegion;
+                vectorField hitNormal;
+                surfaces.findNearestRegion
+                (
+                    unzonedSurfaces,
+                    localPoints,
+                    sqr(snapDist),
+                    hitSurface,
+                    hitInfo,
+                    hitRegion,
+                    hitNormal
+                );
+
+                forAll(hitInfo, pointI)
+                {
+                    if (hitInfo[pointI].hit())
+                    {
+                        nearestPoint[pointI] = hitInfo[pointI].hitPoint();
+                        nearestNormal[pointI] = hitNormal[pointI];
+                    }
+                }
+            }
+            else
+            {
+                surfaces.findNearest
+                (
+                    unzonedSurfaces,
+                    localPoints,
+                    sqr(snapDist),        // sqr of attract distance
+                    hitSurface,
+                    hitInfo
+                );
+            }
 
             forAll(hitInfo, pointI)
             {
@@ -888,14 +1541,43 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
             // Find nearest for points both on faceZone and pp.
             List<pointIndexHit> hitInfo;
             labelList hitSurface;
-            surfaces.findNearest
-            (
-                labelList(1, zoneSurfI),
-                pointField(localPoints, zonePointIndices),
-                sqr(scalarField(minSnapDist, zonePointIndices)),
-                hitSurface,
-                hitInfo
-            );
+
+            if (nearestNormal.size() == localPoints.size())
+            {
+                labelList hitRegion;
+                vectorField hitNormal;
+                surfaces.findNearestRegion
+                (
+                    surfacesToTest,
+                    pointField(localPoints, zonePointIndices),
+                    sqr(scalarField(minSnapDist, zonePointIndices)),
+                    hitSurface,
+                    hitInfo,
+                    hitRegion,
+                    hitNormal
+                );
+
+                forAll(hitInfo, i)
+                {
+                    if (hitInfo[i].hit())
+                    {
+                        label pointI = zonePointIndices[i];
+                        nearestPoint[pointI] = hitInfo[i].hitPoint();
+                        nearestNormal[pointI] = hitNormal[i];
+                    }
+                }
+            }
+            else
+            {
+                surfaces.findNearest
+                (
+                    surfacesToTest,
+                    pointField(localPoints, zonePointIndices),
+                    sqr(scalarField(minSnapDist, zonePointIndices)),
+                    hitSurface,
+                    hitInfo
+                );
+            }
 
             forAll(hitInfo, i)
             {
@@ -1529,6 +2211,7 @@ void Foam::autoSnapDriver::doSnap
     const dictionary& snapDict,
     const dictionary& motionDict,
     const scalar featureCos,
+    const scalar planarAngle,
     const snapParameters& snapParams
 )
 {
@@ -1791,13 +2474,40 @@ void Foam::autoSnapDriver::doSnap
 
             // Calculate displacement at every patch point. Insert into
             // meshMover.
+            // Calculate displacement at every patch point
+            pointField nearestPoint;
+            vectorField nearestNormal;
+
+            if (snapParams.detectNearSurfacesSnap())
+            {
+                nearestPoint.setSize(pp.nPoints(), vector::max);
+                nearestNormal.setSize(pp.nPoints(), vector::zero);
+            }
+
             vectorField disp = calcNearestSurface
             (
                 meshRefiner_,
                 snapDist,
-                meshMover.patch()
+                pp,
+                nearestPoint,
+                nearestNormal
             );
 
+
+            // Override displacement at thin gaps
+            if (snapParams.detectNearSurfacesSnap())
+            {
+                detectNearSurfaces
+                (
+                    Foam::cos(degToRad(planarAngle)),// planar cos for gaps
+                    pp,
+                    nearestPoint,   // surfacepoint from nearest test
+                    nearestNormal,  // surfacenormal from nearest test
+
+                    disp
+                );
+            }
+
             // Override displacement with feature edge attempt
             if (doFeatures)
             {
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index c0a0adcde2..77327c0689 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -148,6 +148,7 @@ class autoSnapDriver
                 (
                     const label iter,
                     const indirectPrimitivePatch& pp,
+                    const scalarField& faceSnapDist,
                     vectorField& faceDisp,
                     vectorField& faceSurfaceNormal,
                     labelList& faceSurfaceRegion,
@@ -432,6 +433,23 @@ public:
                 const word& zoneName
             );
 
+            //- Helper: calculate average cell centre per point
+            static tmp<pointField> avgCellCentres
+            (
+                const fvMesh& mesh,
+                const indirectPrimitivePatch&
+            );
+
+            //- Per patch point override displacement if in gap situation
+            void detectNearSurfaces
+            (
+                const scalar planarCos,
+                const indirectPrimitivePatch&,
+                const pointField& nearestPoint,
+                const vectorField& nearestNormal,
+                vectorField& disp
+            ) const;
+
             //- Per patch point calculate point on nearest surface. Set as
             //  boundary conditions of motionSmoother displacement field. Return
             //  displacement of patch points.
@@ -439,7 +457,9 @@ public:
             (
                 const meshRefinement& meshRefiner,
                 const scalarField& snapDist,
-                const indirectPrimitivePatch&
+                const indirectPrimitivePatch&,
+                pointField& nearestPoint,
+                vectorField& nearestNormal
             );
 
             ////- Per patch point calculate point on nearest surface. Set as
@@ -483,6 +503,7 @@ public:
                 const dictionary& snapDict,
                 const dictionary& motionDict,
                 const scalar featureCos,
+                const scalar planarAngle,
                 const snapParameters& snapParams
             );
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index ad2ad93a91..887674ea7a 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -319,6 +319,7 @@ void Foam::autoSnapDriver::calcNearestFace
 (
     const label iter,
     const indirectPrimitivePatch& pp,
+    const scalarField& faceSnapDist,
     vectorField& faceDisp,
     vectorField& faceSurfaceNormal,
     labelList& faceSurfaceGlobalRegion,
@@ -409,7 +410,7 @@ void Foam::autoSnapDriver::calcNearestFace
         (
             labelList(1, zoneSurfI),
             fc,
-            sqr(scalarField(fc.size(), GREAT)),// sqr of attract dist
+            sqr(faceSnapDist),// sqr of attract dist
             hitSurface,
             hitInfo,
             hitRegion,
@@ -429,6 +430,14 @@ void Foam::autoSnapDriver::calcNearestFace
                     hitRegion[hitI]
                 );
             }
+            else
+            {
+                WarningIn
+                (
+                    "autoSnapDriver::calcNearestFace(..)"
+                )   << "Did not find surface near face centre " << fc[hitI]
+                    << endl;
+            }
         }
     }
 
@@ -467,7 +476,7 @@ void Foam::autoSnapDriver::calcNearestFace
     (
         unzonedSurfaces,
         fc,
-        sqr(scalarField(fc.size(), GREAT)),// sqr of attract dist
+        sqr(faceSnapDist),// sqr of attract dist
         hitSurface,
         hitInfo,
         hitRegion,
@@ -487,6 +496,14 @@ void Foam::autoSnapDriver::calcNearestFace
                 hitRegion[hitI]
             );
         }
+        else
+        {
+            WarningIn
+            (
+                "autoSnapDriver::calcNearestFace(..)"
+            )   << "Did not find surface near face centre " << fc[hitI]
+                << endl;
+        }
     }
 
 
@@ -560,24 +577,44 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
     forAll(pp.pointFaces(), pointI)
     {
         const labelList& pFaces = pp.pointFaces()[pointI];
-        List<point>& pNormals = pointFaceSurfNormals[pointI];
-        pNormals.setSize(pFaces.size());
-        List<point>& pDisp = pointFaceDisp[pointI];
-        pDisp.setSize(pFaces.size());
-        List<point>& pFc = pointFaceCentres[pointI];
-        pFc.setSize(pFaces.size());
-        labelList& pFid = pointFacePatchID[pointI];
-        pFid.setSize(pFaces.size());
 
+        // Count valid face normals
+        label nFaces = 0;
         forAll(pFaces, i)
         {
             label faceI = pFaces[i];
-            pNormals[i] = faceSurfaceNormal[faceI];
-            pDisp[i] = faceDisp[faceI];
-            pFc[i] = pp.faceCentres()[faceI];
-            //label meshFaceI = pp.addressing()[faceI];
-            //pFid[i] = mesh.boundaryMesh().whichPatch(meshFaceI);
-            pFid[i] = globalToMasterPatch_[faceSurfaceGlobalRegion[faceI]];
+            if (faceSurfaceGlobalRegion[faceI] != -1)
+            {
+                nFaces++;
+            }
+        }
+
+
+        List<point>& pNormals = pointFaceSurfNormals[pointI];
+        pNormals.setSize(nFaces);
+        List<point>& pDisp = pointFaceDisp[pointI];
+        pDisp.setSize(nFaces);
+        List<point>& pFc = pointFaceCentres[pointI];
+        pFc.setSize(nFaces);
+        labelList& pFid = pointFacePatchID[pointI];
+        pFid.setSize(nFaces);
+
+        nFaces = 0;
+        forAll(pFaces, i)
+        {
+            label faceI = pFaces[i];
+            label globalRegionI = faceSurfaceGlobalRegion[faceI];
+
+            if (globalRegionI != -1)
+            {
+                pNormals[nFaces] = faceSurfaceNormal[faceI];
+                pDisp[nFaces] = faceDisp[faceI];
+                pFc[nFaces] = pp.faceCentres()[faceI];
+                //label meshFaceI = pp.addressing()[faceI];
+                //pFid[nFaces] = mesh.boundaryMesh().whichPatch(meshFaceI);
+                pFid[nFaces] = globalToMasterPatch_[globalRegionI];
+                nFaces++;
+            }
         }
     }
 
@@ -2123,9 +2160,85 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
         bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
     }
 
+    // Collect candidate points for attraction
+    DynamicList<label> attractPoints(pp.nPoints());
+    {
+        const fvMesh& mesh = meshRefiner_.mesh();
+
+        boolList isFeatureEdgeOrPoint(pp.nPoints(), false);
+        label nFeats = 0;
+        forAll(rawPatchConstraints, pointI)
+        {
+            if (rawPatchConstraints[pointI].first() >= 2)
+            {
+                isFeatureEdgeOrPoint[pointI] = true;
+                nFeats++;
+            }
+        }
+
+        Info<< "Intially selected " << returnReduce(nFeats, sumOp<label>())
+            << " points out of " << returnReduce(pp.nPoints(), sumOp<label>())
+            << " for reverse attraction." << endl;
+
+        // Make sure is synchronised (note: check if constraint is already
+        // synced in which case this is not needed here)
+        syncTools::syncPointList
+        (
+            mesh,
+            pp.meshPoints(),
+            isFeatureEdgeOrPoint,
+            orEqOp<bool>(),         // combine op
+            false
+        );
+
+        for (label nGrow = 0; nGrow < 1; nGrow++)
+        {
+            forAll(pp.localFaces(), faceI)
+            {
+                const face& f = pp.localFaces()[faceI];
+
+                forAll(f, fp)
+                {
+                    if (isFeatureEdgeOrPoint[f[fp]])
+                    {
+                        // Mark all points on face
+                        forAll(f, fp)
+                        {
+                            isFeatureEdgeOrPoint[f[fp]] = true;
+                        }
+                        break;
+                    }
+                }
+            }
+            syncTools::syncPointList
+            (
+                mesh,
+                pp.meshPoints(),
+                isFeatureEdgeOrPoint,
+                orEqOp<bool>(),         // combine op
+                false
+            );
+        }
+
+
+        // Collect attractPoints
+        forAll(isFeatureEdgeOrPoint, pointI)
+        {
+            if (isFeatureEdgeOrPoint[pointI])
+            {
+                attractPoints.append(pointI);
+            }
+        }
+
+        Info<< "Selected " << returnReduce(attractPoints.size(), sumOp<label>())
+            << " points out of " << returnReduce(pp.nPoints(), sumOp<label>())
+            << " for reverse attraction." << endl;
+    }
+
+
     indexedOctree<treeDataPoint> ppTree
     (
-        treeDataPoint(pp.localPoints()),
+        treeDataPoint(pp.localPoints(), attractPoints),
         bb,                             // overall search domain
         8,                              // maxLevel
         10,                             // leafsize
@@ -2159,7 +2272,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
 
                 if (nearInfo.hit())
                 {
-                    label pointI = nearInfo.index();
+                    label pointI =
+                        ppTree.shapes().pointLabels()[nearInfo.index()];
                     const point attraction = featPt-pp.localPoints()[pointI];
 
                     // Check if this point is already being attracted. If so
@@ -2220,7 +2334,9 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
 
                 if (nearInfo.hit())
                 {
-                    label pointI = nearInfo.index();
+                    label pointI =
+                        ppTree.shapes().pointLabels()[nearInfo.index()];
+
                     const point& pt = pp.localPoints()[pointI];
                     const point attraction = featPt-pt;
 
@@ -2681,6 +2797,19 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     const fvMesh& mesh = meshRefiner_.mesh();
 
 
+    // Calculate attraction distance per face (from the attraction distance
+    // per point)
+    scalarField faceSnapDist(pp.size(), -GREAT);
+    forAll(pp.localFaces(), faceI)
+    {
+        const face& f = pp.localFaces()[faceI];
+        forAll(f, fp)
+        {
+            faceSnapDist[faceI] = max(faceSnapDist[faceI], snapDist[f[fp]]);
+        }
+    }
+
+
     // Displacement and orientation per pp face
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -2695,6 +2824,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     (
         iter,
         pp,
+        faceSnapDist,
         faceDisp,
         faceSurfaceNormal,
         faceSurfaceGlobalRegion,
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C
index 77bff60187..364a6aae49 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-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,11 +40,12 @@ Foam::snapParameters::snapParameters(const dictionary& dict)
     multiRegionFeatureSnap_
     (
         dict.lookupOrDefault("multiRegionFeatureSnap", false)
+    ),
+    detectNearSurfacesSnap_
+    (
+        dict.lookupOrDefault("detectNearSurfacesSnap", true)
     )
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H
index 294d0de55d..5dc4ca34e8 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-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -70,6 +70,8 @@ class snapParameters
 
         const Switch multiRegionFeatureSnap_;
 
+        const Switch detectNearSurfacesSnap_;
+
 
     // Private Member Functions
 
@@ -141,6 +143,11 @@ public:
                 return multiRegionFeatureSnap_;
             }
 
+            Switch detectNearSurfacesSnap() const
+            {
+                return detectNearSurfacesSnap_;
+            }
+
 };
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 499c6e6c0a..b07f64b51a 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -289,16 +289,6 @@ private:
                 label& nRefine
             ) const;
 
-            //- Is local topology a small gap?
-            bool isGap
-            (
-                const scalar,
-                const vector&,
-                const vector&,
-                const vector&,
-                const vector&
-            ) const;
-
             //- Mark cell if local a gap topology or
             bool checkProximity
             (
@@ -664,6 +654,26 @@ public:
 
         // Refinement
 
+            //- Is local topology a small gap?
+            bool isGap
+            (
+                const scalar,
+                const vector&,
+                const vector&,
+                const vector&,
+                const vector&
+            ) const;
+
+            //- Is local topology a small gap normal to the test vector
+            bool isNormalGap
+            (
+                const scalar,
+                const vector&,
+                const vector&,
+                const vector&,
+                const vector&
+            ) const;
+
             //- Calculate list of cells to refine.
             labelList refineCandidates
             (
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
index 8d7e38512b..a1c1ed0f58 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
@@ -1061,13 +1061,17 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
             meshMover
         );
 
+        pointField nearestPoint;
+        vectorField nearestNormal;
         const vectorField disp
         (
             autoSnapDriver::calcNearestSurface
             (
                 *this,
                 snapDist,   // attraction
-                pp
+                pp,
+                nearestPoint,
+                nearestNormal
             )
         );
 
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index 9667750cf5..99597db78d 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -1214,7 +1214,6 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
 }
 
 
-// Mark small gaps
 bool Foam::meshRefinement::isGap
 (
     const scalar planarCos,
@@ -1225,12 +1224,53 @@ bool Foam::meshRefinement::isGap
     const vector& normal1
 ) const
 {
-    ////- hits differ and angles are oppositeish
-    //return
-    //    (mag(point0-point1) > mergeDistance())
-    // && ((normal0 & normal1) < (-1+planarCos));
+    //- hits differ and angles are oppositeish and
+    //  hits have a normal distance
+    vector d = point1-point0;
+    scalar magD = mag(d);
+
+    if (magD > mergeDistance())
+    {
+        scalar cosAngle = (normal0 & normal1);
+
+        vector avg = vector::zero;
+        if (cosAngle < (-1+planarCos))
+        {
+            // Opposite normals
+            avg = 0.5*(normal0-normal1);
+        }
+        else if (cosAngle > (1-planarCos))
+        {
+            avg = 0.5*(normal0+normal1);
+        }
+
+        if (avg != vector::zero)
+        {
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+    else
+    {
+        return false;
+    }
+}
 
 
+// Mark small gaps
+bool Foam::meshRefinement::isNormalGap
+(
+    const scalar planarCos,
+    const vector& point0,
+    const vector& normal0,
+
+    const vector& point1,
+    const vector& normal1
+) const
+{
     //- hits differ and angles are oppositeish and
     //  hits have a normal distance
     vector d = point1-point0;
@@ -1316,7 +1356,7 @@ bool Foam::meshRefinement::checkProximity
             //  - different location
             //  - opposite surface
 
-            bool closeSurfaces = isGap
+            bool closeSurfaces = isNormalGap
             (
                 planarCos,
                 cellMaxLocation,
@@ -1565,7 +1605,7 @@ Foam::label Foam::meshRefinement::markProximityRefinement
             // Have valid data on both sides. Check planarCos.
             if
             (
-                isGap
+                isNormalGap
                 (
                     planarCos,
                     cellMaxLocation[own],
@@ -1635,7 +1675,7 @@ Foam::label Foam::meshRefinement::markProximityRefinement
             // Have valid data on both sides. Check planarCos.
             if
             (
-                isGap
+                isNormalGap
                 (
                     planarCos,
                     cellMaxLocation[own],