diff --git a/src/autoMesh/Make/files b/src/autoMesh/Make/files
index 96141445be..fecec5db5d 100644
--- a/src/autoMesh/Make/files
+++ b/src/autoMesh/Make/files
@@ -15,6 +15,7 @@ $(autoHexMeshDriver)/pointData/pointData.C
 $(autoHexMesh)/meshRefinement/meshRefinementBaffles.C
 $(autoHexMesh)/meshRefinement/meshRefinement.C
 $(autoHexMesh)/meshRefinement/meshRefinementMerge.C
+$(autoHexMesh)/meshRefinement/meshRefinementProblemCells.C
 $(autoHexMesh)/meshRefinement/meshRefinementRefine.C
 $(autoHexMesh)/refinementSurfaces/refinementSurfaces.C
 $(autoHexMesh)/shellSurfaces/shellSurfaces.C
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index ed9fd05808..fc70f98ccc 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -518,7 +518,9 @@ void Foam::autoRefineDriver::baffleAndSplitMesh
     // be like boundary face from now on so not coupled anymore.
     meshRefiner_.baffleAndSplitMesh
     (
-        handleSnapProblems,
+        handleSnapProblems,             // detect&remove potential snap problem
+        false,                          // perpendicular edge connected cells
+        scalarField(0),                 // per region perpendicular angle
         !handleSnapProblems,            // merge free standing baffles?
         motionDict,
         const_cast<Time&>(mesh.time()),
@@ -592,10 +594,14 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
         const_cast<Time&>(mesh.time())++;
     }
 
+    const scalarField& perpAngle = meshRefiner_.surfaces().perpendicularAngle();
+
     meshRefiner_.baffleAndSplitMesh
     (
         handleSnapProblems,
-        false,                  // merge free standing baffles?
+        handleSnapProblems,                 // remove perp edge connected cells
+        perpAngle,                          // perp angle
+        false,                              // merge free standing baffles?
         motionDict,
         const_cast<Time&>(mesh.time()),
         globalToPatch_,
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index da533c4d13..5895f5ae3a 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -36,6 +36,7 @@ SourceFiles
     meshRefinement.C
     meshRefinementBaffles.C
     meshRefinementMerge.C
+    meshRefinementProblemCells.C
     meshRefinementRefine.C
 
 \*---------------------------------------------------------------------------*/
@@ -51,6 +52,7 @@ SourceFiles
 #include "indirectPrimitivePatch.H"
 #include "pointFieldsFwd.H"
 #include "Tuple2.H"
+#include "pointIndexHit.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -345,6 +347,8 @@ private:
                 polyTopoChange& meshMod
             ) const;
 
+        // Problem cell handling
+
             //- Helper function to mark face as being on 'boundary'. Used by
             //  markFacesOnProblemCells
             void markBoundaryFace
@@ -355,15 +359,32 @@ private:
                 boolList& isBoundaryPoint
             ) const;
 
-            //- Returns list with for every internal face -1 or the patch
-            //  they should be baffled into.
-            labelList markFacesOnProblemCells
+            void findNearest
             (
-                const labelList& globalToPatch
+                const labelList& meshFaces,
+                List<pointIndexHit>& nearestInfo,
+                labelList& nearestSurface,
+                labelList& nearestRegion,
+                vectorField& nearestNormal
+            ) const;
+
+            Map<label> findEdgeConnectedProblemCells
+            (
+                const scalarField& perpendicularAngle,
+                const labelList&
             ) const;
 
             //- Returns list with for every internal face -1 or the patch
-            //  they should be baffled into.
+            //  they should be baffled into. If removeEdgeConnectedCells is set
+            //  removes cells based on perpendicularAngle.
+            labelList markFacesOnProblemCells
+            (
+                const bool removeEdgeConnectedCells,
+                const scalarField& perpendicularAngle,
+                const labelList& globalToPatch
+            ) const;
+
+            //- Initial test of marking faces using geometric information.
             labelList markFacesOnProblemCellsGeometric
             (
                 const dictionary& motionDict,
@@ -589,6 +610,8 @@ public:
             void baffleAndSplitMesh
             (
                 const bool handleSnapProblems,
+                const bool removeEdgeConnectedCells,
+                const scalarField& perpendicularAngle,
                 const bool mergeFreeStanding,
                 const dictionary& motionDict,
                 Time& runTime,
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 950cc7251b..6c5315a12a 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -44,10 +44,6 @@ License
 #include "OFstream.H"
 #include "regionSplit.H"
 #include "removeCells.H"
-#include "motionSmoother.H"
-#include "polyMeshGeometry.H"
-#include "IOmanip.H"
-#include "cellSet.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -230,7 +226,7 @@ void Foam::meshRefinement::getBafflePatches
     label vertI = 0;
     if (debug&OBJINTERSECTIONS)
     {
-        str.reset(new OFstream(mesh_.time().path()/"intersections.obj"));
+        str.reset(new OFstream(mesh_.time().timePath()/"intersections.obj"));
 
         Pout<< "getBafflePatches : Writing surface intersections to file "
             << str().name() << nl << endl;
@@ -510,642 +506,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
 }
 
 
-void Foam::meshRefinement::markBoundaryFace
-(
-    const label faceI,
-    boolList& isBoundaryFace,
-    boolList& isBoundaryEdge,
-    boolList& isBoundaryPoint
-) const
-{
-    isBoundaryFace[faceI] = true;
-
-    const labelList& fEdges = mesh_.faceEdges(faceI);
-
-    forAll(fEdges, fp)
-    {
-        isBoundaryEdge[fEdges[fp]] = true;
-    }
-
-    const face& f = mesh_.faces()[faceI];
-
-    forAll(f, fp)
-    {
-        isBoundaryPoint[f[fp]] = true;
-    }
-}
-
-
-// Returns list with for every internal face -1 or the patch they should
-// be baffled into. Gets run after createBaffles so all the surface
-// intersections have already been turned into baffles. Used to remove cells
-// by baffling all their faces and have the splitMeshRegions chuck away non
-// used regions.
-Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
-(
-    const labelList& globalToPatch
-) const
-{
-    const labelList& cellLevel = meshCutter_.cellLevel();
-    const labelList& pointLevel = meshCutter_.pointLevel();
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-
-    // Per internal face (boundary faces not used) the patch that the
-    // baffle should get (or -1)
-    labelList facePatch(mesh_.nFaces(), -1);
-
-    // Mark all points and edges on baffle patches (so not on any inlets,
-    // outlets etc.)
-    boolList isBoundaryPoint(mesh_.nPoints(), false);
-    boolList isBoundaryEdge(mesh_.nEdges(), false);
-    boolList isBoundaryFace(mesh_.nFaces(), false);
-
-    // Fill boundary data. All elements on meshed patches get marked.
-    // Get the labels of added patches.
-    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
-
-    forAll(adaptPatchIDs, i)
-    {
-        label patchI = adaptPatchIDs[i];
-
-        const polyPatch& pp = patches[patchI];
-
-        label faceI = pp.start();
-
-        forAll(pp, j)
-        {
-            markBoundaryFace
-            (
-                faceI,
-                isBoundaryFace,
-                isBoundaryEdge,
-                isBoundaryPoint
-            );
-
-            faceI++;
-        }
-    }
-
-    syncTools::syncPointList
-    (
-        mesh_,
-        isBoundaryPoint,
-        orEqOp<bool>(),
-        false,              // null value
-        false               // no separation
-    );
-
-    syncTools::syncEdgeList
-    (
-        mesh_,
-        isBoundaryEdge,
-        orEqOp<bool>(),
-        false,              // null value
-        false               // no separation
-    );
-
-    syncTools::syncFaceList
-    (
-        mesh_,
-        isBoundaryFace,
-        orEqOp<bool>(),
-        false               // no separation
-    );
-
-
-    // For each cell count the number of anchor points that are on
-    // the boundary:
-    // 8 : check the number of (baffle) boundary faces. If 3 or more block
-    //     off the cell since the cell would get squeezed down to a diamond
-    //     (probably; if the 3 or more faces are unrefined (only use the
-    //      anchor points))
-    // 7 : store. Used to check later on whether there are points with
-    //     3 or more of these cells. (note that on a flat surface a boundary
-    //     point will only have 4 cells connected to it)
-
-    // Does cell have exactly 7 of its 8 anchor points on the boundary?
-    PackedList<1> hasSevenBoundaryAnchorPoints(mesh_.nCells(), 0u);
-    // If so what is the remaining non-boundary anchor point?
-    labelHashSet nonBoundaryAnchors(mesh_.nCells()/10000);
-
-    // On-the-fly addressing storage.
-    DynamicList<label> dynFEdges;
-    DynamicList<label> dynCPoints;
-
-    // Count of faces marked for baffling
-    label nBaffleFaces = 0;
-
-    forAll(cellLevel, cellI)
-    {
-        const labelList& cPoints = mesh_.cellPoints(cellI, dynCPoints);
-
-        // Get number of anchor points (pointLevel == cellLevel)
-
-        label nBoundaryAnchors = 0;
-        label nNonAnchorBoundary = 0;
-        label nonBoundaryAnchor = -1;
-
-        forAll(cPoints, i)
-        {
-            label pointI = cPoints[i];
-
-            if (pointLevel[pointI] <= cellLevel[cellI])
-            {
-                // Anchor point
-                if (isBoundaryPoint[pointI])
-                {
-                    nBoundaryAnchors++;
-                }
-                else
-                {
-                    // Anchor point which is not on the surface
-                    nonBoundaryAnchor = pointI;
-                }
-            }
-            else if (isBoundaryPoint[pointI])
-            {
-                nNonAnchorBoundary++;
-            }
-        }
-
-        if (nBoundaryAnchors == 8)
-        {
-            const cell& cFaces = mesh_.cells()[cellI];
-
-            // Count boundary faces.
-            label nBfaces = 0;
-
-            forAll(cFaces, cFaceI)
-            {
-                if (isBoundaryFace[cFaces[cFaceI]])
-                {
-                    nBfaces++;
-                }
-            }
-
-            // If nBfaces > 1 make all non-boundary non-baffle faces baffles.
-            // We assume that this situation is where there is a single
-            // cell sticking out which would get flattened.
-
-            // Eugene: delete cell no matter what.
-            //if (nBfaces > 1)
-            {
-                forAll(cFaces, cf)
-                {
-                    label faceI = cFaces[cf];
-
-                    if (facePatch[faceI] == -1 && mesh_.isInternalFace(faceI))
-                    {
-                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
-                        nBaffleFaces++;
-
-                        // Mark face as a 'boundary'
-                        markBoundaryFace
-                        (
-                            faceI,
-                            isBoundaryFace,
-                            isBoundaryEdge,
-                            isBoundaryPoint
-                        );
-                    }
-                }
-            }
-        }
-        else if (nBoundaryAnchors == 7)
-        {
-            // Mark the cell. Store the (single!) non-boundary anchor point.
-            hasSevenBoundaryAnchorPoints.set(cellI, 1u);
-            nonBoundaryAnchors.insert(nonBoundaryAnchor);
-        }
-    }
-
-
-    // Loop over all points. If a point is connected to 4 or more cells
-    // with 7 anchor points on the boundary set those cell's non-boundary faces
-    // to baffles
-
-    DynamicList<label> dynPCells;
-
-    forAllConstIter(labelHashSet, nonBoundaryAnchors, iter)
-    {
-        label pointI = iter.key();
-
-        const labelList& pCells = mesh_.pointCells(pointI, dynPCells);
-
-        // Count number of 'hasSevenBoundaryAnchorPoints' cells.
-        label n = 0;
-
-        forAll(pCells, i)
-        {
-            if (hasSevenBoundaryAnchorPoints.get(pCells[i]) == 1u)
-            {
-                n++;
-            }
-        }
-
-        if (n > 3)
-        {
-            // Point in danger of being what? Remove all 7-cells.
-            forAll(pCells, i)
-            {
-                label cellI = pCells[i];
-
-                if (hasSevenBoundaryAnchorPoints.get(cellI) == 1u)
-                {
-                    const cell& cFaces = mesh_.cells()[cellI];
-
-                    forAll(cFaces, cf)
-                    {
-                        label faceI = cFaces[cf];
-
-                        if
-                        (
-                            facePatch[faceI] == -1
-                         && mesh_.isInternalFace(faceI)
-                        )
-                        {
-                            facePatch[faceI] = getBafflePatch(facePatch, faceI);
-                            nBaffleFaces++;
-
-                            // Mark face as a 'boundary'
-                            markBoundaryFace
-                            (
-                                faceI,
-                                isBoundaryFace,
-                                isBoundaryEdge,
-                                isBoundaryPoint
-                            );
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-
-    // Sync all. (note that pointdata and facedata not used anymore but sync
-    // anyway)
-
-    syncTools::syncPointList
-    (
-        mesh_,
-        isBoundaryPoint,
-        orEqOp<bool>(),
-        false,              // null value
-        false               // no separation
-    );
-
-    syncTools::syncEdgeList
-    (
-        mesh_,
-        isBoundaryEdge,
-        orEqOp<bool>(),
-        false,              // null value
-        false               // no separation
-    );
-
-    syncTools::syncFaceList
-    (
-        mesh_,
-        isBoundaryFace,
-        orEqOp<bool>(),
-        false               // no separation
-    );
-
-
-    // Find faces with all edges on the boundary and make them baffles
-    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
-    {
-        if (facePatch[faceI] == -1)
-        {
-            const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
-            label nFaceBoundaryEdges = 0;
-
-            forAll(fEdges, fe)
-            {
-                if (isBoundaryEdge[fEdges[fe]])
-                {
-                    nFaceBoundaryEdges++;
-                }
-            }
-
-            if (nFaceBoundaryEdges == fEdges.size())
-            {
-                facePatch[faceI] = getBafflePatch(facePatch, faceI);
-                nBaffleFaces++;
-
-                // Do NOT update boundary data since this would grow blocked
-                // faces across gaps.
-            }
-        }
-    }
-
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (pp.coupled())
-        {
-            label faceI = pp.start();
-
-            forAll(pp, i)
-            {
-                if (facePatch[faceI] == -1)
-                {
-                    const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
-                    label nFaceBoundaryEdges = 0;
-
-                    forAll(fEdges, fe)
-                    {
-                        if (isBoundaryEdge[fEdges[fe]])
-                        {
-                            nFaceBoundaryEdges++;
-                        }
-                    }
-
-                    if (nFaceBoundaryEdges == fEdges.size())
-                    {
-                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
-                        nBaffleFaces++;
-
-                        // Do NOT update boundary data since this would grow
-                        // blocked faces across gaps.
-                    }
-                }
-
-                faceI++;
-            }
-        }
-    }
-
-    Info<< "markFacesOnProblemCells : marked "
-        << returnReduce(nBaffleFaces, sumOp<label>())
-        << " additional internal faces to be converted into baffles."
-        << endl;
-
-    return facePatch;
-}
-//XXXXXXXXXXXXXX
-// 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 labelList& globalToPatch
-) const
-{
-    // Get the labels of added patches.
-    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
-
-    // Construct addressing engine.
-    autoPtr<indirectPrimitivePatch> ppPtr
-    (
-        meshRefinement::makePatch
-        (
-            mesh_,
-            adaptPatchIDs
-        )
-    );
-    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;
-
-
-//    // Sync position? Or not since same face on both side so just sync
-//    // result of baffle.
-//
-//    const scalar minArea(readScalar(motionDict.lookup("minArea")));
-//
-//    Pout<< "markFacesOnProblemCellsGeometric : Comparing to minArea:"
-//        << minArea << endl;
-//
-//    pointField facePoints;
-//    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
-//    {
-//        const face& f = mesh_.faces()[faceI];
-//
-//        bool usesPatchPoint = false;
-//
-//        facePoints.setSize(f.size());
-//        forAll(f, fp)
-//        {
-//            Map<label>::const_iterator iter = pp.meshPointMap().find(f[fp]);
-//
-//            if (iter != pp.meshPointMap().end())
-//            {
-//                facePoints[fp] = newPosition[iter()];
-//                usesPatchPoint = true;
-//            }
-//            else
-//            {
-//                facePoints[fp] = mesh_.points()[f[fp]];
-//            }
-//        }
-//
-//        if (usesPatchPoint)
-//        {
-//            // Check area of face wrt original area
-//            face identFace(identity(f.size()));
-//
-//            if (identFace.mag(facePoints) < minArea)
-//            {
-//                facePatch[faceI] = getBafflePatch(facePatch, faceI);
-//                nBaffleFaces++;
-//            }
-//        }
-//    }
-//
-//
-//    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-//    forAll(patches, patchI)
-//    {
-//        const polyPatch& pp = patches[patchI];
-//
-//        if (pp.coupled())
-//        {
-//            forAll(pp, i)
-//            {
-//                label faceI = pp.start()+i;
-//
-//                const face& f = mesh_.faces()[faceI];
-//
-//                bool usesPatchPoint = false;
-//
-//                facePoints.setSize(f.size());
-//                forAll(f, fp)
-//                {
-//                    Map<label>::const_iterator iter =
-//                        pp.meshPointMap().find(f[fp]);
-//
-//                    if (iter != pp.meshPointMap().end())
-//                    {
-//                        facePoints[fp] = newPosition[iter()];
-//                        usesPatchPoint = true;
-//                    }
-//                    else
-//                    {
-//                        facePoints[fp] = mesh_.points()[f[fp]];
-//                    }
-//                }
-//
-//                if (usesPatchPoint)
-//                {
-//                    // Check area of face wrt original area
-//                    face identFace(identity(f.size()));
-//
-//                    if (identFace.mag(facePoints) < minArea)
-//                    {
-//                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
-//                        nBaffleFaces++;
-//                    }
-//                }
-//            }
-//        }
-//    }
-
-    {
-        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>(),
-        false               // no separation
-    );
-
-    return facePatch;
-}
-//XXXXXXXX
-
-
 // Return a list of coupled face pairs, i.e. faces that use the same vertices.
 // (this information is recalculated instead of maintained since would be too
 // hard across splitMeshRegions).
@@ -1855,13 +1215,13 @@ void Foam::meshRefinement::findCellZoneTopo
             break;
         }
 
-	// Synchronise regionToCellZone.
-	// Note:
-	// - region numbers are identical on all processors
-	// - keepRegion is identical ,,
-	// - cellZones are identical ,,
-	Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
-	Pstream::listCombineScatter(regionToCellZone);
+        // Synchronise regionToCellZone.
+        // Note:
+        // - region numbers are identical on all processors
+        // - keepRegion is identical ,,
+        // - cellZones are identical ,,
+        Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
+        Pstream::listCombineScatter(regionToCellZone);
     }
 
 
@@ -1904,6 +1264,8 @@ void Foam::meshRefinement::findCellZoneTopo
 void Foam::meshRefinement::baffleAndSplitMesh
 (
     const bool handleSnapProblems,
+    const bool removeEdgeConnectedCells,
+    const scalarField& perpendicularAngle,
     const bool mergeFreeStanding,
     const dictionary& motionDict,
     Time& runTime,
@@ -1981,6 +1343,8 @@ void Foam::meshRefinement::baffleAndSplitMesh
         (
             markFacesOnProblemCells
             (
+                removeEdgeConnectedCells,
+                perpendicularAngle,
                 globalToPatch
             )
             //markFacesOnProblemCellsGeometric
@@ -2024,6 +1388,8 @@ void Foam::meshRefinement::baffleAndSplitMesh
             (    
                 markFacesOnProblemCells
                 (
+                    removeEdgeConnectedCells,
+                    perpendicularAngle,
                     globalToPatch
                 )
             );
@@ -2099,7 +1465,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
     {
         Pout<< "Writing subsetted mesh to time " << mesh_.time().timeName()
             << endl;
-        write(debug, runTime.path()/runTime.timeName());
+        write(debug, runTime.timePath());
         Pout<< "Dumped debug data in = "
             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
     }
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
new file mode 100644
index 0000000000..8a027fa654
--- /dev/null
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
@@ -0,0 +1,952 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "meshRefinement.H"
+#include "fvMesh.H"
+#include "syncTools.H"
+#include "Time.H"
+#include "refinementSurfaces.H"
+#include "pointSet.H"
+#include "faceSet.H"
+#include "indirectPrimitivePatch.H"
+#include "OFstream.H"
+#include "cellSet.H"
+#include "searchableSurfaces.H"
+#include "polyMeshGeometry.H"
+#include "IOmanip.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::meshRefinement::markBoundaryFace
+(
+    const label faceI,
+    boolList& isBoundaryFace,
+    boolList& isBoundaryEdge,
+    boolList& isBoundaryPoint
+) const
+{
+    isBoundaryFace[faceI] = true;
+
+    const labelList& fEdges = mesh_.faceEdges(faceI);
+
+    forAll(fEdges, fp)
+    {
+        isBoundaryEdge[fEdges[fp]] = true;
+    }
+
+    const face& f = mesh_.faces()[faceI];
+
+    forAll(f, fp)
+    {
+        isBoundaryPoint[f[fp]] = true;
+    }
+}
+
+
+void Foam::meshRefinement::findNearest
+(
+    const labelList& meshFaces,
+    List<pointIndexHit>& nearestInfo,
+    labelList& nearestSurface,
+    labelList& nearestRegion,
+    vectorField& nearestNormal
+) const
+{
+    pointField fc(meshFaces.size());
+    forAll(meshFaces, i)
+    {
+        fc[i] = mesh_.faceCentres()[meshFaces[i]];
+    }
+
+    const labelList allSurfaces(identity(surfaces_.surfaces().size()));
+
+    surfaces_.findNearest
+    (
+        allSurfaces,
+        fc,
+        scalarField(fc.size(), sqr(GREAT)),    // sqr of attraction
+        nearestSurface,
+        nearestInfo
+    );
+
+    // Do normal testing per surface.
+    nearestNormal.setSize(nearestInfo.size());
+    nearestRegion.setSize(nearestInfo.size());
+
+    forAll(allSurfaces, surfI)
+    {
+        DynamicList<pointIndexHit> localHits;
+
+        forAll(nearestSurface, i)
+        {
+            if (nearestSurface[i] == surfI)
+            {
+                localHits.append(nearestInfo[i]);
+            }
+        }
+
+        label geomI = surfaces_.surfaces()[surfI];
+
+        pointField localNormals;
+        surfaces_.geometry()[geomI].getNormal(localHits, localNormals);
+
+        labelList localRegion;
+        surfaces_.geometry()[geomI].getRegion(localHits, localRegion);
+
+        label localI = 0;
+        forAll(nearestSurface, i)
+        {
+            if (nearestSurface[i] == surfI)
+            {
+                nearestNormal[i] = localNormals[localI];
+                nearestRegion[i] = localRegion[localI];
+                localI++;
+            }
+        }
+    }
+}
+
+
+Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
+(
+    const scalarField& perpendicularAngle,
+    const labelList& globalToPatch
+) const
+{
+    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
+
+    // Construct addressing engine.
+    autoPtr<indirectPrimitivePatch> ppPtr
+    (
+        meshRefinement::makePatch
+        (
+            mesh_,
+            adaptPatchIDs
+        )
+    );
+    const indirectPrimitivePatch& pp = ppPtr();
+
+
+    // 1. Collect faces to test
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    DynamicList<label> candidateFaces(pp.size()/20);
+
+    const labelListList& edgeFaces = pp.edgeFaces();
+
+    const labelList& cellLevel = meshCutter_.cellLevel();
+
+    forAll(edgeFaces, edgeI)
+    {
+        const labelList& eFaces = edgeFaces[edgeI];
+
+        if (eFaces.size() == 2)
+        {
+            label face0 = pp.addressing()[eFaces[0]];
+            label face1 = pp.addressing()[eFaces[1]];
+
+            label cell0 = mesh_.faceOwner()[face0];
+            label cell1 = mesh_.faceOwner()[face1];
+
+            if (cellLevel[cell0] > cellLevel[cell1])
+            {
+                // cell0 smaller.
+                const vector& n0 = pp.faceNormals()[eFaces[0]];
+                const vector& n1 = pp.faceNormals()[eFaces[1]];
+
+                if (mag(n0 & n1) < 0.1)
+                {
+                    candidateFaces.append(face0);
+                }
+            }
+            else if (cellLevel[cell1] > cellLevel[cell0])
+            {
+                // cell1 smaller.
+                const vector& n0 = pp.faceNormals()[eFaces[0]];
+                const vector& n1 = pp.faceNormals()[eFaces[1]];
+
+                if (mag(n0 & n1) < 0.1)
+                {
+                    candidateFaces.append(face1);
+                }
+            }
+        }
+    }
+    candidateFaces.shrink();
+
+    Info<< "Testing " << returnReduce(candidateFaces.size(), sumOp<label>())
+        << " faces on edge-connected cells of differing level."
+        << endl;
+
+    if (debug)
+    {
+        faceSet fSet(mesh_, "edgeConnectedFaces", candidateFaces);
+        Pout<< "Writing " << fSet.size()
+            << " with problematic topology to faceSet "
+            << fSet.objectPath() << endl;
+        fSet.write();
+    }
+
+
+    // 2. Find nearest surface on candidate faces
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    List<pointIndexHit> nearestInfo;
+    labelList nearestSurface;
+    labelList nearestRegion;
+    vectorField nearestNormal;
+    findNearest
+    (
+        candidateFaces,
+        nearestInfo,
+        nearestSurface,
+        nearestRegion,
+        nearestNormal
+    );
+
+
+    // 3. Test angle to surface
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    Map<label> candidateCells(candidateFaces.size());
+
+    faceSet perpFaces(mesh_, "perpendicularFaces", pp.size()/100);
+
+    forAll(candidateFaces, i)
+    {
+        label faceI = candidateFaces[i];
+
+        vector n = mesh_.faceAreas()[faceI];
+        n /= mag(n);
+
+        label region = surfaces_.globalRegion
+        (
+            nearestSurface[i],
+            nearestRegion[i]
+        );
+
+        scalar angle =
+            perpendicularAngle[region]
+          / 180.0
+          * mathematicalConstant::pi;
+
+        if (angle >= 0)
+        {
+            if (mag(n & nearestNormal[i]) < Foam::sin(angle))
+            {
+                perpFaces.insert(faceI);
+                candidateCells.insert
+                (
+                    mesh_.faceOwner()[faceI],
+                    globalToPatch[region]
+                );
+            }
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "Writing " << perpFaces.size()
+            << " faces that are perpendicular to the surface to set "
+            << perpFaces.objectPath() << endl;
+        perpFaces.write();
+    }
+    return candidateCells;
+}
+
+
+// Returns list with for every internal face -1 or the patch they should
+// be baffled into. Gets run after createBaffles so all the surface
+// intersections have already been turned into baffles. Used to remove cells
+// by baffling all their faces and have the splitMeshRegions chuck away non
+// used regions.
+Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
+(
+    const bool removeEdgeConnectedCells,
+    const scalarField& perpendicularAngle,
+    const labelList& globalToPatch
+) const
+{
+    const labelList& cellLevel = meshCutter_.cellLevel();
+    const labelList& pointLevel = meshCutter_.pointLevel();
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+
+    // Per internal face (boundary faces not used) the patch that the
+    // baffle should get (or -1)
+    labelList facePatch(mesh_.nFaces(), -1);
+
+    // Mark all points and edges on baffle patches (so not on any inlets,
+    // outlets etc.)
+    boolList isBoundaryPoint(mesh_.nPoints(), false);
+    boolList isBoundaryEdge(mesh_.nEdges(), false);
+    boolList isBoundaryFace(mesh_.nFaces(), false);
+
+    // Fill boundary data. All elements on meshed patches get marked.
+    // Get the labels of added patches.
+    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
+
+    forAll(adaptPatchIDs, i)
+    {
+        label patchI = adaptPatchIDs[i];
+
+        const polyPatch& pp = patches[patchI];
+
+        label faceI = pp.start();
+
+        forAll(pp, j)
+        {
+            markBoundaryFace
+            (
+                faceI,
+                isBoundaryFace,
+                isBoundaryEdge,
+                isBoundaryPoint
+            );
+
+            faceI++;
+        }
+    }
+
+    // Count of faces marked for baffling
+    label nBaffleFaces = 0;
+
+    if (removeEdgeConnectedCells && max(perpendicularAngle) >= 0)
+    {
+        Info<< "markFacesOnProblemCells :"
+            << " Checking for edge-connected cells of highly differing sizes."
+            << endl;
+
+        // Pick up the cells that need to be removed and (a guess for)
+        // the patch they should be patched with.
+        Map<label> problemCells
+        (
+            findEdgeConnectedProblemCells
+            (
+                perpendicularAngle,
+                globalToPatch
+            )
+        );
+
+        // Baffle all faces of cells that need to be removed
+        forAllConstIter(Map<label>, problemCells, iter)
+        {
+            const cell& cFaces = mesh_.cells()[iter.key()];
+
+            forAll(cFaces, i)
+            {
+                label faceI = cFaces[i];
+
+                if (facePatch[faceI] == -1 && mesh_.isInternalFace(faceI))
+                {
+                    facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                    nBaffleFaces++;
+
+                    // Mark face as a 'boundary'
+                    markBoundaryFace
+                    (
+                        faceI,
+                        isBoundaryFace,
+                        isBoundaryEdge,
+                        isBoundaryPoint
+                    );
+                }
+            }
+        }
+        Info<< "markFacesOnProblemCells : Marked "
+            << returnReduce(nBaffleFaces, sumOp<label>())
+            << " additional internal faces to be converted into baffles"
+            << " due to "
+            << returnReduce(problemCells.size(), sumOp<label>())
+            << " cells edge-connected to lower level cells." << endl;
+
+        if (debug)
+        {
+            cellSet problemCellSet(mesh_, "problemCells", problemCells.toc());
+            Pout<< "Writing " << problemCellSet.size()
+                << " cells that are edge connected to coarser cell to set "
+                << problemCellSet.objectPath() << endl;
+            problemCellSet.write();
+        }
+    }
+
+    syncTools::syncPointList
+    (
+        mesh_,
+        isBoundaryPoint,
+        orEqOp<bool>(),
+        false,              // null value
+        false               // no separation
+    );
+
+    syncTools::syncEdgeList
+    (
+        mesh_,
+        isBoundaryEdge,
+        orEqOp<bool>(),
+        false,              // null value
+        false               // no separation
+    );
+
+    syncTools::syncFaceList
+    (
+        mesh_,
+        isBoundaryFace,
+        orEqOp<bool>(),
+        false               // no separation
+    );
+
+
+    // For each cell count the number of anchor points that are on
+    // the boundary:
+    // 8 : check the number of (baffle) boundary faces. If 3 or more block
+    //     off the cell since the cell would get squeezed down to a diamond
+    //     (probably; if the 3 or more faces are unrefined (only use the
+    //      anchor points))
+    // 7 : store. Used to check later on whether there are points with
+    //     3 or more of these cells. (note that on a flat surface a boundary
+    //     point will only have 4 cells connected to it)
+
+    // Does cell have exactly 7 of its 8 anchor points on the boundary?
+    PackedList<1> hasSevenBoundaryAnchorPoints(mesh_.nCells(), 0u);
+    // If so what is the remaining non-boundary anchor point?
+    labelHashSet nonBoundaryAnchors(mesh_.nCells()/10000);
+
+    // On-the-fly addressing storage.
+    DynamicList<label> dynFEdges;
+    DynamicList<label> dynCPoints;
+
+    forAll(cellLevel, cellI)
+    {
+        const labelList& cPoints = mesh_.cellPoints(cellI, dynCPoints);
+
+        // Get number of anchor points (pointLevel == cellLevel)
+
+        label nBoundaryAnchors = 0;
+        label nNonAnchorBoundary = 0;
+        label nonBoundaryAnchor = -1;
+
+        forAll(cPoints, i)
+        {
+            label pointI = cPoints[i];
+
+            if (pointLevel[pointI] <= cellLevel[cellI])
+            {
+                // Anchor point
+                if (isBoundaryPoint[pointI])
+                {
+                    nBoundaryAnchors++;
+                }
+                else
+                {
+                    // Anchor point which is not on the surface
+                    nonBoundaryAnchor = pointI;
+                }
+            }
+            else if (isBoundaryPoint[pointI])
+            {
+                nNonAnchorBoundary++;
+            }
+        }
+
+        if (nBoundaryAnchors == 8)
+        {
+            const cell& cFaces = mesh_.cells()[cellI];
+
+            // Count boundary faces.
+            label nBfaces = 0;
+
+            forAll(cFaces, cFaceI)
+            {
+                if (isBoundaryFace[cFaces[cFaceI]])
+                {
+                    nBfaces++;
+                }
+            }
+
+            // If nBfaces > 1 make all non-boundary non-baffle faces baffles.
+            // We assume that this situation is where there is a single
+            // cell sticking out which would get flattened.
+
+            // Eugene: delete cell no matter what.
+            //if (nBfaces > 1)
+            {
+                forAll(cFaces, cf)
+                {
+                    label faceI = cFaces[cf];
+
+                    if (facePatch[faceI] == -1 && mesh_.isInternalFace(faceI))
+                    {
+                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                        nBaffleFaces++;
+
+                        // Mark face as a 'boundary'
+                        markBoundaryFace
+                        (
+                            faceI,
+                            isBoundaryFace,
+                            isBoundaryEdge,
+                            isBoundaryPoint
+                        );
+                    }
+                }
+            }
+        }
+        else if (nBoundaryAnchors == 7)
+        {
+            // Mark the cell. Store the (single!) non-boundary anchor point.
+            hasSevenBoundaryAnchorPoints.set(cellI, 1u);
+            nonBoundaryAnchors.insert(nonBoundaryAnchor);
+        }
+    }
+
+
+    // Loop over all points. If a point is connected to 4 or more cells
+    // with 7 anchor points on the boundary set those cell's non-boundary faces
+    // to baffles
+
+    DynamicList<label> dynPCells;
+
+    forAllConstIter(labelHashSet, nonBoundaryAnchors, iter)
+    {
+        label pointI = iter.key();
+
+        const labelList& pCells = mesh_.pointCells(pointI, dynPCells);
+
+        // Count number of 'hasSevenBoundaryAnchorPoints' cells.
+        label n = 0;
+
+        forAll(pCells, i)
+        {
+            if (hasSevenBoundaryAnchorPoints.get(pCells[i]) == 1u)
+            {
+                n++;
+            }
+        }
+
+        if (n > 3)
+        {
+            // Point in danger of being what? Remove all 7-cells.
+            forAll(pCells, i)
+            {
+                label cellI = pCells[i];
+
+                if (hasSevenBoundaryAnchorPoints.get(cellI) == 1u)
+                {
+                    const cell& cFaces = mesh_.cells()[cellI];
+
+                    forAll(cFaces, cf)
+                    {
+                        label faceI = cFaces[cf];
+
+                        if
+                        (
+                            facePatch[faceI] == -1
+                         && mesh_.isInternalFace(faceI)
+                        )
+                        {
+                            facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                            nBaffleFaces++;
+
+                            // Mark face as a 'boundary'
+                            markBoundaryFace
+                            (
+                                faceI,
+                                isBoundaryFace,
+                                isBoundaryEdge,
+                                isBoundaryPoint
+                            );
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+
+    // Sync all. (note that pointdata and facedata not used anymore but sync
+    // anyway)
+
+    syncTools::syncPointList
+    (
+        mesh_,
+        isBoundaryPoint,
+        orEqOp<bool>(),
+        false,              // null value
+        false               // no separation
+    );
+
+    syncTools::syncEdgeList
+    (
+        mesh_,
+        isBoundaryEdge,
+        orEqOp<bool>(),
+        false,              // null value
+        false               // no separation
+    );
+
+    syncTools::syncFaceList
+    (
+        mesh_,
+        isBoundaryFace,
+        orEqOp<bool>(),
+        false               // no separation
+    );
+
+
+    // Find faces with all edges on the boundary and make them baffles
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        if (facePatch[faceI] == -1)
+        {
+            const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
+            label nFaceBoundaryEdges = 0;
+
+            forAll(fEdges, fe)
+            {
+                if (isBoundaryEdge[fEdges[fe]])
+                {
+                    nFaceBoundaryEdges++;
+                }
+            }
+
+            if (nFaceBoundaryEdges == fEdges.size())
+            {
+                facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                nBaffleFaces++;
+
+                // Do NOT update boundary data since this would grow blocked
+                // faces across gaps.
+            }
+        }
+    }
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+
+            forAll(pp, i)
+            {
+                if (facePatch[faceI] == -1)
+                {
+                    const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
+                    label nFaceBoundaryEdges = 0;
+
+                    forAll(fEdges, fe)
+                    {
+                        if (isBoundaryEdge[fEdges[fe]])
+                        {
+                            nFaceBoundaryEdges++;
+                        }
+                    }
+
+                    if (nFaceBoundaryEdges == fEdges.size())
+                    {
+                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                        nBaffleFaces++;
+
+                        // Do NOT update boundary data since this would grow
+                        // blocked faces across gaps.
+                    }
+                }
+
+                faceI++;
+            }
+        }
+    }
+
+    Info<< "markFacesOnProblemCells : marked "
+        << returnReduce(nBaffleFaces, sumOp<label>())
+        << " additional internal faces to be converted into baffles."
+        << endl;
+
+    return facePatch;
+}
+//XXXXXXXXXXXXXX
+// 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 labelList& globalToPatch
+) const
+{
+    // Get the labels of added patches.
+    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
+
+    // Construct addressing engine.
+    autoPtr<indirectPrimitivePatch> ppPtr
+    (
+        meshRefinement::makePatch
+        (
+            mesh_,
+            adaptPatchIDs
+        )
+    );
+    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;
+
+
+//    // Sync position? Or not since same face on both side so just sync
+//    // result of baffle.
+//
+//    const scalar minArea(readScalar(motionDict.lookup("minArea")));
+//
+//    Pout<< "markFacesOnProblemCellsGeometric : Comparing to minArea:"
+//        << minArea << endl;
+//
+//    pointField facePoints;
+//    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
+//    {
+//        const face& f = mesh_.faces()[faceI];
+//
+//        bool usesPatchPoint = false;
+//
+//        facePoints.setSize(f.size());
+//        forAll(f, fp)
+//        {
+//            Map<label>::const_iterator iter = pp.meshPointMap().find(f[fp]);
+//
+//            if (iter != pp.meshPointMap().end())
+//            {
+//                facePoints[fp] = newPosition[iter()];
+//                usesPatchPoint = true;
+//            }
+//            else
+//            {
+//                facePoints[fp] = mesh_.points()[f[fp]];
+//            }
+//        }
+//
+//        if (usesPatchPoint)
+//        {
+//            // Check area of face wrt original area
+//            face identFace(identity(f.size()));
+//
+//            if (identFace.mag(facePoints) < minArea)
+//            {
+//                facePatch[faceI] = getBafflePatch(facePatch, faceI);
+//                nBaffleFaces++;
+//            }
+//        }
+//    }
+//
+//
+//    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+//    forAll(patches, patchI)
+//    {
+//        const polyPatch& pp = patches[patchI];
+//
+//        if (pp.coupled())
+//        {
+//            forAll(pp, i)
+//            {
+//                label faceI = pp.start()+i;
+//
+//                const face& f = mesh_.faces()[faceI];
+//
+//                bool usesPatchPoint = false;
+//
+//                facePoints.setSize(f.size());
+//                forAll(f, fp)
+//                {
+//                    Map<label>::const_iterator iter =
+//                        pp.meshPointMap().find(f[fp]);
+//
+//                    if (iter != pp.meshPointMap().end())
+//                    {
+//                        facePoints[fp] = newPosition[iter()];
+//                        usesPatchPoint = true;
+//                    }
+//                    else
+//                    {
+//                        facePoints[fp] = mesh_.points()[f[fp]];
+//                    }
+//                }
+//
+//                if (usesPatchPoint)
+//                {
+//                    // Check area of face wrt original area
+//                    face identFace(identity(f.size()));
+//
+//                    if (identFace.mag(facePoints) < minArea)
+//                    {
+//                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
+//                        nBaffleFaces++;
+//                    }
+//                }
+//            }
+//        }
+//    }
+
+    {
+        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>(),
+        false               // no separation
+    );
+
+    return facePatch;
+}
+//XXXXXXXX
+
+
+// ************************************************************************* //
diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index 2f23e9706c..276aa2fd01 100644
--- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -53,8 +53,10 @@ Foam::refinementSurfaces::refinementSurfaces
 {
     labelList globalMinLevel(surfaceDicts.size(), 0);
     labelList globalMaxLevel(surfaceDicts.size(), 0);
+    scalarField globalAngle(surfaceDicts.size(), -GREAT);
     List<Map<label> > regionMinLevel(surfaceDicts.size());
     List<Map<label> > regionMaxLevel(surfaceDicts.size());
+    List<Map<scalar> > regionAngle(surfaceDicts.size());
 
     //wordList globalPatchType(surfaceDicts.size());
     //List<HashTable<word> > regionPatchType(surfaceDicts.size());
@@ -80,6 +82,12 @@ Foam::refinementSurfaces::refinementSurfaces
             dict.lookup("zoneInside") >> zoneInside_[surfI];
         }
 
+        // Global perpendicular angle
+        if (dict.found("perpendicularAngle"))
+        {
+            globalAngle[surfI] = readScalar(dict.lookup("perpendicularAngle"));
+        }
+
         //// Global patch name per surface
         //if (dict.found("patchType"))
         //{
@@ -130,6 +138,15 @@ Foam::refinementSurfaces::refinementSurfaces
                         << exit(FatalError);
                 }
                 regionMaxLevel[surfI].insert(regionI, max);
+
+                if (regionDict.found("perpendicularAngle"))
+                {
+                    regionAngle[surfI].insert
+                    (
+                        regionI,
+                        readScalar(regionDict.lookup("perpendicularAngle"))
+                    );
+                }
             }
         }
     }
@@ -170,6 +187,8 @@ Foam::refinementSurfaces::refinementSurfaces
     minLevel_ = 0;
     maxLevel_.setSize(nRegions);
     maxLevel_ = 0;
+    perpendicularAngle_.setSize(nRegions);
+    perpendicularAngle_ = -GREAT;
     //patchName_.setSize(nRegions);
     //patchType_.setSize(nRegions);
 
@@ -182,6 +201,7 @@ Foam::refinementSurfaces::refinementSurfaces
         {
             minLevel_[regionOffset_[surfI] + i] = globalMinLevel[surfI];
             maxLevel_[regionOffset_[surfI] + i] = globalMaxLevel[surfI];
+            perpendicularAngle_[regionOffset_[surfI] + i] = globalAngle[surfI];
         }
 
         // Overwrite with region specific information
@@ -191,6 +211,7 @@ Foam::refinementSurfaces::refinementSurfaces
 
             minLevel_[globalRegionI] = iter();
             maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
+            perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
 
             // Check validity
             if
@@ -240,8 +261,10 @@ Foam::refinementSurfaces::refinementSurfaces
 {
     labelList globalMinLevel(surfacesDict.size(), 0);
     labelList globalMaxLevel(surfacesDict.size(), 0);
+    scalarField globalAngle(surfacesDict.size(), -GREAT);
     List<Map<label> > regionMinLevel(surfacesDict.size());
     List<Map<label> > regionMaxLevel(surfacesDict.size());
+    List<Map<scalar> > regionAngle(surfacesDict.size());
 
     label surfI = 0;
     forAllConstIter(dictionary, surfacesDict, iter)
@@ -274,6 +297,12 @@ Foam::refinementSurfaces::refinementSurfaces
             dict.lookup("zoneInside") >> zoneInside_[surfI];
         }
 
+        // Global perpendicular angle
+        if (dict.found("perpendicularAngle"))
+        {
+            globalAngle[surfI] = readScalar(dict.lookup("perpendicularAngle"));
+        }
+
         if (dict.found("regions"))
         {
             const dictionary& regionsDict = dict.subDict("regions");
@@ -306,6 +335,15 @@ Foam::refinementSurfaces::refinementSurfaces
 
                     regionMinLevel[surfI].insert(regionI, refLevel[0]);
                     regionMaxLevel[surfI].insert(regionI, refLevel[1]);
+
+                    if (regionDict.found("perpendicularAngle"))
+                    {
+                        regionAngle[surfI].insert
+                        (
+                            regionI,
+                            readScalar(regionDict.lookup("perpendicularAngle"))
+                        );
+                    }
                 }
             }
         }
@@ -326,6 +364,8 @@ Foam::refinementSurfaces::refinementSurfaces
     minLevel_ = 0;
     maxLevel_.setSize(nRegions);
     maxLevel_ = 0;
+    perpendicularAngle_.setSize(nRegions);
+    perpendicularAngle_ = -GREAT;
 
 
     forAll(globalMinLevel, surfI)
@@ -337,6 +377,7 @@ Foam::refinementSurfaces::refinementSurfaces
         {
             minLevel_[regionOffset_[surfI] + i] = globalMinLevel[surfI];
             maxLevel_[regionOffset_[surfI] + i] = globalMaxLevel[surfI];
+            perpendicularAngle_[regionOffset_[surfI] + i] = globalAngle[surfI];
         }
 
         // Overwrite with region specific information
@@ -346,6 +387,7 @@ Foam::refinementSurfaces::refinementSurfaces
 
             minLevel_[globalRegionI] = iter();
             maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
+            perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
 
             // Check validity
             if
@@ -454,8 +496,6 @@ void Foam::refinementSurfaces::setMinLevelFields
     const shellSurfaces& shells
 )
 {
-    //minLevelFields_.setSize(surfaces_.size());
-
     forAll(surfaces_, surfI)
     {
         const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
index a20df19f29..c727303173 100644
--- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -90,6 +90,9 @@ class refinementSurfaces
         //- From global region number to refinement level
         labelList maxLevel_;
 
+        //- From global region number to perpendicular angle
+        scalarField perpendicularAngle_;
+
 
     // Private Member Functions
 
@@ -178,6 +181,12 @@ public:
                 return maxLevel_;
             }
 
+            //- From global region number to perpendicular angle
+            const scalarField& perpendicularAngle() const
+            {
+                return perpendicularAngle_;
+            }
+
 
         // Helper