diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index 3787156a53..6c1925f064 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -722,10 +722,15 @@ void Foam::autoSnapDriver::preSmoothPatch
     if (debug)
     {
         const_cast<Time&>(mesh.time())++;
-        Pout<< "Writing patch smoothed mesh to time " << meshRefiner_.timeName()
-            << endl;
-
-        mesh.write();
+        Pout<< "Writing patch smoothed mesh to time "
+            << meshRefiner_.timeName() << '.' << endl;
+        meshRefiner_.write
+        (
+            debug,
+            mesh.time().path()/meshRefiner_.timeName()
+        );
+        Pout<< "Dumped mesh in = "
+            << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
     }
 
     Info<< "Patch points smoothed in = "
@@ -956,25 +961,6 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
         vector(GREAT, GREAT, GREAT)     // null value (note: cannot use VGREAT)
     );
 
-
-    // Check for displacement being outwards.
-    outwardsDisplacement(pp, patchDisp);
-
-    // Set initial distribution of displacement field (on patches) from
-    // patchDisp and make displacement consistent with b.c. on displacement
-    // pointVectorField.
-    meshMover.setDisplacement(patchDisp);
-
-    if (debug)
-    {
-        dumpMove
-        (
-            mesh.time().path()/"patchDisplacement.obj",
-            pp.localPoints(),
-            pp.localPoints() + patchDisp
-        );
-    }
-
     return patchDisp;
 }
 
@@ -1128,8 +1114,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
     indirectPrimitivePatch& pp = ppPtr();
 
     // Divide surfaces into zoned and unzoned
-    labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
-    labelList unzonedSurfaces = meshRefiner_.surfaces().getUnnamedSurfaces();
+    labelList zonedSurfaces = surfaces.getNamedSurfaces();
+    labelList unzonedSurfaces = surfaces.getUnnamedSurfaces();
 
 
     // Faces that do not move
@@ -1344,19 +1330,6 @@ void Foam::autoSnapDriver::doSnap
         // Pre-smooth patch vertices (so before determining nearest)
         preSmoothPatch(snapParams, nInitErrors, baffles, meshMover);
 
-        if (debug)
-        {
-            Pout<< "Writing patch smoothed mesh to time "
-                << meshRefiner_.timeName() << '.' << endl;
-            meshRefiner_.write
-            (
-                debug,
-                mesh.time().path()/meshRefiner_.timeName()
-            );
-            Pout<< "Dumped mesh in = "
-                << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
-        }
-
 
         for (label iter = 0; iter < nFeatIter; iter++)
         {
@@ -1382,6 +1355,15 @@ void Foam::autoSnapDriver::doSnap
                 );
             }
 
+            // Check for displacement being outwards.
+            outwardsDisplacement(ppPtr(), disp);
+
+            // Set initial distribution of displacement field (on patches)
+            // from patchDisp and make displacement consistent with b.c.
+            // on displacement pointVectorField.
+            meshMover.setDisplacement(disp);
+
+
             if (debug&meshRefinement::OBJINTERSECTIONS)
             {
                 dumpMove
@@ -1467,7 +1449,7 @@ void Foam::autoSnapDriver::doSnap
         meshRefiner_.write
         (
             debug,
-            mesh.time().path()/meshRefiner_.timeName()
+            meshRefiner_.timeName()
         );
     }
 }
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index 4b0e1a30fd..59cd772a98 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -117,12 +117,12 @@ class autoSnapDriver
                     const List<pointConstraint>& constraints,
                     vectorField& disp
                 ) const;
-                void calcNearest
-                (
-                    const pointField& points,
-                    vectorField& disp,
-                    vectorField& surfaceNormal
-                ) const;
+                //void calcNearest
+                //(
+                //    const pointField& points,
+                //    vectorField& disp,
+                //    vectorField& surfaceNormal
+                //) const;
                 void calcNearestFace
                 (
                     const label iter,
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index cc73efbe51..0227234822 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -135,63 +135,6 @@ void Foam::autoSnapDriver::smoothAndConstrain
 }
 
 
-void Foam::autoSnapDriver::calcNearest
-(
-    const pointField& points,
-    vectorField& disp,
-    vectorField& surfaceNormal
-) const
-{
-    const refinementSurfaces& surfaces = meshRefiner_.surfaces();
-
-    // Displacement and orientation per pp face.
-    disp.setSize(points.size());
-    disp = vector::zero;
-    surfaceNormal.setSize(points.size());
-    surfaceNormal = vector::zero;
-
-    {
-        // Divide surfaces into zoned and unzoned
-        labelList zonedSurfaces =
-            meshRefiner_.surfaces().getNamedSurfaces();
-        labelList unzonedSurfaces =
-            meshRefiner_.surfaces().getUnnamedSurfaces();
-
-        {
-            List<pointIndexHit> hitInfo;
-            labelList hitSurface;
-            labelList hitRegion;
-            surfaces.findNearestRegion
-            (
-                unzonedSurfaces,
-                points,
-                sqr(scalarField(points.size(), GREAT)),// sqr of attract dist
-                hitSurface,
-                hitInfo,
-                hitRegion,
-                surfaceNormal
-            );
-
-            forAll(hitInfo, i)
-            {
-                if (hitInfo[i].hit())
-                {
-                    disp[i] =
-                        hitInfo[i].hitPoint()
-                      - points[i];
-                }
-                else
-                {
-                    WarningIn("calcNearest(..)")
-                        << "Did not hit anything from face:" << i
-                        << " at:" << points[i] << endl;
-                }
-            }
-        }
-    }
-}
-
-
 void Foam::autoSnapDriver::calcNearestFace
 (
     const label iter,
@@ -202,6 +145,7 @@ void Foam::autoSnapDriver::calcNearestFace
 ) const
 {
     const fvMesh& mesh = meshRefiner_.mesh();
+    const refinementSurfaces& surfaces = meshRefiner_.surfaces();
 
     // Displacement and orientation per pp face.
     faceDisp.setSize(pp.size());
@@ -209,7 +153,148 @@ void Foam::autoSnapDriver::calcNearestFace
     faceSurfaceNormal.setSize(pp.size());
     faceSurfaceNormal = vector::zero;
 
-    calcNearest(pp.faceCentres(), faceDisp, faceSurfaceNormal);
+    // Divide surfaces into zoned and unzoned
+    labelList zonedSurfaces = surfaces.getNamedSurfaces();
+    labelList unzonedSurfaces = surfaces.getUnnamedSurfaces();
+
+    // Per pp face the current surface snapped to
+    labelList snapSurf(pp.size(), -1);
+
+
+    // Do zoned surfaces
+    // ~~~~~~~~~~~~~~~~~
+    // Zoned faces only attract to corresponding surface
+
+    // Extract faces per zone
+    const wordList& faceZoneNames = surfaces.faceZoneNames();
+
+    forAll(zonedSurfaces, i)
+    {
+        label zoneSurfI = zonedSurfaces[i];
+
+        // Get indices of faces on pp that are also in zone
+        label zoneI = mesh.faceZones().findZoneID(faceZoneNames[zoneSurfI]);
+        if (zoneI == -1)
+        {
+            FatalErrorIn
+            (
+                "autoSnapDriver::calcNearestFace(..)"
+            )   << "Problem. Cannot find zone " << faceZoneNames[zoneSurfI]
+                << exit(FatalError);
+        }
+        const faceZone& fZone = mesh.faceZones()[zoneI];
+        PackedBoolList isZonedFace(mesh.nFaces());
+        forAll(fZone, i)
+        {
+            isZonedFace[fZone[i]] = 1;
+        }
+
+        DynamicList<label> ppFaces(fZone.size());
+        DynamicList<label> meshFaces(fZone.size());
+        forAll(pp.addressing(), i)
+        {
+            if (isZonedFace[pp.addressing()[i]])
+            {
+                snapSurf[i] = zoneSurfI;
+                ppFaces.append(i);
+                meshFaces.append(pp.addressing()[i]);
+            }
+        }
+
+        //Pout<< "For faceZone " << fZone.name()
+        //    << " found " << ppFaces.size() << " out of " << pp.size()
+        //    << endl;
+
+        pointField fc
+        (
+            indirectPrimitivePatch
+            (
+                IndirectList<face>(mesh.faces(), meshFaces),
+                mesh.points()
+            ).faceCentres()
+        );
+
+        List<pointIndexHit> hitInfo;
+        labelList hitSurface;
+        labelList hitRegion;
+        vectorField hitNormal;
+        surfaces.findNearestRegion
+        (
+            labelList(1, zoneSurfI),
+            fc,
+            sqr(scalarField(fc.size(), GREAT)),// sqr of attract dist
+            hitSurface,
+            hitInfo,
+            hitRegion,
+            hitNormal
+        );
+
+        forAll(hitInfo, hitI)
+        {
+            if (hitInfo[hitI].hit())
+            {
+                label faceI = ppFaces[hitI];
+                faceDisp[faceI] = hitInfo[hitI].hitPoint() - fc[hitI];
+                faceSurfaceNormal[faceI] = hitNormal[hitI];
+            }
+        }
+    }
+
+
+    // Do unzoned surfaces
+    // ~~~~~~~~~~~~~~~~~~~
+    // Unzoned faces attract to any unzoned surface
+
+    DynamicList<label> ppFaces(pp.size());
+    DynamicList<label> meshFaces(pp.size());
+    forAll(pp.addressing(), i)
+    {
+        if (snapSurf[i] == -1)
+        {
+            ppFaces.append(i);
+            meshFaces.append(pp.addressing()[i]);
+        }
+    }
+    //Pout<< "Found " << ppFaces.size() << " unzoned faces out of " << pp.size()
+    //    << endl;
+
+    pointField fc
+    (
+        indirectPrimitivePatch
+        (
+            IndirectList<face>(mesh.faces(), meshFaces),
+            mesh.points()
+        ).faceCentres()
+    );
+
+    List<pointIndexHit> hitInfo;
+    labelList hitSurface;
+    labelList hitRegion;
+    vectorField hitNormal;
+    surfaces.findNearestRegion
+    (
+        unzonedSurfaces,
+        fc,
+        sqr(scalarField(fc.size(), GREAT)),// sqr of attract dist
+        hitSurface,
+        hitInfo,
+        hitRegion,
+        hitNormal
+    );
+
+    forAll(hitInfo, hitI)
+    {
+        if (hitInfo[hitI].hit())
+        {
+            label faceI = ppFaces[hitI];
+            faceDisp[faceI] = hitInfo[hitI].hitPoint() - fc[hitI];
+            faceSurfaceNormal[faceI] = hitNormal[hitI];
+        }
+    }
+
+
+    // Determine rotation
+    // ~~~~~~~~~~~~~~~~~~
 
     // Determine rotation axis
     faceRotation.setSize(pp.size());
@@ -243,104 +328,6 @@ void Foam::autoSnapDriver::calcNearestFace
 }
 
 
-void Foam::autoSnapDriver::interpolateFaceToPoint
-(
-    const label iter,
-    const indirectPrimitivePatch& pp,
-    const vectorField& faceSurfaceNormal,
-    const vectorField& faceDisp,
-    const vectorField& faceRotation,
-
-    vectorField& patchDisp,
-    vectorField& patchRotationDisp
-) const
-{
-    const fvMesh& mesh = meshRefiner_.mesh();
-
-    // Displacement due to normal snapping
-    patchDisp.setSize(pp.nPoints());
-    patchDisp = vector::zero;
-    // Displacement due to rotation
-    patchRotationDisp.setSize(pp.nPoints());
-    patchRotationDisp = vector::zero;
-
-
-    // Unweighted interpolation
-    // ~~~~~~~~~~~~~~~~~~~~~~~~
-
-    labelList nPatchDisp(pp.nPoints(), 0.0);
-
-    forAll(pp, faceI)
-    {
-        const face& f = pp.localFaces()[faceI];
-        const point& fc = pp.faceCentres()[faceI];
-
-        forAll(f, fp)
-        {
-            label pointI = f[fp];
-            const point& pt = pp.localPoints()[pointI];
-
-            patchDisp[pointI] += faceDisp[faceI];
-            patchRotationDisp[pointI] += (faceRotation[faceI]^(pt-fc));
-            nPatchDisp[pointI]++;
-        }
-    }
-
-    syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
-        patchDisp,
-        plusEqOp<vector>(),
-        vector::zero,
-        mapDistribute::transform()
-    );
-    syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
-        patchRotationDisp,
-        plusEqOp<vector>(),
-        vector::zero,
-        mapDistribute::transform()
-    );
-    syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
-        nPatchDisp,
-        plusEqOp<label>(),
-        0,
-        mapDistribute::transform()
-    );
-
-    forAll(patchDisp, pointI)
-    {
-        patchDisp[pointI] /= nPatchDisp[pointI];
-        patchRotationDisp[pointI] /= nPatchDisp[pointI];
-    }
-
-
-    if (debug&meshRefinement::OBJINTERSECTIONS)
-    {
-        dumpMove
-        (
-            mesh.time().path()
-          / "patchDisp_" + name(iter) + ".obj",
-            pp.localPoints(),
-            pp.localPoints() + patchDisp
-        );
-        dumpMove
-        (
-            mesh.time().path()
-          / "patchRotationDisp_" + name(iter) + ".obj",
-            pp.localPoints(),
-            pp.localPoints() + patchRotationDisp
-        );
-    }
-}
-
-
 // Gets passed in offset to nearest point on feature edge. Calculates
 // if the point has a different number of faces on either side of the feature
 // and if so attracts the point to that non-dominant plane.
@@ -500,8 +487,6 @@ void Foam::autoSnapDriver::binFeatureFaces
 
     const label pointI,
 
-//    const vectorField& faceSurfaceNormal,
-//    const vectorField& faceDisp,
     const List<List<point> >& pointFaceNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
@@ -547,9 +532,6 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
 
-    //const vectorField& faceSurfaceNormal,
-    //const vectorField& faceDisp,
-    //const vectorField& faceRotation,
     const List<List<point> >& pointFaceNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
@@ -609,8 +591,6 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
             snapDist,
             pointI,
 
-            //faceSurfaceNormal,
-            //faceDisp,
             pointFaceNormals,
             pointFaceDisp,
             pointFaceCentres,
@@ -741,9 +721,6 @@ void Foam::autoSnapDriver::determineAllFeatures
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
 
-    //const vectorField& faceSurfaceNormal,
-    //const vectorField& faceDisp,
-    //const vectorField& faceRotation,
     const List<List<point> >& pointFaceNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
@@ -1778,33 +1755,8 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         faceRotation
     );
 
-    //
-    // Interpolate face-wise data to points
-    //
-
-
-    // Displacement due to snapping of face centres
-    pointField patchDisp(localPoints.size(), vector::zero);
-    // Displacement due to rotation
-    pointField patchRotationDisp(localPoints.size(), vector::zero);
-
-    interpolateFaceToPoint
-    (
-        iter,
-        pp,
-        faceSurfaceNormal,
-        faceDisp,
-        faceRotation,
-
-        patchDisp,
-        patchRotationDisp
-    );
-
-    // Override patchDisp with exact nearest
-    {
-        pointField pointSurfaceNormal;
-        calcNearest(pp.localPoints(), patchDisp, pointSurfaceNormal);
-    }
+    // Start off with nearest point on surface.
+    vectorField patchDisp = nearestDisp;
 
 
     // Collect (possibly remote) face-wise data on coupled points.
@@ -1818,7 +1770,6 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 
     List<List<point> > pointFaceNormals(pp.nPoints());
     List<List<point> > pointFaceDisp(pp.nPoints());
-    List<List<point> > pointFaceRotation(pp.nPoints());
     List<List<point> > pointFaceCentres(pp.nPoints());
 
     // Fill local data
@@ -1829,15 +1780,12 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         pNormals.setSize(pFaces.size());
         List<point>& pDisp = pointFaceDisp[pointI];
         pDisp.setSize(pFaces.size());
-        List<point>& pRot = pointFaceRotation[pointI];
-        pRot.setSize(pFaces.size());
         List<point>& pFc = pointFaceCentres[pointI];
         pFc.setSize(pFaces.size());
         forAll(pFaces, i)
         {
             pNormals[i] = faceSurfaceNormal[pFaces[i]];
             pDisp[i] = faceDisp[pFaces[i]];
-            pRot[i] = faceRotation[pFaces[i]];
             pFc[i] = pp.faceCentres()[pFaces[i]];
         }
     }
@@ -1861,15 +1809,6 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         listTransform()
     );
     syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
-        pointFaceRotation,
-        listPlusEqOp(),
-        List<point>(),
-        listTransform()
-    );
-    syncTools::syncPointList
     (
         mesh,
         pp.meshPoints(),
@@ -1881,60 +1820,6 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 
 
 
-//XXXXXXX
-//    // Indices of pp points that are coupled
-//    labelList coupledPointLabels(pp.nPoints());
-//    // Corresponding coupled point index
-//    labelList coupledPointAddr(pp.nPoints());
-//    {
-//        label coupledI = 0;
-//        forAll(pp.meshPoints(), pointI)
-//        {
-//            label meshPointI = pp.meshPoints()[pointI];
-//            Map<label>::const_iterator fnd = coupledPatchMP.find(meshPointI);
-//            if (fnd != coupledPatchMP.end())
-//            {
-//                coupledPointLabels[coupledI] = pointI;
-//                coupledPointAddr[coupledI] = fnd();
-//                coupledI++;
-//            }
-//        }
-//        coupledPointLabels.setSize(coupledI);
-//        coupledPointAddr.setSize(coupledI);
-//    }
-//
-//    List<List<point> > pointFaceNormals(map.constructSize());
-//    List<List<point> > pointFaceDisp(map.constructSize());
-//    List<List<point> > pointFaceRotation(map.constructSize());
-//    List<List<point> > pointFaceCentres(map.constructSize());
-//
-//    // Collect local data.
-//    forAll(coupledPointLabels, coupledI)
-//    {
-//        label pointI = coupledPointLabels[coupledI];
-//        label coupledPointI = coupledPointAddr[coupledI];
-//        const labelList& pFaces = pp.pointFaces()[pointI];
-//        List<point>& pNormals = pointFaceNormals[coupledPointI];
-//        pNormals.setSize(pFaces.size());
-//        List<point>& pDisp = pointFaceDisp[coupledPointI];
-//        pDisp.setSize(pFaces.size());
-//        List<point>& pRot = pointFaceRotation[coupledPointI];
-//        pRot.setSize(pFaces.size());
-//        List<point>& pFc = pointFaceCentres[coupledPointI];
-//        pFc.setSize(pFaces.size());
-//        forAll(pFaces, i)
-//        {
-//            pNormals[i] = faceSurfaceNormal[pFaces[i]];
-//            pDisp[i] = pointFaceDisp[pFaces[i]];
-//            pRot[i] = pointFaceRotation[pFaces[i]];
-//            pFc[i] = pp.faceCentres()[pFaces[i]];
-//        }
-//    }
-//
-//    // Pull remote data into local slots
-//XXXXXXX
-
-
     // Nearest feature
     vectorField patchAttraction(localPoints.size(), vector::zero);
     // Constraints at feature
@@ -1973,9 +1858,6 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         << "     linear   : max:" << gMax(patchDisp)
         << " avg:" << gAverage(patchDisp)
         << endl
-        << "     rotation : max:" << gMax(patchRotationDisp)
-        << " avg:" << gAverage(patchRotationDisp)
-        << endl
         << "     feature  : max:" << gMax(patchAttraction)
         << " avg:" << gAverage(patchAttraction)
         << endl;
@@ -1985,20 +1867,13 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     // - patchDisp          : point movement to go to nearest point on surface
     //                       (either direct or through interpolation of
     //                        face nearest)
-    // - patchRotationDisp  : point movement to align face normal with surface.
-    //
     // - patchAttraction    : direct attraction to features
     // - patchConstraints   : type of features
 
-    // Use any combination of patchDisp, patchRotationDisp and direct feature
+    // Use any combination of patchDisp and direct feature
     // attraction.
 
 
-    // Bit of patchRotation
-    //patchDisp = 0.3*patchRotationDisp+0.7*patchDisp;
-    // Disable patchRotationDisp
-    patchRotationDisp = vector::zero;
-
     // Mix with direct feature attraction
     forAll(patchConstraints, pointI)
     {
@@ -2017,13 +1892,6 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     //    pp.localPoints(),
     //    pp.localPoints() + patchDisp
     //);
-    //dumpMove
-    //(
-    //    mesh.time().path()
-    //  / "rotationPatchDisp_" + name(iter) + ".obj",
-    //    pp.localPoints(),
-    //    pp.localPoints() + patchRotationDisp
-    //);
 
 
 
@@ -2130,15 +1998,6 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         vector(GREAT, GREAT, GREAT)     // null value (note: cannot use VGREAT)
     );
 
-
-    // Check for displacement being outwards.
-    outwardsDisplacement(pp, patchDisp);
-
-    // Set initial distribution of displacement field (on patches) from
-    // patchDisp and make displacement consistent with b.c. on displacement
-    // pointVectorField.
-    meshMover.setDisplacement(patchDisp);
-
     return patchDisp;
 }
 
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementMerge.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementMerge.C
index 39a6833b55..dc16bd760f 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementMerge.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementMerge.C
@@ -349,10 +349,32 @@ Foam::label Foam::meshRefinement::mergePatchFacesUndo
 
         faceCombiner.updateMesh(map);
 
-        updateMesh(map, labelList(0));
+        // Get the kept faces that need to be recalculated.
+        // Merging two boundary faces might shift the cell centre
+        // (unless the faces are absolutely planar)
+        labelHashSet retestFaces(6*allFaceSets.size());
+
+        forAll(allFaceSets, setI)
+        {
+            label oldMasterI = allFaceSets[setI][0];
+
+            label faceI = map().reverseFaceMap()[oldMasterI];
+
+            // faceI is always uncoupled boundary face
+            const cell& cFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
+
+            forAll(cFaces, i)
+            {
+                retestFaces.insert(cFaces[i]);
+            }
+        }
+        updateMesh(map, retestFaces.toc());
 
         if (debug)
         {
+            // Check sync
+            checkData();
+
             Pout<< "Writing initial merged-faces mesh to time "
                 << timeName() << nl << endl;
             write();
@@ -529,11 +551,30 @@ Foam::label Foam::meshRefinement::mergePatchFacesUndo
             inplaceMapKey(map().reverseFaceMap(), restoredFaces);
             inplaceMapKey(map().reverseCellMap(), restoredCells);
 
+
+            // Get the kept faces that need to be recalculated.
+            // Merging two boundary faces might shift the cell centre
+            // (unless the faces are absolutely planar)
+            labelHashSet retestFaces(6*restoredFaces.size());
+
+            forAll(restoredFaces, setI)
+            {
+                label faceI = restoredFaces[setI];
+                // faceI is always uncoupled boundary face
+                const cell& cFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
+
+                forAll(cFaces, i)
+                {
+                    retestFaces.insert(cFaces[i]);
+                }
+            }
+
+
             // Experimental:restore all points/face/cells in maps
             updateMesh
             (
                 map,
-                labelList(0),       // changedFaces
+                retestFaces.toc(),
                 restoredPoints,
                 restoredFaces,
                 restoredCells
@@ -541,6 +582,9 @@ Foam::label Foam::meshRefinement::mergePatchFacesUndo
 
             if (debug)
             {
+                // Check sync
+                checkData();
+
                 Pout<< "Writing merged-faces mesh to time "
                     << timeName() << nl << endl;
                 write();