From 15fed04026425aaad7ecf37306661f0de8300aef Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 3 Jun 2013 13:33:02 +0100
Subject: [PATCH] ENH: autoHexMesh: added geometric problem cell removal

---
 .../autoHexMeshDriver/autoRefineDriver.C      |  47 ++-
 .../autoHexMeshDriver/autoRefineDriver.H      |   5 +
 .../autoHexMeshDriver/autoSnapDriver.C        | 315 +++++++++++++-
 .../autoHexMeshDriver/autoSnapDriver.H        |  38 +-
 .../refinementParameters.C                    |  14 +-
 .../refinementParameters.H                    |  11 +
 .../meshRefinement/meshRefinement.H           |  19 +
 .../meshRefinement/meshRefinementBaffles.C    | 124 ++----
 .../meshRefinementProblemCells.C              | 397 +++++++++++-------
 .../refinementSurfaces/refinementSurfaces.C   |   7 +-
 .../refinementSurfaces/refinementSurfaces.H   |   3 +-
 11 files changed, 677 insertions(+), 303 deletions(-)

diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index d78866bcd7..82bfa13522 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -35,6 +35,7 @@ License
 #include "shellSurfaces.H"
 #include "mapDistributePolyMesh.H"
 #include "unitConversion.H"
+#include "snapParameters.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -610,9 +611,16 @@ Foam::label Foam::autoRefineDriver::danglingCellRefine
             << " cells (out of " << mesh.globalData().nTotalCells()
             << ')' << endl;
 
-        // Stop when no cells to refine. No checking of minRefineCells since
-        // too few cells
-        if (nCellsToRefine == 0)
+        // Stop when no cells to refine. After a few iterations check if too
+        // few cells
+        if
+        (
+            nCellsToRefine == 0
+         || (
+                iter >= 1
+             && nCellsToRefine <= refineParams.minRefineCells()
+            )
+        )
         {
             Info<< "Stopping refining since too few cells selected."
                 << nl << endl;
@@ -870,6 +878,7 @@ Foam::label Foam::autoRefineDriver::shellRefine
 void Foam::autoRefineDriver::baffleAndSplitMesh
 (
     const refinementParameters& refineParams,
+    const snapParameters& snapParams,
     const bool handleSnapProblems,
     const dictionary& motionDict
 )
@@ -887,10 +896,17 @@ void Foam::autoRefineDriver::baffleAndSplitMesh
     meshRefiner_.baffleAndSplitMesh
     (
         handleSnapProblems,             // detect&remove potential snap problem
+
+        // Snap problem cell detection
+        snapParams,
+        refineParams.useTopologicalSnapDetection(),
         false,                          // perpendicular edge connected cells
         scalarField(0),                 // per region perpendicular angle
+
+        // Free standing baffles
         !handleSnapProblems,            // merge free standing baffles?
         refineParams.planarAngle(),
+
         motionDict,
         const_cast<Time&>(mesh.time()),
         globalToMasterPatch_,
@@ -951,6 +967,7 @@ void Foam::autoRefineDriver::zonify
 void Foam::autoRefineDriver::splitAndMergeBaffles
 (
     const refinementParameters& refineParams,
+    const snapParameters& snapParams,
     const bool handleSnapProblems,
     const dictionary& motionDict
 )
@@ -973,10 +990,17 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
     meshRefiner_.baffleAndSplitMesh
     (
         handleSnapProblems,
+
+        // Snap problem cell detection
+        snapParams,
+        refineParams.useTopologicalSnapDetection(),
         handleSnapProblems,                 // remove perp edge connected cells
         perpAngle,                          // perp angle
+
+        // Free standing baffles
         true,                               // merge free standing baffles?
         refineParams.planarAngle(),         // planar angle
+
         motionDict,
         const_cast<Time&>(mesh.time()),
         globalToMasterPatch_,
@@ -1081,6 +1105,7 @@ void Foam::autoRefineDriver::doRefine
 (
     const dictionary& refineDict,
     const refinementParameters& refineParams,
+    const snapParameters& snapParams,
     const bool prepareForSnapping,
     const dictionary& motionDict
 )
@@ -1146,13 +1171,25 @@ void Foam::autoRefineDriver::doRefine
 
     // Introduce baffles at surface intersections. Remove sections unreachable
     // from keepPoint.
-    baffleAndSplitMesh(refineParams, prepareForSnapping, motionDict);
+    baffleAndSplitMesh
+    (
+        refineParams,
+        snapParams,
+        prepareForSnapping,
+        motionDict
+    );
 
     // Mesh is at its finest. Do optional zoning.
     zonify(refineParams);
 
     // Pull baffles apart
-    splitAndMergeBaffles(refineParams, prepareForSnapping, motionDict);
+    splitAndMergeBaffles
+    (
+        refineParams,
+        snapParams,
+        prepareForSnapping,
+        motionDict
+    );
 
     // Do something about cells with refined faces on the boundary
     if (prepareForSnapping)
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
index c44a14a794..a72f3c0796 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
@@ -43,6 +43,8 @@ namespace Foam
 
 // Forward declaration of classes
 class refinementParameters;
+class snapParameters;
+
 class meshRefinement;
 class decompositionMethod;
 class fvMeshDistribute;
@@ -121,6 +123,7 @@ class autoRefineDriver
         void baffleAndSplitMesh
         (
             const refinementParameters& refineParams,
+            const snapParameters& snapParams,
             const bool handleSnapProblems,
             const dictionary& motionDict
         );
@@ -131,6 +134,7 @@ class autoRefineDriver
         void splitAndMergeBaffles
         (
             const refinementParameters& refineParams,
+            const snapParameters& snapParams,
             const bool handleSnapProblems,
             const dictionary& motionDict
         );
@@ -176,6 +180,7 @@ public:
         (
             const dictionary& refineDict,
             const refinementParameters& refineParams,
+            const snapParameters& snapParams,
             const bool prepareForSnapping,
             const dictionary& motionDict
         );
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index 4c262bd132..3927305177 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -121,7 +121,7 @@ Foam::pointField Foam::autoSnapDriver::smoothPatchDisplacement
 (
     const motionSmoother& meshMover,
     const List<labelPair>& baffles
-) const
+)
 {
     const indirectPrimitivePatch& pp = meshMover.patch();
 
@@ -615,14 +615,14 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::mergeZoneBaffles
 
 Foam::scalarField Foam::autoSnapDriver::calcSnapDistance
 (
+    const fvMesh& mesh,
     const snapParameters& snapParams,
     const indirectPrimitivePatch& pp
-) const
+)
 {
     const edgeList& edges = pp.edges();
     const labelListList& pointEdges = pp.pointEdges();
     const pointField& localPoints = pp.localPoints();
-    const fvMesh& mesh = meshRefiner_.mesh();
 
     scalarField maxEdgeLen(localPoints.size(), -GREAT);
 
@@ -655,13 +655,14 @@ Foam::scalarField Foam::autoSnapDriver::calcSnapDistance
 
 void Foam::autoSnapDriver::preSmoothPatch
 (
+    const meshRefinement& meshRefiner,
     const snapParameters& snapParams,
     const label nInitErrors,
     const List<labelPair>& baffles,
     motionSmoother& meshMover
-) const
+)
 {
-    const fvMesh& mesh = meshRefiner_.mesh();
+    const fvMesh& mesh = meshRefiner.mesh();
 
     labelList checkFaces;
 
@@ -724,11 +725,11 @@ void Foam::autoSnapDriver::preSmoothPatch
     {
         const_cast<Time&>(mesh.time())++;
         Info<< "Writing patch smoothed mesh to time "
-            << meshRefiner_.timeName() << '.' << endl;
-        meshRefiner_.write
+            << meshRefiner.timeName() << '.' << endl;
+        meshRefiner.write
         (
             debug,
-            mesh.time().path()/meshRefiner_.timeName()
+            mesh.time().path()/meshRefiner.timeName()
         );
         Info<< "Dumped mesh in = "
             << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
@@ -742,12 +743,11 @@ void Foam::autoSnapDriver::preSmoothPatch
 // Get (pp-local) indices of points that are both on zone and on patched surface
 Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
 (
+    const fvMesh& mesh,
     const indirectPrimitivePatch& pp,
     const word& zoneName
-) const
+)
 {
-    const fvMesh& mesh = meshRefiner_.mesh();
-
     label zoneI = mesh.faceZones().findZoneID(zoneName);
 
     if (zoneI == -1)
@@ -755,7 +755,7 @@ Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
         FatalErrorIn
         (
             "autoSnapDriver::getZoneSurfacePoints"
-            "(const indirectPrimitivePatch&, const word&)"
+            "(const fvMesh&, const indirectPrimitivePatch&, const word&)"
         )   << "Cannot find zone " << zoneName
             << exit(FatalError);
     }
@@ -793,17 +793,17 @@ Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
 
 Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
 (
+    const meshRefinement& meshRefiner,
     const scalarField& snapDist,
-    motionSmoother& meshMover
-) const
+    const indirectPrimitivePatch& pp
+)
 {
     Info<< "Calculating patchDisplacement as distance to nearest surface"
         << " point ..." << endl;
 
-    const indirectPrimitivePatch& pp = meshMover.patch();
     const pointField& localPoints = pp.localPoints();
-    const refinementSurfaces& surfaces = meshRefiner_.surfaces();
-    const fvMesh& mesh = meshRefiner_.mesh();
+    const refinementSurfaces& surfaces = meshRefiner.surfaces();
+    const fvMesh& mesh = meshRefiner.mesh();
 
     // Displacement per patch point
     vectorField patchDisp(localPoints.size(), vector::zero);
@@ -815,9 +815,9 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
 
         // Divide surfaces into zoned and unzoned
         labelList zonedSurfaces =
-            meshRefiner_.surfaces().getNamedSurfaces();
+            meshRefiner.surfaces().getNamedSurfaces();
         labelList unzonedSurfaces =
-            meshRefiner_.surfaces().getUnnamedSurfaces();
+            meshRefiner.surfaces().getUnnamedSurfaces();
 
 
         // 1. All points to non-interface surfaces
@@ -870,6 +870,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
             (
                 getZoneSurfacePoints
                 (
+                    mesh,
                     pp,
                     faceZoneNames[zoneSurfI]
                 )
@@ -966,6 +967,254 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
 }
 
 
+////XXXXXXXXX
+//// Get (pp-local) indices of points that are on both patches
+//Foam::labelList Foam::autoSnapDriver::getPatchSurfacePoints
+//(
+//    const fvMesh& mesh,
+//    const indirectPrimitivePatch& allPp,
+//    const polyPatch& pp
+//)
+//{
+//    // Could use PrimitivePatch & localFaces to extract points but might just
+//    // as well do it ourselves.
+//
+//    boolList pointOnZone(allPp.nPoints(), false);
+//
+//    forAll(pp, i)
+//    {
+//        const face& f = pp[i];
+//
+//        forAll(f, fp)
+//        {
+//            label meshPointI = f[fp];
+//
+//            Map<label>::const_iterator iter =
+//                allPp.meshPointMap().find(meshPointI);
+//
+//            if (iter != allPp.meshPointMap().end())
+//            {
+//                label pointI = iter();
+//                pointOnZone[pointI] = true;
+//            }
+//        }
+//    }
+//
+//    return findIndices(pointOnZone, true);
+//}
+//Foam::vectorField Foam::autoSnapDriver::calcNearestLocalSurface
+//(
+//    const meshRefinement& meshRefiner,
+//    const scalarField& snapDist,
+//    const indirectPrimitivePatch& pp
+//)
+//{
+//    Info<< "Calculating patchDisplacement as distance to nearest"
+//        << " local surface point ..." << endl;
+//
+//    const pointField& localPoints = pp.localPoints();
+//    const refinementSurfaces& surfaces = meshRefiner.surfaces();
+//    const fvMesh& mesh = meshRefiner.mesh();
+//    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+//
+//
+////    // Assume that all patch-internal points get attracted to their surface
+////    // only. So we want to know if point is on multiple regions
+////
+////    labelList minPatch(mesh.nPoints(), labelMax);
+////    labelList maxPatch(mesh.nPoints(), labelMin);
+////
+////    forAll(meshMover.adaptPatchIDs(), i)
+////    {
+////        label patchI = meshMover.adaptPatchIDs()[i];
+////        const labelList& meshPoints = pbm[patchI].meshPoints();
+////
+////        forAll(meshPoints, meshPointI)
+////        {
+////            label meshPointI = meshPoints[meshPointI];
+////            minPatch[meshPointI] = min(minPatch[meshPointI], patchI);
+////            maxPatch[meshPointI] = max(maxPatch[meshPointI], patchI);
+////        }
+////    }
+////
+////    syncTools::syncPointList
+////    (
+////        mesh,
+////        minPatch,
+////        minEqOp<label>(),   // combine op
+////        labelMax            // null value
+////    );
+////    syncTools::syncPointList
+////    (
+////        mesh,
+////        maxPatch,
+////        maxEqOp<label>(),   // combine op
+////        labelMin            // null value
+////    );
+//
+//    // Now all points with minPatch != maxPatch will be on the outside of
+//    // the patch.
+//
+//    // Displacement per patch point
+//    vectorField patchDisp(localPoints.size(), vector::zero);
+//    // Current best snap distance
+//    scalarField minSnapDist(snapDist);
+//    // Current surface snapped to
+//    labelList snapSurf(localPoints.size(), -1);
+//
+//    const labelList& surfaceGeometry = surfaces.surfaces();
+//    forAll(surfaceGeometry, surfI)
+//    {
+//        label geomI = surfaceGeometry[surfI];
+//        const wordList& regNames = allGeometry.regionNames()[geomI];
+//        forAll(regNames, regionI)
+//        {
+//            label globalRegionI = surfaces.globalRegion(surfI, regionI);
+//            // Collect master patch points
+//            label masterPatchI = globalToMasterPatch_[globalRegionI];
+//            label slavePatchI = globalToSlavePatch_[globalRegionI];
+//
+//            labelList patchPointIndices
+//            (
+//                getPatchSurfacePoints
+//                (
+//                    mesh,
+//                    pp,
+//                    pbm[masterPatchI]
+//                )
+//            );
+//
+//            // Find nearest for points both on faceZone and pp.
+//            List<pointIndexHit> hitInfo;
+//            surfaces.findNearest
+//            (
+//                surfI,
+//                regionI,
+//                pointField(localPoints, patchPointIndices),
+//                sqr(scalarField(minSnapDist, patchPointIndices)),
+//                hitInfo
+//            );
+//
+//            forAll(hitInfo, i)
+//            {
+//                label pointI = patchPointIndices[i];
+//
+//                if (hitInfo[i].hit())
+//                {
+//                    const point& pt = hitInfo[i].hitPoint();
+//                    patchDisp[pointI] = pt-localPoints[pointI];
+//                    minSnapDist[pointI] = min
+//                    (
+//                        minSnapDist[pointI],
+//                        mag(patchDisp[pointI])
+//                    );
+//                    snapSurf[pointI] = surfI;
+//                }
+//            }
+//
+//            // Slave patch
+//            if (slavePatchI != masterPatchI)
+//            {
+//                labelList patchPointIndices
+//                (
+//                    getPatchSurfacePoints
+//                    (
+//                        mesh,
+//                        pp,
+//                        pbm[slavePatchI]
+//                    )
+//                );
+//
+//                // Find nearest for points both on faceZone and pp.
+//                List<pointIndexHit> hitInfo;
+//                surfaces.findNearest
+//                (
+//                    surfI,
+//                    regionI,
+//                    pointField(localPoints, patchPointIndices),
+//                    sqr(scalarField(minSnapDist, patchPointIndices)),
+//                    hitInfo
+//                );
+//
+//                forAll(hitInfo, i)
+//                {
+//                    label pointI = patchPointIndices[i];
+//
+//                    if (hitInfo[i].hit())
+//                    {
+//                        const point& pt = hitInfo[i].hitPoint();
+//                        patchDisp[pointI] = pt-localPoints[pointI];
+//                        minSnapDist[pointI] = min
+//                        (
+//                            minSnapDist[pointI],
+//                            mag(patchDisp[pointI])
+//                        );
+//                        snapSurf[pointI] = surfI;
+//                    }
+//                }
+//            }
+//        }
+//    }
+//
+//
+//    // Check if all points are being snapped
+//    forAll(snapSurf, pointI)
+//    {
+//        if (snapSurf[pointI] == -1)
+//        {
+//            WarningIn("autoSnapDriver::calcNearestLocalSurface(..)")
+//                << "For point:" << pointI
+//                << " coordinate:" << localPoints[pointI]
+//                << " did not find any surface within:"
+//                << minSnapDist[pointI]
+//                << " metre." << endl;
+//        }
+//    }
+//
+//    {
+//        scalarField magDisp(mag(patchDisp));
+//
+//        Info<< "Wanted displacement : average:"
+//            << gSum(magDisp)/returnReduce(patchDisp.size(), sumOp<label>())
+//            << " min:" << gMin(magDisp)
+//            << " max:" << gMax(magDisp) << endl;
+//    }
+//
+//
+//    Info<< "Calculated surface displacement in = "
+//        << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
+//
+//
+//    // Limit amount of movement.
+//    forAll(patchDisp, patchPointI)
+//    {
+//        scalar magDisp = mag(patchDisp[patchPointI]);
+//
+//        if (magDisp > snapDist[patchPointI])
+//        {
+//            patchDisp[patchPointI] *= snapDist[patchPointI] / magDisp;
+//
+//            Pout<< "Limiting displacement for " << patchPointI
+//                << " from " << magDisp << " to " << snapDist[patchPointI]
+//                << endl;
+//        }
+//    }
+//
+//    // Points on zones in one domain but only present as point on other
+//    // will not do condition 2 on all. Sync explicitly.
+//    syncTools::syncPointList
+//    (
+//        mesh,
+//        pp.meshPoints(),
+//        patchDisp,
+//        minMagSqrEqOp<point>(),    // combine op
+//        vector(GREAT, GREAT, GREAT)// null value (note: cannot use VGREAT)
+//    );
+//
+//    return patchDisp;
+//}
+////XXXXXXXXX
+
 void Foam::autoSnapDriver::smoothDisplacement
 (
     const snapParameters& snapParams,
@@ -1153,7 +1402,15 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
         scalarField faceSnapDist(pp.size(), -GREAT);
         {
             // Distance to attract to nearest feature on surface
-            const scalarField snapDist(calcSnapDistance(snapParams, pp));
+            const scalarField snapDist
+            (
+                calcSnapDistance
+                (
+                    mesh,
+                    snapParams,
+                    pp
+                )
+            );
 
             const faceList& localFaces = pp.localFaces();
 
@@ -1429,7 +1686,7 @@ void Foam::autoSnapDriver::doSnap
         indirectPrimitivePatch& pp = ppPtr();
 
         // Distance to attract to nearest feature on surface
-        const scalarField snapDist(calcSnapDistance(snapParams, pp));
+        const scalarField snapDist(calcSnapDistance(mesh, snapParams, pp));
 
 
         // Construct iterative mesh mover.
@@ -1467,7 +1724,14 @@ void Foam::autoSnapDriver::doSnap
             << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
 
         // Pre-smooth patch vertices (so before determining nearest)
-        preSmoothPatch(snapParams, nInitErrors, baffles, meshMover);
+        preSmoothPatch
+        (
+            meshRefiner_,
+            snapParams,
+            nInitErrors,
+            baffles,
+            meshMover
+        );
 
 
         for (label iter = 0; iter < nFeatIter; iter++)
@@ -1478,7 +1742,12 @@ void Foam::autoSnapDriver::doSnap
 
             // Calculate displacement at every patch point. Insert into
             // meshMover.
-            vectorField disp = calcNearestSurface(snapDist, meshMover);
+            vectorField disp = calcNearestSurface
+            (
+                meshRefiner_,
+                snapDist,
+                meshMover.patch()
+            );
 
             // 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 e6763bc084..f93ef62e08 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.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
@@ -81,11 +81,11 @@ class autoSnapDriver
 
             //- Calculate displacement per patch point to smooth out patch.
             //  Quite complicated in determining which points to move where.
-            pointField smoothPatchDisplacement
+            static pointField smoothPatchDisplacement
             (
                 const motionSmoother&,
                 const List<labelPair>&
-            ) const;
+            );
 
             //- Check that face zones are synced
             void checkCoupledFaceZones() const;
@@ -406,37 +406,51 @@ public:
             autoPtr<mapPolyMesh> mergeZoneBaffles(const List<labelPair>&);
 
             //- Calculate edge length per patch point.
-            scalarField calcSnapDistance
+            static scalarField calcSnapDistance
             (
+                const fvMesh& mesh,
                 const snapParameters& snapParams,
                 const indirectPrimitivePatch&
-            ) const;
+            );
 
             //- Smooth the mesh (patch and internal) to increase visibility
             //  of surface points (on castellated mesh) w.r.t. surface.
-            void preSmoothPatch
+            static void preSmoothPatch
             (
+                const meshRefinement& meshRefiner,
                 const snapParameters& snapParams,
                 const label nInitErrors,
                 const List<labelPair>& baffles,
                 motionSmoother&
-            ) const;
+            );
 
             //- Get points both on patch and facezone.
-            labelList getZoneSurfacePoints
+            static labelList getZoneSurfacePoints
             (
+                const fvMesh& mesh,
                 const indirectPrimitivePatch&,
                 const word& zoneName
-            ) const;
+            );
 
             //- Per patch point calculate point on nearest surface. Set as
             //  boundary conditions of motionSmoother displacement field. Return
             //  displacement of patch points.
-            vectorField calcNearestSurface
+            static vectorField calcNearestSurface
             (
+                const meshRefinement& meshRefiner,
                 const scalarField& snapDist,
-                motionSmoother& meshMover
-            ) const;
+                const indirectPrimitivePatch&
+            );
+
+            ////- Per patch point calculate point on nearest surface. Set as
+            ////  boundary conditions of motionSmoother displacement field.
+            ////  Return displacement of patch points.
+            //static vectorField calcNearestLocalSurface
+            //(
+            //    const meshRefinement& meshRefiner,
+            //    const scalarField& snapDist,
+            //    const indirectPrimitivePatch&
+            //);
 
             //- Smooth the displacement field to the internal.
             void smoothDisplacement
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
index 7bf8d8cf6c..0e2fcd07ce 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
@@ -45,7 +45,11 @@ Foam::refinementParameters::refinementParameters
     nBufferLayers_(readLabel(dict.lookup("nBufferLayers"))),
     keepPoints_(dict.lookup("keepPoints")),
     allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")),
-    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0))
+    useTopologicalSnapDetection_
+    (
+        dict.lookupOrDefault<bool>("useTopologicalSnapDetection", true)
+    ),
+    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance", 0))
 {}
 
 
@@ -65,7 +69,11 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
     nBufferLayers_(readLabel(dict.lookup("nCellsBetweenLevels"))),
     keepPoints_(pointField(1, dict.lookup("locationInMesh"))),
     allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")),
-    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0))
+    useTopologicalSnapDetection_
+    (
+        dict.lookupOrDefault<bool>("useTopologicalSnapDetection", true)
+    ),
+    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance", 0))
 {
     scalar featAngle(readScalar(dict.lookup("resolveFeatureAngle")));
 
@@ -83,7 +91,7 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::labelList Foam::refinementParameters::findCells(const polyMesh& mesh)
- const
+const
 {
     // Force calculation of tet-diag decomposition (for use in findCell)
     (void)mesh.tetBasePtIs();
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H
index 4fb05644e8..427fbd64b1 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H
@@ -80,6 +80,10 @@ class refinementParameters
         //  cellZone?
         Switch allowFreeStandingZoneFaces_;
 
+        //- Use old topology based problem-cell removal (cells with 8 points
+        //  on surface)
+        Switch useTopologicalSnapDetection_;
+
         //- Allowed load unbalance
         scalar maxLoadUnbalance_;
 
@@ -157,6 +161,13 @@ public:
                 return allowFreeStandingZoneFaces_;
             }
 
+            //- Use old topology based problem-cell removal
+            //  (cells with 8 points on surface)
+            bool useTopologicalSnapDetection() const
+            {
+                return useTopologicalSnapDetection_;
+            }
+
             //- Allowed load unbalance
             scalar maxLoadUnbalance() const
             {
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 7487e8a2bf..09708a3d1a 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -72,6 +72,8 @@ class globalIndex;
 class removePoints;
 class localPointRegion;
 
+class snapParameters;
+
 /*---------------------------------------------------------------------------*\
                            Class meshRefinement Declaration
 \*---------------------------------------------------------------------------*/
@@ -470,6 +472,14 @@ private:
                 const labelList& globalToMasterPatch
             ) const;
 
+            //- Returns list with for every internal face -1 or the patch
+            //  they should be baffled into.
+            labelList markFacesOnProblemCellsGeometric
+            (
+                const snapParameters& snapParams,
+                const dictionary& motionDict
+            ) const;
+
 
         // Baffle merging
 
@@ -536,6 +546,8 @@ private:
             //- Remove any loose standing cells
             void handleSnapProblems
             (
+                const snapParameters& snapParams,
+                const bool useTopologicalSnapDetection,
                 const bool removeEdgeConnectedCells,
                 const scalarField& perpendicularAngle,
                 const dictionary& motionDict,
@@ -760,10 +772,17 @@ public:
             void baffleAndSplitMesh
             (
                 const bool handleSnapProblems,
+
+                // How to remove problem snaps
+                const snapParameters& snapParams,
+                const bool useTopologicalSnapDetection,
                 const bool removeEdgeConnectedCells,
                 const scalarField& perpendicularAngle,
+
+                // How to handle free-standing baffles
                 const bool mergeFreeStanding,
                 const scalar freeStandingAngle,
+
                 const dictionary& motionDict,
                 Time& runTime,
                 const labelList& globalToMasterPatch,
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 628dd7898a..6b9d1074c5 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -1872,6 +1872,8 @@ void Foam::meshRefinement::makeConsistentFaceIndex
 
 void Foam::meshRefinement::handleSnapProblems
 (
+    const snapParameters& snapParams,
+    const bool useTopologicalSnapDetection,
     const bool removeEdgeConnectedCells,
     const scalarField& perpendicularAngle,
     const dictionary& motionDict,
@@ -1885,44 +1887,39 @@ void Foam::meshRefinement::handleSnapProblems
         << "----------------------------------------------" << nl
         << endl;
 
-    labelList facePatch
-    (
-        markFacesOnProblemCells
+    labelList facePatch;
+    if (useTopologicalSnapDetection)
+    {
+        facePatch = markFacesOnProblemCells
         (
             motionDict,
             removeEdgeConnectedCells,
             perpendicularAngle,
             globalToMasterPatch
-        )
-    );
+        );
+    }
+    else
+    {
+        facePatch = markFacesOnProblemCellsGeometric(snapParams, motionDict);
+    }
     Info<< "Analyzed problem cells in = "
         << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
 
     if (debug&meshRefinement::MESH)
     {
-        faceSet problemTopo(mesh_, "problemFacesTopo", 100);
+        faceSet problemFaces(mesh_, "problemFaces", 100);
 
-        const labelList facePatchTopo
-        (
-            markFacesOnProblemCells
-            (
-                motionDict,
-                removeEdgeConnectedCells,
-                perpendicularAngle,
-                globalToMasterPatch
-            )
-        );
-        forAll(facePatchTopo, faceI)
+        forAll(facePatch, faceI)
         {
-            if (facePatchTopo[faceI] != -1)
+            if (facePatch[faceI] != -1)
             {
-                problemTopo.insert(faceI);
+                problemFaces.insert(faceI);
             }
         }
-        problemTopo.instance() = timeName();
-        Pout<< "Dumping " << problemTopo.size()
-            << " problem faces to " << problemTopo.objectPath() << endl;
-        problemTopo.write();
+        problemFaces.instance() = timeName();
+        Pout<< "Dumping " << problemFaces.size()
+            << " problem faces to " << problemFaces.objectPath() << endl;
+        problemFaces.write();
     }
 
     Info<< "Introducing baffles to delete problem cells." << nl << endl;
@@ -1962,6 +1959,8 @@ void Foam::meshRefinement::handleSnapProblems
 void Foam::meshRefinement::baffleAndSplitMesh
 (
     const bool doHandleSnapProblems,
+    const snapParameters& snapParams,
+    const bool useTopologicalSnapDetection,
     const bool removeEdgeConnectedCells,
     const scalarField& perpendicularAngle,
     const bool mergeFreeStanding,
@@ -2029,83 +2028,10 @@ void Foam::meshRefinement::baffleAndSplitMesh
 
     if (doHandleSnapProblems)
     {
-        //Info<< nl
-        //    << "Introducing baffles to block off problem cells" << nl
-        //    << "----------------------------------------------" << nl
-        //    << endl;
-        //
-        //labelList facePatch
-        //(
-        //    markFacesOnProblemCells
-        //    (
-        //        motionDict,
-        //        removeEdgeConnectedCells,
-        //        perpendicularAngle,
-        //        globalToMasterPatch
-        //    )
-        //    //markFacesOnProblemCellsGeometric(motionDict)
-        //);
-        //Info<< "Analyzed problem cells in = "
-        //    << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
-        //
-        //if (debug&meshRefinement::MESH)
-        //{
-        //    faceSet problemTopo(mesh_, "problemFacesTopo", 100);
-        //
-        //    const labelList facePatchTopo
-        //    (
-        //        markFacesOnProblemCells
-        //        (
-        //            motionDict,
-        //            removeEdgeConnectedCells,
-        //            perpendicularAngle,
-        //            globalToMasterPatch
-        //        )
-        //    );
-        //    forAll(facePatchTopo, faceI)
-        //    {
-        //        if (facePatchTopo[faceI] != -1)
-        //        {
-        //            problemTopo.insert(faceI);
-        //        }
-        //    }
-        //    problemTopo.instance() = timeName();
-        //    Pout<< "Dumping " << problemTopo.size()
-        //        << " problem faces to " << problemTopo.objectPath() << endl;
-        //    problemTopo.write();
-        //}
-        //
-        //Info<< "Introducing baffles to delete problem cells." << nl << endl;
-        //
-        //if (debug)
-        //{
-        //    runTime++;
-        //}
-        //
-        //// Create baffles with same owner and neighbour for now.
-        //createBaffles(facePatch, facePatch);
-        //
-        //if (debug)
-        //{
-        //    // Debug:test all is still synced across proc patches
-        //    checkData();
-        //}
-        //Info<< "Created baffles in = "
-        //    << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
-        //
-        //printMeshInfo(debug, "After introducing baffles");
-        //
-        //if (debug&meshRefinement::MESH)
-        //{
-        //    Pout<< "Writing extra baffled mesh to time "
-        //        << timeName() << endl;
-        //    write(debug, runTime.path()/"extraBaffles");
-        //    Pout<< "Dumped debug data in = "
-        //        << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
-        //}
-
         handleSnapProblems
         (
+            snapParams,
+            useTopologicalSnapDetection,
             removeEdgeConnectedCells,
             perpendicularAngle,
             motionDict,
@@ -2192,6 +2118,8 @@ void Foam::meshRefinement::baffleAndSplitMesh
             // and delete them
             handleSnapProblems
             (
+                snapParams,
+                useTopologicalSnapDetection,
                 removeEdgeConnectedCells,
                 perpendicularAngle,
                 motionDict,
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
index 8676847fd4..dcc268c238 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
@@ -36,6 +36,10 @@ License
 #include "polyMeshGeometry.H"
 #include "IOmanip.H"
 #include "unitConversion.H"
+#include "autoSnapDriver.H"
+
+#include "snapParameters.H"
+#include "motionSmoother.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -968,164 +972,241 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 }
 
 
-//// Mark faces to be baffled to prevent snapping problems. Does
-//// test to find nearest surface and checks which faces would get squashed.
-//Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
-//(
-//    const dictionary& motionDict
-//) const
-//{
-//    // Construct addressing engine.
-//    autoPtr<indirectPrimitivePatch> ppPtr
-//    (
-//        meshRefinement::makePatch
-//        (
-//            mesh_,
-//            meshedPatches()
-//        )
-//    );
-//    const indirectPrimitivePatch& pp = ppPtr();
-//    const pointField& localPoints = pp.localPoints();
-//    const labelList& meshPoints = pp.meshPoints();
-//
-//    // Find nearest (non-baffle) surface
-//    pointField newPoints(mesh_.points());
-//    {
-//        List<pointIndexHit> hitInfo;
-//        labelList hitSurface;
-//        surfaces_.findNearest
-//        (
-//            surfaces_.getUnnamedSurfaces(),
-//            localPoints,
-//            scalarField(localPoints.size(), sqr(GREAT)),// sqr of attraction
-//            hitSurface,
-//            hitInfo
-//        );
-//
-//        forAll(hitInfo, i)
-//        {
-//            if (hitInfo[i].hit())
-//            {
-//                //label pointI = meshPoints[i];
-//                //Pout<< "   " << pointI << " moved from "
-//                //    << mesh_.points()[pointI] << " by "
-//                //    << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
-//                //    << endl;
-//                newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
-//            }
-//        }
-//    }
-//
-//    // Per face (internal or coupled!) the patch that the
-//    // baffle should get (or -1).
-//    labelList facePatch(mesh_.nFaces(), -1);
-//    // Count of baffled faces
-//    label nBaffleFaces = 0;
-//
-//    {
-//        pointField oldPoints(mesh_.points());
-//        mesh_.movePoints(newPoints);
-//        faceSet wrongFaces(mesh_, "wrongFaces", 100);
-//        {
-//            //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
-//
-//            // Just check the errors from squashing
-//            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-//
-//            const labelList allFaces(identity(mesh_.nFaces()));
-//            label nWrongFaces = 0;
-//
-//            scalar minArea(readScalar(motionDict.lookup("minArea")));
-//            if (minArea > -SMALL)
-//            {
-//                polyMeshGeometry::checkFaceArea
-//                (
-//                    false,
-//                    minArea,
-//                    mesh_,
-//                    mesh_.faceAreas(),
-//                    allFaces,
-//                    &wrongFaces
-//                );
-//
-//                label nNewWrongFaces = returnReduce
-//                (
-//                    wrongFaces.size(),
-//                    sumOp<label>()
-//                );
-//
-//                Info<< "    faces with area < "
-//                    << setw(5) << minArea
-//                    << " m^2                            : "
-//                    << nNewWrongFaces-nWrongFaces << endl;
-//
-//                nWrongFaces = nNewWrongFaces;
-//            }
-//
-////            scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
-//            scalar minDet = 0.01;
-//            if (minDet > -1)
-//            {
-//                polyMeshGeometry::checkCellDeterminant
-//                (
-//                    false,
-//                    minDet,
-//                    mesh_,
-//                    mesh_.faceAreas(),
-//                    allFaces,
-//                    polyMeshGeometry::affectedCells(mesh_, allFaces),
-//                    &wrongFaces
-//                );
-//
-//                label nNewWrongFaces = returnReduce
-//                (
-//                    wrongFaces.size(),
-//                    sumOp<label>()
-//                );
-//
-//                Info<< "    faces on cells with determinant < "
-//                    << setw(5) << minDet << "                : "
-//                    << nNewWrongFaces-nWrongFaces << endl;
-//
-//                nWrongFaces = nNewWrongFaces;
-//            }
-//        }
-//
-//
-//        forAllConstIter(faceSet, wrongFaces, iter)
-//        {
-//            label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
-//
-//            if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
-//            {
-//                facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
-//                nBaffleFaces++;
-//
-//                //Pout<< "    " << iter.key()
-//                //    //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
-//                //    << " is destined for patch " << facePatch[iter.key()]
-//                //    << endl;
-//            }
-//        }
-//        // Restore points.
-//        mesh_.movePoints(oldPoints);
-//    }
-//
-//
-//    Info<< "markFacesOnProblemCellsGeometric : marked "
-//        << returnReduce(nBaffleFaces, sumOp<label>())
-//        << " additional internal and coupled faces"
-//        << " to be converted into baffles." << endl;
-//
-//    syncTools::syncFaceList
-//    (
-//        mesh_,
-//        facePatch,
-//        maxEqOp<label>()
-//    );
-//
-//    return facePatch;
-//}
+// Mark faces to be baffled to prevent snapping problems. Does
+// test to find nearest surface and checks which faces would get squashed.
+Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
+(
+    const snapParameters& snapParams,
+    const dictionary& motionDict
+) const
+{
+    pointField oldPoints(mesh_.points());
+
+    // Repeat (most of) autoSnapDriver::doSnap
+    {
+        labelList adaptPatchIDs(meshedPatches());
+
+        // Construct addressing engine.
+        autoPtr<indirectPrimitivePatch> ppPtr
+        (
+            meshRefinement::makePatch
+            (
+                mesh_,
+                adaptPatchIDs
+            )
+        );
+        indirectPrimitivePatch& pp = ppPtr();
+
+        // Distance to attract to nearest feature on surface
+        const scalarField snapDist
+        (
+            autoSnapDriver::calcSnapDistance(mesh_, snapParams, pp)
+        );
+
+
+        // Construct iterative mesh mover.
+        Info<< "Constructing mesh displacer ..." << endl;
+        Info<< "Using mesh parameters " << motionDict << nl << endl;
+
+        const pointMesh& pMesh = pointMesh::New(mesh_);
+
+        motionSmoother meshMover
+        (
+            mesh_,
+            pp,
+            adaptPatchIDs,
+            meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs)(),
+            motionDict
+        );
+
+
+        // Check initial mesh
+        Info<< "Checking initial mesh ..." << endl;
+        labelHashSet wrongFaces(mesh_.nFaces()/100);
+        motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
+        const label nInitErrors = returnReduce
+        (
+            wrongFaces.size(),
+            sumOp<label>()
+        );
+
+        Info<< "Detected " << nInitErrors << " illegal faces"
+            << " (concave, zero area or negative cell pyramid volume)"
+            << endl;
+
+
+        Info<< "Checked initial mesh in = "
+            << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
+
+        // Pre-smooth patch vertices (so before determining nearest)
+        autoSnapDriver::preSmoothPatch
+        (
+            *this,
+            snapParams,
+            nInitErrors,
+            List<labelPair>(0), //baffles
+            meshMover
+        );
+
+        const vectorField disp
+        (
+            autoSnapDriver::calcNearestSurface
+            (
+                *this,
+                snapDist,   // attraction
+                pp
+            )
+        );
+
+        const labelList& meshPoints = pp.meshPoints();
+
+        pointField newPoints(mesh_.points());
+        forAll(meshPoints, i)
+        {
+            newPoints[meshPoints[i]] += disp[i];
+        }
+        mesh_.movePoints(newPoints);
+    }
+
+
+    // Per face (internal or coupled!) the patch that the
+    // baffle should get (or -1).
+    labelList facePatch(mesh_.nFaces(), -1);
+    // Count of baffled faces
+    label nBaffleFaces = 0;
+
+    {
+        faceSet wrongFaces(mesh_, "wrongFaces", 100);
+        {
+            //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
+
+            // Just check the errors from squashing
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            const labelList allFaces(identity(mesh_.nFaces()));
+            label nWrongFaces = 0;
+
+            //const scalar minV(readScalar(motionDict.lookup("minVol", true)));
+            //if (minV > -GREAT)
+            //{
+            //    polyMeshGeometry::checkFacePyramids
+            //    (
+            //        false,
+            //        minV,
+            //        mesh_,
+            //        mesh_.cellCentres(),
+            //        mesh_.points(),
+            //        allFaces,
+            //        List<labelPair>(0),
+            //        &wrongFaces
+            //    );
+            //
+            //    label nNewWrongFaces = returnReduce
+            //    (
+            //        wrongFaces.size(),
+            //        sumOp<label>()
+            //    );
+            //
+            //    Info<< "    faces with pyramid volume < "
+            //        << setw(5) << minV
+            //        << " m^3                  : "
+            //        << nNewWrongFaces-nWrongFaces << endl;
+            //
+            //    nWrongFaces = nNewWrongFaces;
+            //}
+
+            scalar minArea(readScalar(motionDict.lookup("minArea")));
+            if (minArea > -SMALL)
+            {
+                polyMeshGeometry::checkFaceArea
+                (
+                    false,
+                    minArea,
+                    mesh_,
+                    mesh_.faceAreas(),
+                    allFaces,
+                    &wrongFaces
+                );
+
+                label nNewWrongFaces = returnReduce
+                (
+                    wrongFaces.size(),
+                    sumOp<label>()
+                );
+
+                Info<< "    faces with area < "
+                    << setw(5) << minArea
+                    << " m^2                            : "
+                    << nNewWrongFaces-nWrongFaces << endl;
+
+                nWrongFaces = nNewWrongFaces;
+            }
+
+            scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
+            if (minDet > -1)
+            {
+                polyMeshGeometry::checkCellDeterminant
+                (
+                    false,
+                    minDet,
+                    mesh_,
+                    mesh_.faceAreas(),
+                    allFaces,
+                    polyMeshGeometry::affectedCells(mesh_, allFaces),
+                    &wrongFaces
+                );
+
+                label nNewWrongFaces = returnReduce
+                (
+                    wrongFaces.size(),
+                    sumOp<label>()
+                );
+
+                Info<< "    faces on cells with determinant < "
+                    << setw(5) << minDet << "                : "
+                    << nNewWrongFaces-nWrongFaces << endl;
+
+                nWrongFaces = nNewWrongFaces;
+            }
+        }
+
+
+        forAllConstIter(faceSet, wrongFaces, iter)
+        {
+            label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
+
+            if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
+            {
+                facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
+                nBaffleFaces++;
+
+                //Pout<< "    " << iter.key()
+                //    //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
+                //    << " is destined for patch " << facePatch[iter.key()]
+                //    << endl;
+            }
+        }
+    }
+
+
+    // Restore points.
+    mesh_.movePoints(oldPoints);
+
+
+    Info<< "markFacesOnProblemCellsGeometric : marked "
+        << returnReduce(nBaffleFaces, sumOp<label>())
+        << " additional internal and coupled faces"
+        << " to be converted into baffles." << endl;
+
+    syncTools::syncFaceList
+    (
+        mesh_,
+        facePatch,
+        maxEqOp<label>()
+    );
+
+    return facePatch;
+}
 
 
 // ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index f1917ee6f9..93c4600f9c 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -77,7 +77,8 @@ const Foam::NamedEnum<Foam::refinementSurfaces::faceZoneType, 3>
 Foam::refinementSurfaces::refinementSurfaces
 (
     const searchableSurfaces& allGeometry,
-    const dictionary& surfacesDict
+    const dictionary& surfacesDict,
+    const label gapLevelIncrement
 )
 :
     allGeometry_(allGeometry),
@@ -143,7 +144,7 @@ Foam::refinementSurfaces::refinementSurfaces
             globalLevelIncr[surfI] = dict.lookupOrDefault
             (
                 "gapLevelIncrement",
-                0
+                gapLevelIncrement
             );
 
             if
@@ -274,7 +275,7 @@ Foam::refinementSurfaces::refinementSurfaces
                         label levelIncr = regionDict.lookupOrDefault
                         (
                             "gapLevelIncrement",
-                            0
+                            gapLevelIncrement
                         );
                         regionLevelIncr[surfI].insert(regionI, levelIncr);
 
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
index c9773398e5..05ecdc8323 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -149,7 +149,8 @@ public:
         refinementSurfaces
         (
             const searchableSurfaces& allGeometry,
-            const dictionary&
+            const dictionary&,
+            const label gapLevelIncrement
         );