diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index e358ba4ffa4d04a5ea96f1571d9286feb790ab31..54a78f477c6a5e7170572da5c8867654c564c543 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -1406,9 +1406,10 @@ void Foam::autoSnapDriver::doSnap
                 adaptPatchIDs
             )
         );
+        indirectPrimitivePatch& pp = ppPtr();
 
         // Distance to attract to nearest feature on surface
-        const scalarField snapDist(calcSnapDistance(snapParams, ppPtr()));
+        const scalarField snapDist(calcSnapDistance(snapParams, pp));
 
 
         // Construct iterative mesh mover.
@@ -1420,7 +1421,7 @@ void Foam::autoSnapDriver::doSnap
         motionSmoother meshMover
         (
             mesh,
-            ppPtr(),
+            pp,
             adaptPatchIDs,
             meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs),
             motionDict
@@ -1475,7 +1476,7 @@ void Foam::autoSnapDriver::doSnap
             }
 
             // Check for displacement being outwards.
-            outwardsDisplacement(ppPtr(), disp);
+            outwardsDisplacement(pp, disp);
 
             // Set initial distribution of displacement field (on patches)
             // from patchDisp and make displacement consistent with b.c.
@@ -1489,8 +1490,8 @@ void Foam::autoSnapDriver::doSnap
                 (
                     mesh.time().path()
                   / "patchDisplacement_" + name(iter) + ".obj",
-                    ppPtr().localPoints(),
-                    ppPtr().localPoints() + disp
+                    pp.localPoints(),
+                    pp.localPoints() + disp
                 );
             }
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index bd5dd86188fa9e64b467d600257b569679e63a77..7cf1e79f6121a963dc5ccdf108ac05e0bd7f75d9 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -141,14 +141,6 @@ class autoSnapDriver
                     vectorField& pointSurfaceNormal,
                     vectorField& pointRotation
                 ) const;
-                //void calcNearestFace
-                //(
-                //    const label iter,
-                //    const indirectPrimitivePatch& pp,
-                //    vectorField& faceDisp,
-                //    vectorField& faceSurfaceNormal,
-                //    vectorField& faceRotation
-                //) const;
                 void calcNearestFace
                 (
                     const label iter,
@@ -158,15 +150,19 @@ class autoSnapDriver
                     labelList& faceSurfaceRegion,
                     vectorField& faceRotation
                 ) const;
-                void interpolateFaceToPoint
+                void calcNearestFacePointProperties
                 (
                     const label iter,
                     const indirectPrimitivePatch& pp,
-                    const List<List<point> >& pointFaceDisp,
-                    const List<List<point> >& pointFaceRotation,
-                    const List<List<point> >& pointFaceCentres,
-                    vectorField& patchDisp
-                    //vectorField& patchRotationDisp
+
+                    const vectorField& faceDisp,
+                    const vectorField& faceSurfaceNormal,
+                    const labelList& faceSurfaceRegion,
+
+                    List<List<point> >& pointFaceSurfNormals,
+                    List<List<point> >& pointFaceDisp,
+                    List<List<point> >& pointFaceCentres,
+                    List<labelList>&    pointFacePatchID
                 ) const;
                 void correctAttraction
                 (
@@ -276,6 +272,7 @@ class autoSnapDriver
                 (
                     const label iter,
                     const scalar featureCos,
+                    const bool multiRegionFeatureSnap,
 
                     const indirectPrimitivePatch&,
                     const scalarField&,
@@ -339,6 +336,7 @@ class autoSnapDriver
                 (
                     const label iter,
                     const scalar featureCos,
+                    const bool multiRegionFeatureSnap,
                     const indirectPrimitivePatch& pp,
                     const scalarField& snapDist,
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index 8bcd9c6270a22ae42ce171a63d8f50000919204e..f32d9ce4bd3a775e76bac478a361c1a96d6aab16 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -321,7 +321,7 @@ void Foam::autoSnapDriver::calcNearestFace
     const indirectPrimitivePatch& pp,
     vectorField& faceDisp,
     vectorField& faceSurfaceNormal,
-    labelList& faceSurfaceRegion,
+    labelList& faceSurfaceGlobalRegion,
     vectorField& faceRotation
 ) const
 {
@@ -333,8 +333,8 @@ void Foam::autoSnapDriver::calcNearestFace
     faceDisp = vector::zero;
     faceSurfaceNormal.setSize(pp.size());
     faceSurfaceNormal = vector::zero;
-    faceSurfaceRegion.setSize(pp.size());
-    faceSurfaceRegion = -1;
+    faceSurfaceGlobalRegion.setSize(pp.size());
+    faceSurfaceGlobalRegion = -1;
 
     // Divide surfaces into zoned and unzoned
     labelList zonedSurfaces = surfaces.getNamedSurfaces();
@@ -419,7 +419,7 @@ void Foam::autoSnapDriver::calcNearestFace
                 label faceI = ppFaces[hitI];
                 faceDisp[faceI] = hitInfo[hitI].hitPoint() - fc[hitI];
                 faceSurfaceNormal[faceI] = hitNormal[hitI];
-                faceSurfaceRegion[faceI] = surfaces.globalRegion
+                faceSurfaceGlobalRegion[faceI] = surfaces.globalRegion
                 (
                     hitSurface[hitI],
                     hitRegion[hitI]
@@ -477,7 +477,11 @@ void Foam::autoSnapDriver::calcNearestFace
             label faceI = ppFaces[hitI];
             faceDisp[faceI] = hitInfo[hitI].hitPoint() - fc[hitI];
             faceSurfaceNormal[faceI] = hitNormal[hitI];
-            faceSurfaceRegion[faceI] = hitRegion[hitI];
+            faceSurfaceGlobalRegion[faceI] = surfaces.globalRegion
+            (
+                hitSurface[hitI],
+                hitRegion[hitI]
+            );
         }
     }
 
@@ -517,6 +521,169 @@ void Foam::autoSnapDriver::calcNearestFace
 }
 
 
+// Collect (possibly remote) per point data of all surrounding faces
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// - faceSurfaceNormal
+// - faceDisp
+// - faceCentres&faceNormal
+void Foam::autoSnapDriver::calcNearestFacePointProperties
+(
+    const label iter,
+    const indirectPrimitivePatch& pp,
+
+    const vectorField& faceDisp,
+    const vectorField& faceSurfaceNormal,
+    const labelList& faceSurfaceGlobalRegion,
+
+    List<List<point> >& pointFaceSurfNormals,
+    List<List<point> >& pointFaceDisp,
+    List<List<point> >& pointFaceCentres,
+    List<labelList>&    pointFacePatchID
+) const
+{
+    const fvMesh& mesh = meshRefiner_.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)
+
+    pointFaceSurfNormals.setSize(pp.nPoints());
+    pointFaceDisp.setSize(pp.nPoints());
+    pointFaceCentres.setSize(pp.nPoints());
+    pointFacePatchID.setSize(pp.nPoints());
+
+    // Fill local data
+    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());
+
+        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] = globalToPatch_[faceSurfaceGlobalRegion[faceI]];
+        }
+    }
+
+
+    // Collect additionally 'normal' boundary faces for boundaryPoints of pp
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // points on the boundary of pp should pick up non-pp normals
+    // as well for the feature-reconstruction to behave correctly.
+    // (the movement is already constrained outside correctly so it
+    //  is only that the unconstrained attraction vector is calculated
+    //  correctly)
+    {
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+        labelList patchID(pbm.patchID());
+
+        // Unmark all non-coupled boundary faces
+        forAll(pbm, patchI)
+        {
+            const polyPatch& pp = pbm[patchI];
+
+            if (pp.coupled() || isA<emptyPolyPatch>(pp))
+            {
+                forAll(pp, i)
+                {
+                    label meshFaceI = pp.start()+i;
+                    patchID[meshFaceI-mesh.nInternalFaces()] = -1;
+                }
+            }
+        }
+
+        // Remove any meshed faces
+        forAll(pp.addressing(), i)
+        {
+            label meshFaceI = pp.addressing()[i];
+            patchID[meshFaceI-mesh.nInternalFaces()] = -1;
+        }
+
+        // See if pp point uses any non-meshed boundary faces
+
+        const labelList& boundaryPoints = pp.boundaryPoints();
+        forAll(boundaryPoints, i)
+        {
+            label pointI = boundaryPoints[i];
+            label meshPointI = pp.meshPoints()[pointI];
+            const point& pt = mesh.points()[meshPointI];
+            const labelList& pFaces = mesh.pointFaces()[meshPointI];
+
+            List<point>& pNormals = pointFaceSurfNormals[pointI];
+            List<point>& pDisp = pointFaceDisp[pointI];
+            List<point>& pFc = pointFaceCentres[pointI];
+            labelList& pFid = pointFacePatchID[pointI];
+
+            forAll(pFaces, i)
+            {
+                label meshFaceI = pFaces[i];
+                if (!mesh.isInternalFace(meshFaceI))
+                {
+                    label patchI = patchID[meshFaceI-mesh.nInternalFaces()];
+
+                    if (patchI != -1)
+                    {
+                        vector fn = mesh.faceAreas()[meshFaceI];
+                        pNormals.append(fn/mag(fn));
+                        pDisp.append(mesh.faceCentres()[meshFaceI]-pt);
+                        pFc.append(mesh.faceCentres()[meshFaceI]);
+                        pFid.append(patchI);
+                    }
+                }
+            }
+        }
+    }
+
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        pointFaceSurfNormals,
+        listPlusEqOp<point>(),
+        List<point>(),
+        listTransform()
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        pointFaceDisp,
+        listPlusEqOp<point>(),
+        List<point>(),
+        listTransform()
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        pointFaceCentres,
+        listPlusEqOp<point>(),
+        List<point>(),
+        listTransform()
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        pointFacePatchID,
+        listPlusEqOp<label>(),
+        List<label>()
+    );
+}
+
+
 // Gets passed in offset to nearest point on feature edge. Calculates
 // if the point has a different number of faces on either side of the feature
 // and if so attracts the point to that non-dominant plane.
@@ -580,56 +747,7 @@ Foam::pointIndexHit Foam::autoSnapDriver::findMultiPatchPoint
     }
     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
 (
@@ -1361,6 +1479,7 @@ void Foam::autoSnapDriver::determineFeatures
 (
     const label iter,
     const scalar featureCos,
+    const bool multiRegionFeatureSnap,
 
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
@@ -1464,69 +1583,72 @@ void Foam::autoSnapDriver::determineFeatures
             if (patchConstraints[pointI].first() == 1)
             {
                 // Flat surface. Check for different patchIDs
-                pointIndexHit multiPatchPt
-                (
-                    findMultiPatchPoint
-                    (
-                        pt,
-                        pointFacePatchID[pointI],
-                        pointFaceCentres[pointI]
-                    )
-                );
-                if (multiPatchPt.hit())
+                if (multiRegionFeatureSnap)
                 {
-                    // 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
+                    pointIndexHit multiPatchPt
                     (
-                        pp,
-                        snapDist,
-                        pointI,
-                        multiPatchPt.hitPoint(),        //estimatedPt
+                        findMultiPatchPoint
+                        (
+                            pt,
+                            pointFacePatchID[pointI],
+                            pointFaceCentres[pointI]
+                        )
+                    );
+                    if (multiPatchPt.hit())
+                    {
+                        // 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,
+                            featI,
+                            edgeAttractors,
+                            edgeConstraints,
 
-                        patchAttraction,
-                        patchConstraints
-                    );
+                            patchAttraction,
+                            patchConstraints
+                        );
 
-                    if (nearInfo.hit())
-                    {
-                        // Dump
-                        if (featureEdgeStr.valid())
+                        if (nearInfo.hit())
                         {
-                            meshTools::writeOBJ(featureEdgeStr(), pt);
-                            featureEdgeVertI++;
-                            meshTools::writeOBJ
-                            (
-                                featureEdgeStr(),
-                                nearInfo.hitPoint()
-                            );
-                            featureEdgeVertI++;
-                            featureEdgeStr()
-                                << "l " << featureEdgeVertI-1 << ' '
-                                << featureEdgeVertI << nl;
+                            // 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())
+                        else
                         {
-                            meshTools::writeOBJ(missedEdgeStr(), pt);
-                            missedVertI++;
-                            meshTools::writeOBJ
-                            (
-                                missedEdgeStr(),
-                                nearInfo.missPoint()
-                            );
-                            missedVertI++;
-                            missedEdgeStr()
-                                << "l " << missedVertI-1 << ' '
-                                << missedVertI << nl;
+                            if (missedEdgeStr.valid())
+                            {
+                                meshTools::writeOBJ(missedEdgeStr(), pt);
+                                missedVertI++;
+                                meshTools::writeOBJ
+                                (
+                                    missedEdgeStr(),
+                                    nearInfo.missPoint()
+                                );
+                                missedVertI++;
+                                missedEdgeStr()
+                                    << "l " << missedVertI-1 << ' '
+                                    << missedVertI << nl;
+                            }
                         }
                     }
                 }
@@ -1645,6 +1767,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
 (
     const label iter,
     const scalar featureCos,
+    const bool multiRegionFeatureSnap,
 
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
@@ -1697,6 +1820,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     (
         iter,
         featureCos,
+        multiRegionFeatureSnap,
 
         pp,
         snapDist,
@@ -2131,6 +2255,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     //MEJ: any faces that have multi-patch points only keep the multi-patch
     //     points. The other points on the face will be dragged along
     //     (hopefully)
+    if (multiRegionFeatureSnap)
     {
         autoPtr<OFstream> multiPatchStr;
         if (debug&meshRefinement::OBJINTERSECTIONS)
@@ -2536,10 +2661,12 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 {
     const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
     const Switch explicitFeatureAttraction = snapParams.explicitFeatureSnap();
+    const Switch multiRegionFeatureSnap = snapParams.multiRegionFeatureSnap();
 
     Info<< "Overriding displacement on features :" << nl
-        << "   implicit features : " << implicitFeatureAttraction << nl
-        << "   explicit features : " << explicitFeatureAttraction << nl
+        << "   implicit features    : " << implicitFeatureAttraction << nl
+        << "   explicit features    : " << explicitFeatureAttraction << nl
+        << "   multi-patch features : " << multiRegionFeatureSnap << nl
         << endl;
 
 
@@ -2547,6 +2674,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     const pointField& localPoints = pp.localPoints();
     const fvMesh& mesh = meshRefiner_.mesh();
 
+
     // Displacement and orientation per pp face
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -2554,7 +2682,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);
+    labelList faceSurfaceGlobalRegion(pp.size(), -1);
     vectorField faceRotation(pp.size(), vector::zero);
 
     calcNearestFace
@@ -2563,7 +2691,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         pp,
         faceDisp,
         faceSurfaceNormal,
-        faceSurfaceRegion,
+        faceSurfaceGlobalRegion,
         faceRotation
     );
 
@@ -2587,158 +2715,25 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // - faceSurfaceNormal
     // - faceDisp
-    // - faceRotation
     // - 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)
-
     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());
 
-    // Fill local data
-    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>& pRot = pointFaceRotation[pointI];
-        //pRot.setSize(pFaces.size());
-        List<point>& pFc = pointFaceCentres[pointI];
-        pFc.setSize(pFaces.size());
-        labelList& pFid = pointFacePatchID[pointI];
-        pFid.setSize(pFaces.size());
-
-        forAll(pFaces, i)
-        {
-            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);
-        }
-    }
-
-
-    // Collect additionally 'normal' boundary faces for boundaryPoints of pp
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // points on the boundary of pp should pick up non-pp normals
-    // as well for the feature-reconstruction to behave correctly.
-    // (the movement is already constrained outside correctly so it
-    //  is only that the unconstrained attraction vector is calculated
-    //  correctly)
-    {
-        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
-        labelList patchID(pbm.patchID());
-
-        // Unmark all non-coupled boundary faces
-        forAll(pbm, patchI)
-        {
-            const polyPatch& pp = pbm[patchI];
-
-            if (pp.coupled() || isA<emptyPolyPatch>(pp))
-            {
-                forAll(pp, i)
-                {
-                    label meshFaceI = pp.start()+i;
-                    patchID[meshFaceI-mesh.nInternalFaces()] = -1;
-                }
-            }
-        }
-
-        // Remove any meshed faces
-        forAll(pp.addressing(), i)
-        {
-            label meshFaceI = pp.addressing()[i];
-            patchID[meshFaceI-mesh.nInternalFaces()] = -1;
-        }
-
-        // See if pp point uses any non-meshed boundary faces
-
-        const labelList& boundaryPoints = pp.boundaryPoints();
-        forAll(boundaryPoints, i)
-        {
-            label pointI = boundaryPoints[i];
-            label meshPointI = pp.meshPoints()[pointI];
-            const point& pt = mesh.points()[meshPointI];
-            const labelList& pFaces = mesh.pointFaces()[meshPointI];
-
-            List<point>& pNormals = pointFaceSurfNormals[pointI];
-            List<point>& pDisp = pointFaceDisp[pointI];
-            List<point>& pFc = pointFaceCentres[pointI];
-            labelList& pFid = pointFacePatchID[pointI];
-
-            forAll(pFaces, i)
-            {
-                label meshFaceI = pFaces[i];
-                if (!mesh.isInternalFace(meshFaceI))
-                {
-                    label patchI = patchID[meshFaceI-mesh.nInternalFaces()];
+    calcNearestFacePointProperties
+    (
+        iter,
+        pp,
 
-                    if (patchI != -1)
-                    {
-                        vector fn = mesh.faceAreas()[meshFaceI];
-                        pNormals.append(fn/mag(fn));
-                        pDisp.append(mesh.faceCentres()[meshFaceI]-pt);
-                        pFc.append(mesh.faceCentres()[meshFaceI]);
-                        pFid.append(patchI);
-                    }
-                }
-            }
-        }
-    }
+        faceDisp,
+        faceSurfaceNormal,
+        faceSurfaceGlobalRegion,
 
-    syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
         pointFaceSurfNormals,
-        listPlusEqOp<point>(),
-        List<point>(),
-        listTransform()
-    );
-    syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
         pointFaceDisp,
-        listPlusEqOp<point>(),
-        List<point>(),
-        listTransform()
-    );
-    //syncTools::syncPointList
-    //(
-    //    mesh,
-    //    pp.meshPoints(),
-    //    pointFaceRotation,
-    //    listPlusEqOp<point>(),
-    //    List<point>(),
-    //    listTransform()
-    //);
-    syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
         pointFaceCentres,
-        listPlusEqOp<point>(),
-        List<point>(),
-        listTransform()
-    );
-    syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
-        pointFacePatchID,
-        listPlusEqOp<label>(),
-        List<label>()
+        pointFacePatchID
     );
 
 
@@ -2793,6 +2788,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         (
             iter,
             featureCos,
+            multiRegionFeatureSnap,
 
             pp,
             snapDist,
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C
index b9d6f3b0c6bcb2771130e4c3c3fcd962c1fde4a0..77bff601870e835f49db1f33c27017047c517484 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C
@@ -36,7 +36,11 @@ Foam::snapParameters::snapParameters(const dictionary& dict)
     nSnap_(readLabel(dict.lookup("nRelaxIter"))),
     nFeatureSnap_(dict.lookupOrDefault("nFeatureSnapIter", -1)),
     explicitFeatureSnap_(dict.lookupOrDefault("explicitFeatureSnap", true)),
-    implicitFeatureSnap_(dict.lookupOrDefault("implicitFeatureSnap", false))
+    implicitFeatureSnap_(dict.lookupOrDefault("implicitFeatureSnap", false)),
+    multiRegionFeatureSnap_
+    (
+        dict.lookupOrDefault("multiRegionFeatureSnap", false)
+    )
 {}
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H
index 8103cd4082dac4c468603a86eff1d4c5d9a320f6..294d0de55d5538df2e9557262be482eb1d9b1157 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H
@@ -68,6 +68,8 @@ class snapParameters
 
         const Switch implicitFeatureSnap_;
 
+        const Switch multiRegionFeatureSnap_;
+
 
     // Private Member Functions
 
@@ -134,6 +136,11 @@ public:
                 return implicitFeatureSnap_;
             }
 
+            Switch multiRegionFeatureSnap() const
+            {
+                return multiRegionFeatureSnap_;
+            }
+
 };