From 4b9f2ae6fcc531241ffef72781a99cbbf3d8516f Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 27 May 2013 11:15:39 +0100
Subject: [PATCH] ENH: snappyHexMesh: added proximity refinement

---
 .../autoHexMeshDriver/autoRefineDriver.C      | 404 +++++++++++++-
 .../autoHexMeshDriver/autoRefineDriver.H      |  17 +-
 .../meshRefinement/meshRefinement.H           |  44 ++
 .../meshRefinement/meshRefinementBaffles.C    |   4 +-
 .../meshRefinement/meshRefinementRefine.C     | 502 ++++++++++++++++++
 .../refinementSurfaces/refinementSurfaces.C   | 176 +++++-
 .../refinementSurfaces/refinementSurfaces.H   |  27 +-
 7 files changed, 1147 insertions(+), 27 deletions(-)

diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index 75ea461102..d78866bcd7 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -46,9 +46,6 @@ defineTypeNameAndDebug(autoRefineDriver, 0);
 } // End namespace Foam
 
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
@@ -97,12 +94,14 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
                 (
                     refineParams.keepPoints()[0],    // For now only use one.
                     refineParams.curvature(),
+                    refineParams.planarAngle(),
 
                     true,               // featureRefinement
                     false,              // featureDistanceRefinement
                     false,              // internalRefinement
                     false,              // surfaceRefinement
                     false,              // curvatureRefinement
+                    false,              // gapRefinement
                     refineParams.maxGlobalCells(),
                     refineParams.maxLocalCells()
                 )
@@ -208,12 +207,14 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
             (
                 refineParams.keepPoints()[0],
                 refineParams.curvature(),
+                refineParams.planarAngle(),
 
                 false,              // featureRefinement
                 false,              // featureDistanceRefinement
                 false,              // internalRefinement
                 true,               // surfaceRefinement
                 true,               // curvatureRefinement
+                false,              // gapRefinement
                 refineParams.maxGlobalCells(),
                 refineParams.maxLocalCells()
             )
@@ -294,6 +295,371 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
 }
 
 
+Foam::label Foam::autoRefineDriver::gapOnlyRefine
+(
+    const refinementParameters& refineParams,
+    const label maxIter
+)
+{
+    const fvMesh& mesh = meshRefiner_.mesh();
+
+    // Determine the maximum refinement level over all surfaces. This
+    // determines the minumum number of surface refinement iterations.
+
+    label maxIncrement = 0;
+    const labelList& maxLevel = meshRefiner_.surfaces().maxLevel();
+    const labelList& gapLevel = meshRefiner_.surfaces().gapLevel();
+
+    forAll(maxLevel, i)
+    {
+        maxIncrement = max(maxIncrement, gapLevel[i]-maxLevel[i]);
+    }
+
+    label iter = 0;
+
+    if (maxIncrement == 0)
+    {
+        return iter;
+    }
+
+    for (iter = 0; iter < maxIter; iter++)
+    {
+        Info<< nl
+            << "Gap refinement iteration " << iter << nl
+            << "--------------------------" << nl
+            << endl;
+
+
+        // Determine cells to refine
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~
+        // Only look at surface intersections (minLevel and surface curvature),
+        // do not do internal refinement (refinementShells)
+
+        labelList candidateCells
+        (
+            meshRefiner_.refineCandidates
+            (
+                refineParams.keepPoints()[0],
+                refineParams.curvature(),
+                refineParams.planarAngle(),
+
+                false,              // featureRefinement
+                false,              // featureDistanceRefinement
+                false,              // internalRefinement
+                false,              // surfaceRefinement
+                false,              // curvatureRefinement
+                true,               // gapRefinement
+                refineParams.maxGlobalCells(),
+                refineParams.maxLocalCells()
+            )
+        );
+
+        if (debug&meshRefinement::MESH)
+        {
+            Pout<< "Dumping " << candidateCells.size()
+                << " cells to cellSet candidateCellsFromGap." << endl;
+            cellSet c(mesh, "candidateCellsFromGap", candidateCells);
+            c.instance() = meshRefiner_.timeName();
+            c.write();
+        }
+
+        // Grow by one layer to make sure we're covering the gap
+        {
+            boolList isCandidateCell(mesh.nCells(), false);
+            forAll(candidateCells, i)
+            {
+                isCandidateCell[candidateCells[i]] = true;
+            }
+
+            for (label i=0; i<1; i++)
+            {
+                boolList newIsCandidateCell(isCandidateCell);
+
+                // Internal faces
+                for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+                {
+                    label own = mesh.faceOwner()[faceI];
+                    label nei = mesh.faceNeighbour()[faceI];
+
+                    if (isCandidateCell[own] != isCandidateCell[nei])
+                    {
+                        newIsCandidateCell[own] = true;
+                        newIsCandidateCell[nei] = true;
+                    }
+                }
+
+                // Get coupled boundary condition values
+                boolList neiIsCandidateCell;
+                syncTools::swapBoundaryCellList
+                (
+                    mesh,
+                    isCandidateCell,
+                    neiIsCandidateCell
+                );
+
+                // Boundary faces
+                for
+                (
+                    label faceI = mesh.nInternalFaces();
+                    faceI < mesh.nFaces();
+                    faceI++
+                )
+                {
+                    label own = mesh.faceOwner()[faceI];
+                    label bFaceI = faceI-mesh.nInternalFaces();
+
+                    if (isCandidateCell[own] != neiIsCandidateCell[bFaceI])
+                    {
+                        newIsCandidateCell[own] = true;
+                    }
+                }
+
+                isCandidateCell.transfer(newIsCandidateCell);
+            }
+
+            label n = 0;
+            forAll(isCandidateCell, cellI)
+            {
+                if (isCandidateCell[cellI])
+                {
+                    n++;
+                }
+            }
+            candidateCells.setSize(n);
+            n = 0;
+            forAll(isCandidateCell, cellI)
+            {
+                if (isCandidateCell[cellI])
+                {
+                    candidateCells[n++] = cellI;
+                }
+            }
+        }
+
+
+        if (debug&meshRefinement::MESH)
+        {
+            Pout<< "Dumping " << candidateCells.size()
+                << " cells to cellSet candidateCellsFromGapPlusBuffer." << endl;
+            cellSet c(mesh, "candidateCellsFromGapPlusBuffer", candidateCells);
+            c.instance() = meshRefiner_.timeName();
+            c.write();
+        }
+
+
+        labelList cellsToRefine
+        (
+            meshRefiner_.meshCutter().consistentRefinement
+            (
+                candidateCells,
+                true
+            )
+        );
+        Info<< "Determined cells to refine in = "
+            << mesh.time().cpuTimeIncrement() << " s" << endl;
+
+
+        label nCellsToRefine = cellsToRefine.size();
+        reduce(nCellsToRefine, sumOp<label>());
+
+        Info<< "Selected for refinement : " << nCellsToRefine
+            << " cells (out of " << mesh.globalData().nTotalCells()
+            << ')' << endl;
+
+        // Stop when no cells to refine or have done minimum necessary
+        // iterations and not enough cells to refine.
+        if
+        (
+            nCellsToRefine == 0
+         || (
+                iter >= maxIncrement
+             && nCellsToRefine <= refineParams.minRefineCells()
+            )
+        )
+        {
+            Info<< "Stopping refining since too few cells selected."
+                << nl << endl;
+            break;
+        }
+
+
+        if (debug)
+        {
+            const_cast<Time&>(mesh.time())++;
+        }
+
+
+        if
+        (
+            returnReduce
+            (
+                (mesh.nCells() >= refineParams.maxLocalCells()),
+                orOp<bool>()
+            )
+        )
+        {
+            meshRefiner_.balanceAndRefine
+            (
+                "gap refinement iteration " + name(iter),
+                decomposer_,
+                distributor_,
+                cellsToRefine,
+                refineParams.maxLoadUnbalance()
+            );
+        }
+        else
+        {
+            meshRefiner_.refineAndBalance
+            (
+                "gap refinement iteration " + name(iter),
+                decomposer_,
+                distributor_,
+                cellsToRefine,
+                refineParams.maxLoadUnbalance()
+            );
+        }
+    }
+    return iter;
+}
+
+
+Foam::label Foam::autoRefineDriver::danglingCellRefine
+(
+    const refinementParameters& refineParams,
+    const label nFaces,
+    const label maxIter
+)
+{
+    const fvMesh& mesh = meshRefiner_.mesh();
+
+    label iter;
+    for (iter = 0; iter < maxIter; iter++)
+    {
+        Info<< nl
+            << "Dangling coarse cells refinement iteration " << iter << nl
+            << "--------------------------------------------" << nl
+            << endl;
+
+
+        // Determine cells to refine
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        const cellList& cells = mesh.cells();
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+        labelList candidateCells;
+        {
+            cellSet candidateCellSet(mesh, "candidateCells", cells.size()/1000);
+
+            forAll(cells, cellI)
+            {
+                const cell& cFaces = cells[cellI];
+
+                label nIntFaces = 0;
+                forAll(cFaces, i)
+                {
+                    label bFaceI = cFaces[i]-mesh.nInternalFaces();
+                    if (bFaceI < 0)
+                    {
+                        nIntFaces++;
+                    }
+                    else
+                    {
+                        label patchI = pbm.patchID()[bFaceI];
+                        if (pbm[patchI].coupled())
+                        {
+                            nIntFaces++;
+                        }
+                    }
+                }
+
+                if (nIntFaces == nFaces)
+                {
+                    candidateCellSet.insert(cellI);
+                }
+            }
+
+            if (debug&meshRefinement::MESH)
+            {
+                Pout<< "Dumping " << candidateCellSet.size()
+                    << " cells to cellSet candidateCellSet." << endl;
+                candidateCellSet.instance() = meshRefiner_.timeName();
+                candidateCellSet.write();
+            }
+            candidateCells = candidateCellSet.toc();
+        }
+
+
+
+        labelList cellsToRefine
+        (
+            meshRefiner_.meshCutter().consistentRefinement
+            (
+                candidateCells,
+                true
+            )
+        );
+        Info<< "Determined cells to refine in = "
+            << mesh.time().cpuTimeIncrement() << " s" << endl;
+
+
+        label nCellsToRefine = cellsToRefine.size();
+        reduce(nCellsToRefine, sumOp<label>());
+
+        Info<< "Selected for refinement : " << nCellsToRefine
+            << " cells (out of " << mesh.globalData().nTotalCells()
+            << ')' << endl;
+
+        // Stop when no cells to refine. No checking of minRefineCells since
+        // too few cells
+        if (nCellsToRefine == 0)
+        {
+            Info<< "Stopping refining since too few cells selected."
+                << nl << endl;
+            break;
+        }
+
+
+        if (debug)
+        {
+            const_cast<Time&>(mesh.time())++;
+        }
+
+
+        if
+        (
+            returnReduce
+            (
+                (mesh.nCells() >= refineParams.maxLocalCells()),
+                orOp<bool>()
+            )
+        )
+        {
+            meshRefiner_.balanceAndRefine
+            (
+                "coarse cell refinement iteration " + name(iter),
+                decomposer_,
+                distributor_,
+                cellsToRefine,
+                refineParams.maxLoadUnbalance()
+            );
+        }
+        else
+        {
+            meshRefiner_.refineAndBalance
+            (
+                "coarse cell refinement iteration " + name(iter),
+                decomposer_,
+                distributor_,
+                cellsToRefine,
+                refineParams.maxLoadUnbalance()
+            );
+        }
+    }
+    return iter;
+}
+
+
 void Foam::autoRefineDriver::removeInsideCells
 (
     const refinementParameters& refineParams,
@@ -371,12 +737,14 @@ Foam::label Foam::autoRefineDriver::shellRefine
             (
                 refineParams.keepPoints()[0],
                 refineParams.curvature(),
+                refineParams.planarAngle(),
 
                 false,              // featureRefinement
                 true,               // featureDistanceRefinement
                 true,               // internalRefinement
                 false,              // surfaceRefinement
                 false,              // curvatureRefinement
+                false,              // gapRefinement
                 refineParams.maxGlobalCells(),
                 refineParams.maxLocalCells()
             )
@@ -653,6 +1021,16 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
             // Debug:test all is still synced across proc patches
             meshRefiner_.checkData();
         }
+
+        // Remove any now dangling parts
+        meshRefiner_.splitMeshRegions(refineParams.keepPoints()[0]);
+
+        if (debug)
+        {
+            // Debug:test all is still synced across proc patches
+            meshRefiner_.checkData();
+        }
+
         Info<< "Merged free-standing baffles in = "
             << mesh.time().cpuTimeIncrement() << " s." << endl;
     }
@@ -732,6 +1110,12 @@ void Foam::autoRefineDriver::doRefine
         100     // maxIter
     );
 
+    gapOnlyRefine
+    (
+        refineParams,
+        100     // maxIter
+    );
+
     // Remove cells (a certain distance) beyond surface intersections
     removeInsideCells
     (
@@ -746,6 +1130,20 @@ void Foam::autoRefineDriver::doRefine
         100    // maxIter
     );
 
+    // Refine any hexes with 5 or 6 faces refined to make smooth edges
+    danglingCellRefine
+    (
+        refineParams,
+        21,     // 1 coarse face + 5 refined faces
+        100     // maxIter
+    );
+    danglingCellRefine
+    (
+        refineParams,
+        24,     // 0 coarse faces + 6 refined faces
+        100     // maxIter
+    );
+
     // Introduce baffles at surface intersections. Remove sections unreachable
     // from keepPoint.
     baffleAndSplitMesh(refineParams, prepareForSnapping, motionDict);
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
index 37a0f6b102..c44a14a794 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.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
@@ -88,6 +88,21 @@ class autoRefineDriver
             const label maxIter
         );
 
+        //- Refine all cells in small gaps
+        label gapOnlyRefine
+        (
+            const refinementParameters& refineParams,
+            const label maxIter
+        );
+
+        //- Refine cells with almost all sides refined
+        label danglingCellRefine
+        (
+            const refinementParameters& refineParams,
+            const label nFaces,
+            const label maxIter
+        );
+
         //- Remove all cells within intersected region
         void removeInsideCells
         (
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 60e6261ccc..7487e8a2bf 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -335,6 +335,48 @@ private:
                 label& nRefine
             ) const;
 
+            //- Is local topology a small gap?
+            bool isGap
+            (
+                const scalar,
+                const vector&,
+                const vector&,
+                const vector&,
+                const vector&
+            ) const;
+
+            //- Mark cell if local a gap topology or
+            bool checkProximity
+            (
+                const scalar planarCos,
+                const label nAllowRefine,
+
+                const label surfaceLevel,
+                const vector& surfaceLocation,
+                const vector& surfaceNormal,
+
+                const label cellI,
+
+                label& cellMaxLevel,
+                vector& cellMaxLocation,
+                vector& cellMaxNormal,
+
+                labelList& refineCell,
+                label& nRefine
+            ) const;
+
+            //- Mark cells for surface proximity based refinement.
+            label markProximityRefinement
+            (
+                const scalar curvature,
+                const label nAllowRefine,
+                const labelList& neiLevel,
+                const pointField& neiCc,
+
+                labelList& refineCell,
+                label& nRefine
+            ) const;
+
         // Baffle handling
 
             //- Get faces to repatch. Returns map from face to patch.
@@ -676,12 +718,14 @@ public:
             (
                 const point& keepPoint,
                 const scalar curvature,
+                const scalar planarAngle,
 
                 const bool featureRefinement,
                 const bool featureDistanceRefinement,
                 const bool internalRefinement,
                 const bool surfaceRefinement,
                 const bool curvatureRefinement,
+                const bool gapRefinement,
                 const label maxGlobalCells,
                 const label maxLocalCells
             ) const;
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 7a4f87753d..628dd7898a 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -1965,7 +1965,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
     const bool removeEdgeConnectedCells,
     const scalarField& perpendicularAngle,
     const bool mergeFreeStanding,
-    const scalar freeStandingAngle,
+    const scalar planarAngle,
     const dictionary& motionDict,
     Time& runTime,
     const labelList& globalToMasterPatch,
@@ -2172,7 +2172,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
                     identity(mesh_.nFaces()-mesh_.nInternalFaces())
                    +mesh_.nInternalFaces()
                 ),
-                freeStandingAngle
+                planarAngle
             )
         );
 
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index ee734bcebc..9667750cf5 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -38,6 +38,7 @@ License
 #include "featureEdgeMesh.H"
 #include "Cloud.H"
 //#include "globalIndex.H"
+//#include "OBJstream.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -1022,6 +1023,9 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
             end,
             minLevel,           // max level of surface has to be bigger
                                 // than min level of neighbouring cells
+
+            surfaces_.maxLevel(),
+
             surfaceNormal,
             surfaceLevel
         );
@@ -1210,6 +1214,473 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
 }
 
 
+// Mark small gaps
+bool Foam::meshRefinement::isGap
+(
+    const scalar planarCos,
+    const vector& point0,
+    const vector& normal0,
+
+    const vector& point1,
+    const vector& normal1
+) const
+{
+    ////- hits differ and angles are oppositeish
+    //return
+    //    (mag(point0-point1) > mergeDistance())
+    // && ((normal0 & normal1) < (-1+planarCos));
+
+
+    //- hits differ and angles are oppositeish and
+    //  hits have a normal distance
+    vector d = point1-point0;
+    scalar magD = mag(d);
+
+    if (magD > mergeDistance())
+    {
+        scalar cosAngle = (normal0 & normal1);
+
+        vector avg = vector::zero;
+        if (cosAngle < (-1+planarCos))
+        {
+            // Opposite normals
+            avg = 0.5*(normal0-normal1);
+        }
+        else if (cosAngle > (1-planarCos))
+        {
+            avg = 0.5*(normal0+normal1);
+        }
+
+        if (avg != vector::zero)
+        {
+            avg /= mag(avg);
+            d /= magD;
+
+            // Check average normal with respect to intersection locations
+            if (mag(avg&d) > Foam::cos(degToRad(45)))
+            {
+                return true;
+            }
+            else
+            {
+                return false;
+            }
+        }
+        else
+        {
+            return false;
+        }
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+bool Foam::meshRefinement::checkProximity
+(
+    const scalar planarCos,
+    const label nAllowRefine,
+
+    const label surfaceLevel,       // current intersection max level
+    const vector& surfaceLocation,  // current intersection location
+    const vector& surfaceNormal,    // current intersection normal
+
+    const label cellI,
+
+    label& cellMaxLevel,        // cached max surface level for this cell
+    vector& cellMaxLocation,    // cached surface normal for this cell
+    vector& cellMaxNormal,      // cached surface normal for this cell
+
+    labelList& refineCell,
+    label& nRefine
+) const
+{
+    const labelList& cellLevel = meshCutter_.cellLevel();
+
+    // Test if surface applicable
+    if (surfaceLevel > cellLevel[cellI])
+    {
+        if (cellMaxLevel == -1)
+        {
+            // First visit of cell. Store
+            cellMaxLevel = surfaceLevel;
+            cellMaxLocation = surfaceLocation;
+            cellMaxNormal = surfaceNormal;
+        }
+        else
+        {
+            // Second or more visit.
+            // Check if
+            //  - different location
+            //  - opposite surface
+
+            bool closeSurfaces = isGap
+            (
+                planarCos,
+                cellMaxLocation,
+                cellMaxNormal,
+                surfaceLocation,
+                surfaceNormal
+            );
+
+            // Set normal to that of highest surface. Not really necessary
+            // over here but we reuse cellMax info when doing coupled faces.
+            if (surfaceLevel > cellMaxLevel)
+            {
+                cellMaxLevel = surfaceLevel;
+                cellMaxLocation = surfaceLocation;
+                cellMaxNormal = surfaceNormal;
+            }
+
+
+            if (closeSurfaces)
+            {
+                //Pout<< "Found gap:" << nl
+                //    << "    location:" << surfaceLocation
+                //    << "\tnormal  :" << surfaceNormal << nl
+                ///    << "    location:" << cellMaxLocation
+                //    << "\tnormal  :" << cellMaxNormal << nl
+                //    << "\tcos     :" << (surfaceNormal&cellMaxNormal) << nl
+                //    << endl;
+
+                return markForRefine
+                (
+                    surfaceLevel,   // mark with any non-neg number.
+                    nAllowRefine,
+                    refineCell[cellI],
+                    nRefine
+                );
+            }
+        }
+    }
+
+    // Did not reach refinement limit.
+    return true;
+}
+
+
+Foam::label Foam::meshRefinement::markProximityRefinement
+(
+    const scalar planarCos,
+    const label nAllowRefine,
+    const labelList& neiLevel,
+    const pointField& neiCc,
+
+    labelList& refineCell,
+    label& nRefine
+) const
+{
+    const labelList& cellLevel = meshCutter_.cellLevel();
+    const pointField& cellCentres = mesh_.cellCentres();
+
+    label oldNRefine = nRefine;
+
+    // 1. local test: any cell on more than one surface gets refined
+    // (if its current level is < max of the surface max level)
+
+    // 2. 'global' test: any cell on only one surface with a neighbour
+    // on a different surface gets refined (if its current level etc.)
+
+
+    // Collect candidate faces (i.e. intersecting any surface and
+    // owner/neighbour not yet refined.
+    labelList testFaces(getRefineCandidateFaces(refineCell));
+
+    // Collect segments
+    pointField start(testFaces.size());
+    pointField end(testFaces.size());
+    labelList minLevel(testFaces.size());
+
+    forAll(testFaces, i)
+    {
+        label faceI = testFaces[i];
+
+        label own = mesh_.faceOwner()[faceI];
+
+        if (mesh_.isInternalFace(faceI))
+        {
+            label nei = mesh_.faceNeighbour()[faceI];
+
+            start[i] = cellCentres[own];
+            end[i] = cellCentres[nei];
+            minLevel[i] = min(cellLevel[own], cellLevel[nei]);
+        }
+        else
+        {
+            label bFaceI = faceI - mesh_.nInternalFaces();
+
+            start[i] = cellCentres[own];
+            end[i] = neiCc[bFaceI];
+            minLevel[i] = min(cellLevel[own], neiLevel[bFaceI]);
+        }
+    }
+
+    // Extend segments a bit
+    {
+        const vectorField smallVec(Foam::sqrt(SMALL)*(end-start));
+        start -= smallVec;
+        end += smallVec;
+    }
+
+
+    // Test for all intersections (with surfaces of higher gap level than
+    // minLevel) and cache per cell the max surface level and the local normal
+    // on that surface.
+    labelList cellMaxLevel(mesh_.nCells(), -1);
+    vectorField cellMaxNormal(mesh_.nCells(), vector::zero);
+    pointField cellMaxLocation(mesh_.nCells(), vector::zero);
+
+    {
+        // Per segment the normals of the surfaces
+        List<vectorList> surfaceLocation;
+        List<vectorList> surfaceNormal;
+        // Per segment the list of levels of the surfaces
+        labelListList surfaceLevel;
+
+        surfaces_.findAllHigherIntersections
+        (
+            start,
+            end,
+            minLevel,           // gap level of surface has to be bigger
+                                // than min level of neighbouring cells
+
+            surfaces_.gapLevel(),
+
+            surfaceLocation,
+            surfaceNormal,
+            surfaceLevel
+        );
+        // Clear out unnecessary data
+        start.clear();
+        end.clear();
+        minLevel.clear();
+
+        //// Extract per cell information on the surface with the highest max
+        //OBJstream str
+        //(
+        //    mesh_.time().path()
+        //  / "findAllHigherIntersections_"
+        //  + mesh_.time().timeName()
+        //  + ".obj"
+        //);
+        //// All intersections
+        //OBJstream str2
+        //(
+        //    mesh_.time().path()
+        //  / "findAllHigherIntersections2_"
+        //  + mesh_.time().timeName()
+        //  + ".obj"
+        //);
+
+        forAll(testFaces, i)
+        {
+            label faceI = testFaces[i];
+            label own = mesh_.faceOwner()[faceI];
+
+            const labelList& fLevels = surfaceLevel[i];
+            const vectorList& fPoints = surfaceLocation[i];
+            const vectorList& fNormals = surfaceNormal[i];
+
+            forAll(fLevels, hitI)
+            {
+                checkProximity
+                (
+                    planarCos,
+                    nAllowRefine,
+
+                    fLevels[hitI],
+                    fPoints[hitI],
+                    fNormals[hitI],
+
+                    own,
+                    cellMaxLevel[own],
+                    cellMaxLocation[own],
+                    cellMaxNormal[own],
+
+                    refineCell,
+                    nRefine
+                );
+            }
+
+            if (mesh_.isInternalFace(faceI))
+            {
+                label nei = mesh_.faceNeighbour()[faceI];
+
+                forAll(fLevels, hitI)
+                {
+                    checkProximity
+                    (
+                        planarCos,
+                        nAllowRefine,
+
+                        fLevels[hitI],
+                        fPoints[hitI],
+                        fNormals[hitI],
+
+                        nei,
+                        cellMaxLevel[nei],
+                        cellMaxLocation[nei],
+                        cellMaxNormal[nei],
+
+                        refineCell,
+                        nRefine
+                    );
+                }
+            }
+        }
+    }
+
+    // 2. Find out a measure of surface curvature and region edges.
+    // Send over surface region and surface normal to neighbour cell.
+
+    labelList neiBndMaxLevel(mesh_.nFaces()-mesh_.nInternalFaces());
+    pointField neiBndMaxLocation(mesh_.nFaces()-mesh_.nInternalFaces());
+    vectorField neiBndMaxNormal(mesh_.nFaces()-mesh_.nInternalFaces());
+
+    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
+    {
+        label bFaceI = faceI-mesh_.nInternalFaces();
+        label own = mesh_.faceOwner()[faceI];
+
+        neiBndMaxLevel[bFaceI] = cellMaxLevel[own];
+        neiBndMaxLocation[bFaceI] = cellMaxLocation[own];
+        neiBndMaxNormal[bFaceI] = cellMaxNormal[own];
+    }
+    syncTools::swapBoundaryFaceList(mesh_, neiBndMaxLevel);
+    syncTools::swapBoundaryFaceList(mesh_, neiBndMaxLocation);
+    syncTools::swapBoundaryFaceList(mesh_, neiBndMaxNormal);
+
+    // Loop over all faces. Could only be checkFaces.. except if they're coupled
+
+    // Internal faces
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        label own = mesh_.faceOwner()[faceI];
+        label nei = mesh_.faceNeighbour()[faceI];
+
+        if (cellMaxLevel[own] != -1 && cellMaxLevel[nei] != -1)
+        {
+            // Have valid data on both sides. Check planarCos.
+            if
+            (
+                isGap
+                (
+                    planarCos,
+                    cellMaxLocation[own],
+                    cellMaxNormal[own],
+                    cellMaxLocation[nei],
+                    cellMaxNormal[nei]
+                )
+            )
+            {
+                // See which side to refine
+                if (cellLevel[own] < cellMaxLevel[own])
+                {
+                    if
+                    (
+                        !markForRefine
+                        (
+                            cellMaxLevel[own],
+                            nAllowRefine,
+                            refineCell[own],
+                            nRefine
+                        )
+                    )
+                    {
+                        if (debug)
+                        {
+                            Pout<< "Stopped refining since reaching my cell"
+                                << " limit of " << mesh_.nCells()+7*nRefine
+                                << endl;
+                        }
+                        break;
+                    }
+                }
+
+                if (cellLevel[nei] < cellMaxLevel[nei])
+                {
+                    if
+                    (
+                        !markForRefine
+                        (
+                            cellMaxLevel[nei],
+                            nAllowRefine,
+                            refineCell[nei],
+                            nRefine
+                        )
+                    )
+                    {
+                        if (debug)
+                        {
+                            Pout<< "Stopped refining since reaching my cell"
+                                << " limit of " << mesh_.nCells()+7*nRefine
+                                << endl;
+                        }
+                        break;
+                    }
+                }
+            }
+        }
+    }
+    // Boundary faces
+    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
+    {
+        label own = mesh_.faceOwner()[faceI];
+        label bFaceI = faceI - mesh_.nInternalFaces();
+
+        if (cellLevel[own] < cellMaxLevel[own] && neiBndMaxLevel[bFaceI] != -1)
+        {
+            // Have valid data on both sides. Check planarCos.
+            if
+            (
+                isGap
+                (
+                    planarCos,
+                    cellMaxLocation[own],
+                    cellMaxNormal[own],
+                    neiBndMaxLocation[bFaceI],
+                    neiBndMaxNormal[bFaceI]
+                )
+            )
+            {
+                if
+                (
+                    !markForRefine
+                    (
+                        cellMaxLevel[own],
+                        nAllowRefine,
+                        refineCell[own],
+                        nRefine
+                    )
+                )
+                {
+                    if (debug)
+                    {
+                        Pout<< "Stopped refining since reaching my cell"
+                            << " limit of " << mesh_.nCells()+7*nRefine
+                            << endl;
+                    }
+                    break;
+                }
+            }
+        }
+    }
+
+    if
+    (
+        returnReduce(nRefine, sumOp<label>())
+      > returnReduce(nAllowRefine, sumOp<label>())
+    )
+    {
+        Info<< "Reached refinement limit." << endl;
+    }
+
+    return returnReduce(nRefine-oldNRefine, sumOp<label>());
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 // Calculate list of cells to refine. Gets for any edge (start - end)
@@ -1221,12 +1692,14 @@ Foam::labelList Foam::meshRefinement::refineCandidates
 (
     const point& keepPoint,
     const scalar curvature,
+    const scalar planarAngle,
 
     const bool featureRefinement,
     const bool featureDistanceRefinement,
     const bool internalRefinement,
     const bool surfaceRefinement,
     const bool curvatureRefinement,
+    const bool gapRefinement,
     const label maxGlobalCells,
     const label maxLocalCells
 ) const
@@ -1370,6 +1843,35 @@ Foam::labelList Foam::meshRefinement::refineCandidates
                 << ": " << nCurv << " cells."  << endl;
         }
 
+
+        const scalar planarCos = Foam::cos(degToRad(planarAngle));
+
+        if
+        (
+            gapRefinement
+         && (planarCos >= -1 && planarCos <= 1)
+         && (max(surfaces_.gapLevel()) > -1)
+        )
+        {
+            Info<< "Specified gap level : " << max(surfaces_.gapLevel())
+                << ", planar angle " << planarAngle << endl;
+
+            label nGap = markProximityRefinement
+            (
+                planarCos,
+                nAllowRefine,
+                neiLevel,
+                neiCc,
+
+                refineCell,
+                nRefine
+            );
+            Info<< "Marked for refinement due to close opposite surfaces       "
+                << ": " << nGap << " cells."  << endl;
+        }
+
+
+
         // Pack cells-to-refine
         // ~~~~~~~~~~~~~~~~~~~~
 
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index f53968532c..f1917ee6f9 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -116,10 +116,12 @@ Foam::refinementSurfaces::refinementSurfaces
 
     labelList globalMinLevel(surfI, 0);
     labelList globalMaxLevel(surfI, 0);
+    labelList globalLevelIncr(surfI, 0);
     scalarField globalAngle(surfI, -GREAT);
     PtrList<dictionary> globalPatchInfo(surfI);
     List<Map<label> > regionMinLevel(surfI);
     List<Map<label> > regionMaxLevel(surfI);
+    List<Map<label> > regionLevelIncr(surfI);
     List<Map<scalar> > regionAngle(surfI);
     List<Map<autoPtr<dictionary> > > regionPatchInfo(surfI);
 
@@ -138,6 +140,31 @@ Foam::refinementSurfaces::refinementSurfaces
             const labelPair refLevel(dict.lookup("level"));
             globalMinLevel[surfI] = refLevel[0];
             globalMaxLevel[surfI] = refLevel[1];
+            globalLevelIncr[surfI] = dict.lookupOrDefault
+            (
+                "gapLevelIncrement",
+                0
+            );
+
+            if
+            (
+                globalMinLevel[surfI] < 0
+             || globalMaxLevel[surfI] < globalMinLevel[surfI]
+             || globalLevelIncr[surfI] < 0
+            )
+            {
+                FatalErrorIn
+                (
+                    "refinementSurfaces::refinementSurfaces"
+                    "(const searchableSurfaces&, const dictionary>&"
+                )   << "Illegal level specification for surface "
+                    << names_[surfI]
+                    << " : minLevel:" << globalMinLevel[surfI]
+                    << " maxLevel:" << globalMaxLevel[surfI]
+                    << " levelIncrement:" << globalLevelIncr[surfI]
+                    << exit(FatalError);
+            }
+
 
             // Global zone names per surface
             if (dict.readIfPresent("faceZone", faceZoneNames_[surfI]))
@@ -244,6 +271,34 @@ Foam::refinementSurfaces::refinementSurfaces
 
                         regionMinLevel[surfI].insert(regionI, refLevel[0]);
                         regionMaxLevel[surfI].insert(regionI, refLevel[1]);
+                        label levelIncr = regionDict.lookupOrDefault
+                        (
+                            "gapLevelIncrement",
+                            0
+                        );
+                        regionLevelIncr[surfI].insert(regionI, levelIncr);
+
+                        if
+                        (
+                            refLevel[0] < 0
+                         || refLevel[1] < refLevel[0]
+                         || levelIncr < 0
+                        )
+                        {
+                            FatalErrorIn
+                            (
+                                "refinementSurfaces::refinementSurfaces"
+                                "(const searchableSurfaces&, const dictionary>&"
+                            )   << "Illegal level specification for surface "
+                                << names_[surfI] << " region "
+                                << regionNames[regionI]
+                                << " : minLevel:" << refLevel[0]
+                                << " maxLevel:" << refLevel[1]
+                                << " levelIncrement:" << levelIncr
+                                << exit(FatalError);
+                        }
+
+
 
                         if (regionDict.found("perpendicularAngle"))
                         {
@@ -286,6 +341,8 @@ Foam::refinementSurfaces::refinementSurfaces
     minLevel_ = 0;
     maxLevel_.setSize(nRegions);
     maxLevel_ = 0;
+    gapLevel_.setSize(nRegions);
+    gapLevel_ = -1;
     perpendicularAngle_.setSize(nRegions);
     perpendicularAngle_ = -GREAT;
     patchInfo_.setSize(nRegions);
@@ -301,6 +358,10 @@ Foam::refinementSurfaces::refinementSurfaces
             label globalRegionI = regionOffset_[surfI] + i;
             minLevel_[globalRegionI] = globalMinLevel[surfI];
             maxLevel_[globalRegionI] = globalMaxLevel[surfI];
+            gapLevel_[globalRegionI] =
+                maxLevel_[globalRegionI]
+              + globalLevelIncr[surfI];
+
             perpendicularAngle_[globalRegionI] = globalAngle[surfI];
             if (globalPatchInfo.set(surfI))
             {
@@ -319,24 +380,9 @@ Foam::refinementSurfaces::refinementSurfaces
 
             minLevel_[globalRegionI] = iter();
             maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
-
-            // Check validity
-            if
-            (
-                minLevel_[globalRegionI] < 0
-             || maxLevel_[globalRegionI] < minLevel_[globalRegionI]
-            )
-            {
-                FatalErrorIn
-                (
-                    "refinementSurfaces::refinementSurfaces"
-                    "(const searchableSurfaces&, const dictionary>&"
-                )   << "Illegal level or layer specification for surface "
-                    << names_[surfI]
-                    << " : minLevel:" << minLevel_[globalRegionI]
-                    << " maxLevel:" << maxLevel_[globalRegionI]
-                    << exit(FatalError);
-            }
+            gapLevel_[globalRegionI] =
+                maxLevel_[globalRegionI]
+              + regionLevelIncr[surfI][iter.key()];
         }
         forAllConstIter(Map<scalar>, regionAngle[surfI], iter)
         {
@@ -714,6 +760,8 @@ void Foam::refinementSurfaces::findAllHigherIntersections
     const pointField& end,
     const labelList& currentLevel,   // current cell refinement level
 
+    const labelList& globalRegionLevel,
+
     List<vectorList>& surfaceNormal,
     labelListList& surfaceLevel
 ) const
@@ -777,7 +825,7 @@ void Foam::refinementSurfaces::findAllHigherIntersections
             label region = globalRegion(surfI, surfRegion[i]);
             label pointI = pointMap[i];
 
-            if (maxLevel_[region] > currentLevel[pointI])
+            if (globalRegionLevel[region] > currentLevel[pointI])
             {
                 // Append to pointI info
                 label sz = surfaceNormal[pointI].size();
@@ -785,7 +833,95 @@ void Foam::refinementSurfaces::findAllHigherIntersections
                 surfaceNormal[pointI][sz] = surfNormal[i];
 
                 surfaceLevel[pointI].setSize(sz+1);
-                surfaceLevel[pointI][sz] = maxLevel_[region];
+                surfaceLevel[pointI][sz] = globalRegionLevel[region];
+            }
+        }
+    }
+}
+
+
+void Foam::refinementSurfaces::findAllHigherIntersections
+(
+    const pointField& start,
+    const pointField& end,
+    const labelList& currentLevel,   // current cell refinement level
+
+    const labelList& globalRegionLevel,
+
+    List<pointList>& surfaceLocation,
+    List<vectorList>& surfaceNormal,
+    labelListList& surfaceLevel
+) const
+{
+    surfaceLevel.setSize(start.size());
+    surfaceNormal.setSize(start.size());
+    surfaceLocation.setSize(start.size());
+
+    if (surfaces_.empty())
+    {
+        return;
+    }
+
+    // Work arrays
+    List<List<pointIndexHit> > hitInfo;
+    labelList pRegions;
+    vectorField pNormals;
+
+    forAll(surfaces_, surfI)
+    {
+        allGeometry_[surfaces_[surfI]].findLineAll(start, end, hitInfo);
+
+        // Repack hits for surface into flat list
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // To avoid overhead of calling getRegion for every point
+
+        label n = 0;
+        forAll(hitInfo, pointI)
+        {
+            n += hitInfo[pointI].size();
+        }
+
+        List<pointIndexHit> surfInfo(n);
+        labelList pointMap(n);
+        n = 0;
+
+        forAll(hitInfo, pointI)
+        {
+            const List<pointIndexHit>& pHits = hitInfo[pointI];
+
+            forAll(pHits, i)
+            {
+                surfInfo[n] = pHits[i];
+                pointMap[n] = pointI;
+                n++;
+            }
+        }
+
+        labelList surfRegion(n);
+        vectorField surfNormal(n);
+        allGeometry_[surfaces_[surfI]].getRegion(surfInfo, surfRegion);
+        allGeometry_[surfaces_[surfI]].getNormal(surfInfo, surfNormal);
+
+        // Extract back into pointwise
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        forAll(surfRegion, i)
+        {
+            label region = globalRegion(surfI, surfRegion[i]);
+            label pointI = pointMap[i];
+
+            if (globalRegionLevel[region] > currentLevel[pointI])
+            {
+                // Append to pointI info
+                label sz = surfaceNormal[pointI].size();
+                surfaceLocation[pointI].setSize(sz+1);
+                surfaceLocation[pointI][sz] = surfInfo[i].hitPoint();
+
+                surfaceNormal[pointI].setSize(sz+1);
+                surfaceNormal[pointI][sz] = surfNormal[i];
+
+                surfaceLevel[pointI].setSize(sz+1);
+                surfaceLevel[pointI][sz] = globalRegionLevel[region];
             }
         }
     }
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
index aed21821b5..c9773398e5 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.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
@@ -51,6 +51,8 @@ class searchableSurfaces;
 class shellSurfaces;
 class triSurfaceMesh;
 
+typedef List<point> pointList;
+
 /*---------------------------------------------------------------------------*\
                            Class refinementSurfaces Declaration
 \*---------------------------------------------------------------------------*/
@@ -120,6 +122,9 @@ private:
         //- From global region number to refinement level
         labelList maxLevel_;
 
+        //- From global region number to small-gap level
+        labelList gapLevel_;
+
         //- From global region number to perpendicular angle
         scalarField perpendicularAngle_;
 
@@ -226,6 +231,12 @@ public:
                 return maxLevel_;
             }
 
+            //- From global region number to small gap refinement level
+            const labelList& gapLevel() const
+            {
+                return gapLevel_;
+            }
+
             //- From global region number to perpendicular angle
             const scalarField& perpendicularAngle() const
             {
@@ -295,11 +306,25 @@ public:
                 const pointField& start,
                 const pointField& end,
                 const labelList& currentLevel,  // current cell refinement level
+                const labelList& globalRegionLevel, // level per surfregion
 
                 List<vectorList>& surfaceNormal,
                 labelListList& surfaceLevel
             ) const;
 
+            //- Find all intersections of edge. Unsorted order.
+            void findAllHigherIntersections
+            (
+                const pointField& start,
+                const pointField& end,
+                const labelList& currentLevel,  // current cell refinement level
+                const labelList& globalRegionLevel, // level per surfregion
+
+                List<pointList>& surfaceLocation,
+                List<vectorList>& surfaceNormal,
+                labelListList& surfaceLevel
+            ) const;
+
             //- Find intersection nearest to the endpoints. surface1,2 are
             //  not indices into surfacesToTest but refinement surface indices.
             //  Returns surface, region on surface (so not global surface)