ENH: snappyHexMesh: add automatic gap-level detection and refinement

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C

@@ -1075,6 +1075,26 @@ int main(int argc, char *argv[])
         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
 
 
+    // Optionally read limit shells
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    const dictionary limitDict(refineDict.subOrEmptyDict("limitRegions"));
+
+    if (!limitDict.empty())
+    {
+        Info<< "Reading limit shells." << endl;
+    }
+
+    shellSurfaces limitShells(allGeometry, limitDict);
+
+    if (!limitDict.empty())
+    {
+        Info<< "Read refinement shells in = "
+            << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
+    }
+
+
+
     // Read feature meshes
     // ~~~~~~~~~~~~~~~~~~~
 
@@ -1105,7 +1125,8 @@ int main(int argc, char *argv[])
         overwrite,          // overwrite mesh files?
         surfaces,           // for surface intersection refinement
         features,           // for feature edges/point based refinement
-        shells              // for volume (inside/outside) refinement
+        shells,             // for volume (inside/outside) refinement
+        limitShells         // limit of volume refinement
     );
     Info<< "Calculated surface intersections in = "
         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -95,10 +95,10 @@ castellatedMeshControls
     // actually be a lot less.
     maxGlobalCells 2000000;
 
-    // The surface refinement loop might spend lots of iterations refining just a
+    // The surface refinement loop might spend lots of iterations refining just
-    // few cells. This setting will cause refinement to stop if <= minimumRefine
+    // a few cells. This setting will cause refinement to stop if
-    // are selected for refinement. Note: it will at least do one iteration
+    // <= minimumRefine cells are selected for refinement. Note: it will
-    // (unless the number of cells to refine is 0)
+    // at least do one iteration (unless the number of cells to refine is 0)
     minRefinementCells 0;
 
     // Allow a certain level of imbalance during refining
@@ -249,6 +249,13 @@ castellatedMeshControls
         //}
     }
 
+
+    // Limit refinement in geometric region
+    limitRegions
+    {
+    }
+
+
     // Mesh selection
     // ~~~~~~~~~~~~~~
 

src/OpenFOAM/algorithms/indexedOctree/volumeType.H

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -95,6 +95,12 @@ public:
             t_(t)
         {}
 
+        //- Construct from integer
+        explicit volumeType(const int t)
+        :
+            t_(static_cast<volumeType::type>(t))
+        {}
+
 
     // Member Functions
 

src/mesh/autoMesh/Make/files

@@ -15,6 +15,7 @@ $(autoHexMesh)/meshRefinement/meshRefinement.C
 $(autoHexMesh)/meshRefinement/meshRefinementMerge.C
 $(autoHexMesh)/meshRefinement/meshRefinementProblemCells.C
 $(autoHexMesh)/meshRefinement/meshRefinementRefine.C
+$(autoHexMesh)/meshRefinement/meshRefinementGapRefine.C
 $(autoHexMesh)/meshRefinement/patchFaceOrientation.C
 
 $(autoHexMesh)/refinementFeatures/refinementFeatures.C

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C

@@ -105,7 +105,9 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
                     false,              // internalRefinement
                     false,              // surfaceRefinement
                     false,              // curvatureRefinement
+                    false,              // smallFeatureRefinement
                     false,              // gapRefinement
+                    false,              // bigGapRefinement
                     refineParams.maxGlobalCells(),
                     refineParams.maxLocalCells()
                 )
@@ -179,6 +181,127 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
 }
 
 
+Foam::label Foam::autoRefineDriver::smallFeatureRefine
+(
+    const refinementParameters& refineParams,
+    const label maxIter
+)
+{
+    const fvMesh& mesh = meshRefiner_.mesh();
+
+
+    label iter = 0;
+
+    // See if any surface has an extendedGapLevel
+    labelList surfaceMaxLevel(meshRefiner_.surfaces().maxGapLevel());
+    labelList shellMaxLevel(meshRefiner_.shells().maxGapLevel());
+
+    if (max(surfaceMaxLevel) == 0 && max(shellMaxLevel) == 0)
+    {
+        return iter;
+    }
+
+    for (; iter < maxIter; iter++)
+    {
+        Info<< nl
+            << "Small surface feature refinement iteration " << iter << nl
+            << "--------------------------------------------" << nl
+            << endl;
+
+
+        // Determine cells to refine
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        labelList candidateCells
+        (
+            meshRefiner_.refineCandidates
+            (
+                refineParams.locationsInMesh(),
+                refineParams.curvature(),
+                refineParams.planarAngle(),
+
+                false,              // featureRefinement
+                false,              // featureDistanceRefinement
+                false,              // internalRefinement
+                false,              // surfaceRefinement
+                false,              // curvatureRefinement
+                true,               // smallFeatureRefinement
+                false,              // gapRefinement
+                false,              // bigGapRefinement
+                refineParams.maxGlobalCells(),
+                refineParams.maxLocalCells()
+            )
+        );
+
+        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)
+        {
+            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
+            (
+                "small feature refinement iteration " + name(iter),
+                decomposer_,
+                distributor_,
+                cellsToRefine,
+                refineParams.maxLoadUnbalance()
+            );
+        }
+        else
+        {
+            meshRefiner_.refineAndBalance
+            (
+                "small feature refinement iteration " + name(iter),
+                decomposer_,
+                distributor_,
+                cellsToRefine,
+                refineParams.maxLoadUnbalance()
+            );
+        }
+    }
+    return iter;
+}
+
+
 Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
 (
     const refinementParameters& refineParams,
@@ -218,7 +341,9 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
                 false,              // internalRefinement
                 true,               // surfaceRefinement
                 true,               // curvatureRefinement
+                false,              // smallFeatureRefinement
                 false,              // gapRefinement
+                false,              // bigGapRefinement
                 refineParams.maxGlobalCells(),
                 refineParams.maxLocalCells()
             )
@@ -352,7 +477,9 @@ Foam::label Foam::autoRefineDriver::gapOnlyRefine
                 false,              // internalRefinement
                 false,              // surfaceRefinement
                 false,              // curvatureRefinement
+                false,              // smallFeatureRefinement
                 true,               // gapRefinement
+                false,              // bigGapRefinement
                 refineParams.maxGlobalCells(),
                 refineParams.maxLocalCells()
             )
@@ -527,6 +654,146 @@ Foam::label Foam::autoRefineDriver::gapOnlyRefine
 }
 
 
+Foam::label Foam::autoRefineDriver::bigGapOnlyRefine
+(
+    const refinementParameters& refineParams,
+    const label maxIter
+)
+{
+    const fvMesh& mesh = meshRefiner_.mesh();
+
+    label iter = 0;
+
+    // See if any surface has an extendedGapLevel
+    labelList surfaceMaxLevel(meshRefiner_.surfaces().maxGapLevel());
+    labelList shellMaxLevel(meshRefiner_.shells().maxGapLevel());
+
+    label overallMaxLevel(max(max(surfaceMaxLevel), max(shellMaxLevel)));
+
+    if (overallMaxLevel == 0)
+    {
+        return iter;
+    }
+
+
+    for (; iter < maxIter; iter++)
+    {
+        Info<< nl
+            << "Big gap refinement iteration " << iter << nl
+            << "------------------------------" << nl
+            << endl;
+
+
+        // Determine cells to refine
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        labelList candidateCells
+        (
+            meshRefiner_.refineCandidates
+            (
+                refineParams.locationsInMesh(),
+                refineParams.curvature(),
+                refineParams.planarAngle(),
+
+                false,              // featureRefinement
+                false,              // featureDistanceRefinement
+                false,              // internalRefinement
+                false,              // surfaceRefinement
+                false,              // curvatureRefinement
+                false,              // smallFeatureRefinement
+                false,              // gapRefinement
+                true,               // bigGapRefinement
+                refineParams.maxGlobalCells(),
+                refineParams.maxLocalCells()
+            )
+        );
+
+
+        if (debug&meshRefinement::MESH)
+        {
+            Pout<< "Dumping " << candidateCells.size()
+                << " cells to cellSet candidateCellsFromBigGap." << endl;
+            cellSet c(mesh, "candidateCellsFromBigGap", 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 >= overallMaxLevel
+             && 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
+            (
+                "big gap refinement iteration " + name(iter),
+                decomposer_,
+                distributor_,
+                cellsToRefine,
+                refineParams.maxLoadUnbalance()
+            );
+        }
+        else
+        {
+            meshRefiner_.refineAndBalance
+            (
+                "big gap refinement iteration " + name(iter),
+                decomposer_,
+                distributor_,
+                cellsToRefine,
+                refineParams.maxLoadUnbalance()
+            );
+        }
+    }
+    return iter;
+}
+
+
 Foam::label Foam::autoRefineDriver::danglingCellRefine
 (
     const refinementParameters& refineParams,
@@ -1108,7 +1375,9 @@ Foam::label Foam::autoRefineDriver::shellRefine
                 true,               // internalRefinement
                 false,              // surfaceRefinement
                 false,              // curvatureRefinement
+                false,              // smallFeatureRefinement
                 false,              // gapRefinement
+                false,              // bigGapRefinement
                 refineParams.maxGlobalCells(),
                 refineParams.maxLocalCells()
             )
@@ -1630,6 +1899,13 @@ void Foam::autoRefineDriver::doRefine
         0       // min cells to refine
     );
 
+    // Refine cells that contain a gap
+    smallFeatureRefine
+    (
+        refineParams,
+        100     // maxIter
+    );
+
     // Refine based on surface
     surfaceOnlyRefine
     (
@@ -1650,6 +1926,13 @@ void Foam::autoRefineDriver::doRefine
         1       // nBufferLayers
     );
 
+    // Refine consistently across narrow gaps (a form of shell refinement)
+    bigGapOnlyRefine
+    (
+        refineParams,
+        100     // maxIter
+    );
+
     // Internal mesh refinement
     shellRefine
     (

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H

@@ -84,6 +84,13 @@ class autoRefineDriver
             const label minRefine
         );
 
+        //- Refine all cells containing small surface features
+        label smallFeatureRefine
+        (
+            const refinementParameters& refineParams,
+            const label maxIter
+        );
+
         //- Refine all cells interacting with the surface
         label surfaceOnlyRefine
         (
@@ -98,6 +105,13 @@ class autoRefineDriver
             const label maxIter
         );
 
+        //- Refine all cells in large gaps
+        label bigGapOnlyRefine
+        (
+            const refinementParameters& refineParams,
+            const label maxIter
+        );
+
         //- Refine cells with almost all sides refined
         label danglingCellRefine
         (

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C

@@ -215,6 +215,55 @@ void Foam::meshRefinement::calcNeighbourData
 }
 
 
+void Foam::meshRefinement::calcCellCellRays
+(
+    const pointField& neiCc,
+    const labelList& neiLevel,
+    const labelList& testFaces,
+    pointField& start,
+    pointField& end,
+    labelList& minLevel
+) const
+{
+    const labelList& cellLevel = meshCutter_.cellLevel();
+    const pointField& cellCentres = mesh_.cellCentres();
+
+    start.setSize(testFaces.size());
+    end.setSize(testFaces.size());
+    minLevel.setSize(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(ROOTSMALL*(end-start));
+        start -= smallVec;
+        end += smallVec;
+    }
+}
+
+
 // Find intersections of edges (given by their two endpoints) with surfaces.
 // Returns first intersection if there are more than one.
 void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
@@ -1187,7 +1236,8 @@ Foam::meshRefinement::meshRefinement
     const bool overwrite,
     const refinementSurfaces& surfaces,
     const refinementFeatures& features,
-    const shellSurfaces& shells
+    const shellSurfaces& shells,
+    const shellSurfaces& limitShells
 )
 :
     mesh_(mesh),
@@ -1197,6 +1247,7 @@ Foam::meshRefinement::meshRefinement
     surfaces_(surfaces),
     features_(features),
     shells_(shells),
+    limitShells_(limitShells),
     meshCutter_
     (
         mesh,

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H

@@ -34,6 +34,7 @@ Description
 SourceFiles
     meshRefinement.C
     meshRefinementBaffles.C
+    meshRefinementGapRefine.C
     meshRefinementMerge.C
     meshRefinementProblemCells.C
     meshRefinementRefine.C
@@ -53,6 +54,7 @@ SourceFiles
 #include "pointIndexHit.H"
 #include "wordPairHashTable.H"
 #include "surfaceZonesInfo.H"
+#include "volumeType.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -176,6 +178,9 @@ private:
         //- All shell-refinement interaction
         const shellSurfaces& shells_;
 
+        //- All limit-refinement interaction
+        const shellSurfaces& limitShells_;
+
         //- Refinement engine
         hexRef8 meshCutter_;
 
@@ -226,10 +231,18 @@ private:
         );
 
         //- Calculate coupled boundary end vector and refinement level
-        void calcNeighbourData
+        void calcNeighbourData(labelList& neiLevel, pointField& neiCc) const;
+
+        //- Calculate rays from cell-centre to cell-centre and corresponding
+        //  min cell refinement level
+        void calcCellCellRays
         (
-            labelList& neiLevel,
+            const pointField& neiCc,
-            pointField& neiCc
+            const labelList& neiLevel,
+            const labelList& testFaces,
+            pointField& start,
+            pointField& end,
+            labelList& minLevel
         ) const;
 
         //- Remove cells. Put exposedFaces into exposedPatchIDs.
@@ -241,12 +254,12 @@ private:
             removeCells& cellRemover
         );
 
-        // Get cells which are inside any closed surface. Note that
+        //- Get cells which are inside any closed surface. Note that
-        // all closed surfaces
+        //  all closed surfaces
-        // will have already been oriented to have keepPoint outside.
+        //  will have already been oriented to have keepPoint outside.
         labelList getInsideCells(const word&) const;
 
-        // Do all to remove inside cells
+        //- Do all to remove inside cells
         autoPtr<mapPolyMesh> removeInsideCells
         (
             const string& msg,
@@ -302,6 +315,14 @@ private:
                 label& nRefine
             ) const;
 
+            //- Unmark cells for refinement based on limit-shells. Return number
+            //  of limited cells.
+            label unmarkInternalRefinement
+            (
+                labelList& refineCell,
+                label& nRefine
+            ) const;
+
             //- Collect faces that are intersected and whose neighbours aren't
             //  yet marked  for refinement.
             labelList getRefineCandidateFaces
@@ -319,6 +340,83 @@ private:
                 label& nRefine
             ) const;
 
+            //- Mark cells intersected by the surface if they are inside
+            //  close gaps
+            label markSurfaceGapRefinement
+            (
+                const scalar planarCos,
+                const label nAllowRefine,
+                const labelList& neiLevel,
+                const pointField& neiCc,
+
+                labelList& refineCell,
+                label& nRefine
+            ) const;
+
+            //- Generate single ray for cell to detect gap
+            bool generateRay
+            (
+                const bool useSurfaceNormal,
+
+                const point& nearPoint,
+                const vector& nearNormal,
+                const FixedList<label, 3>& gapInfo,
+                const volumeType& mode,
+
+                const point& cc,
+                const label cLevel,
+
+                point& start,
+                point& end,
+                point& start2,
+                point& end2
+            ) const;
+
+            //- Select candidate cells (cells inside a shell with gapLevel
+            //  specified)
+            void selectGapCandidates
+            (
+                const labelList& refineCell,
+                const label nRefine,
+
+                labelList& cellMap,
+                List<FixedList<label, 3> >& shellGapInfo,
+                List<volumeType>& shellGapMode
+            ) const;
+
+            //- Merge gap information coming from shell and from surface
+            //  (surface wins)
+            void mergeGapInfo
+            (
+                const FixedList<label, 3>& shellGapInfo,
+                const volumeType shellGapMode,
+                const FixedList<label, 3>& surfGapInfo,
+                const volumeType surfGapMode,
+                FixedList<label, 3>& gapInfo,
+                volumeType& gapMode
+            ) const;
+
+            //- Mark cells for non-surface intersection based gap refinement
+            label markInternalGapRefinement
+            (
+                const scalar planarCos,
+                const label nAllowRefine,
+                labelList& refineCell,
+                label& nRefine
+            ) const;
+
+            //- Refine cells containing small gaps
+            label markSmallFeatureRefinement
+            (
+                const scalar planarCos,
+                const label nAllowRefine,
+                const labelList& neiLevel,
+                const pointField& neiCc,
+
+                labelList& refineCell,
+                label& nRefine
+            ) const;
+
             //- Helper: count number of normals1 that are in normals2
             label countMatches
             (
@@ -387,6 +485,7 @@ private:
             void getIntersections
             (
                 const labelList& surfacesToTest,
+                const labelList& neiLevel,
                 const pointField& neiCc,
                 const labelList& testFaces,
 
@@ -567,7 +666,7 @@ private:
                 const labelList& cellToZone,
                 const labelList& neiCellZone,
                 const labelList& faceToZone,
-                const boolList& meshFlipMap,
+                const PackedBoolList& meshFlipMap,
                 polyTopoChange& meshMod
             ) const;
 
@@ -631,7 +730,7 @@ private:
                 const labelList& nMasterFaces,
                 const labelList& faceToZone,
                 const Map<label>& zoneToOrientation,
-                boolList& meshFlipMap
+                PackedBoolList& meshFlipMap
             ) const;
 
 
@@ -657,6 +756,7 @@ public:
             const bool overwrite,
             const refinementSurfaces&,
             const refinementFeatures&,
+            const shellSurfaces&,
             const shellSurfaces&
         );
 
@@ -855,7 +955,9 @@ public:
                 const bool internalRefinement,
                 const bool surfaceRefinement,
                 const bool curvatureRefinement,
+                const bool smallFeatureRefinement,
                 const bool gapRefinement,
+                const bool bigGapRefinement,
                 const label maxGlobalCells,
                 const label maxLocalCells
             ) const;

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C

@@ -175,6 +175,7 @@ Foam::label Foam::meshRefinement::createBaffle
 void Foam::meshRefinement::getIntersections
 (
     const labelList& surfacesToTest,
+    const labelList& neiLevel,
     const pointField& neiCc,
     const labelList& testFaces,
 
@@ -198,8 +199,6 @@ void Foam::meshRefinement::getIntersections
             << str().name() << nl << endl;
     }
 
-    const pointField& cellCentres = mesh_.cellCentres();
-
 
     globalRegion1.setSize(mesh_.nFaces());
     globalRegion1 = -1;
@@ -213,29 +212,17 @@ void Foam::meshRefinement::getIntersections
     pointField start(testFaces.size());
     pointField end(testFaces.size());
 
-    forAll(testFaces, i)
     {
-        label faceI = testFaces[i];
+        labelList minLevel;
-
+        calcCellCellRays
-        label own = mesh_.faceOwner()[faceI];
+        (
-
+            neiCc,
-        if (mesh_.isInternalFace(faceI))
+            neiLevel,
-        {
+            testFaces,
-            start[i] = cellCentres[own];
+            start,
-            end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
+            end,
-        }
+            minLevel
-        else
+        );
-        {
-            start[i] = cellCentres[own];
-            end[i] = neiCc[faceI-mesh_.nInternalFaces()];
-        }
-    }
-
-    // Extend segments a bit
-    {
-        const vectorField smallVec(ROOTSMALL*(end-start));
-        start -= smallVec;
-        end += smallVec;
     }
 
 
@@ -319,6 +306,7 @@ void Foam::meshRefinement::getBafflePatches
         getIntersections
         (
             surfaceZonesInfo::getUnnamedSurfaces(surfaces_.surfZones()),
+            neiLevel,
             neiCc,
             testFaces,
             globalRegion1,
@@ -2608,7 +2596,7 @@ void Foam::meshRefinement::consistentOrientation
     const labelList& nMasterFacesPerEdge,
     const labelList& faceToZone,
     const Map<label>& zoneToOrientation,
-    boolList& meshFlipMap
+    PackedBoolList& meshFlipMap
 ) const
 {
     const polyBoundaryMesh& bm = mesh_.boundaryMesh();
@@ -2858,7 +2846,7 @@ void Foam::meshRefinement::zonify
     const labelList& cellToZone,
     const labelList& neiCellZone,
     const labelList& faceToZone,
-    const boolList& meshFlipMap,
+    const PackedBoolList& meshFlipMap,
     polyTopoChange& meshMod
 ) const
 {
@@ -3902,6 +3890,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
         getIntersections
         (
             identity(surfaces_.surfaces().size()),  // surfacesToTest,
+            neiLevel,
             neiCc,
             intersectedFaces(),     // testFaces
             globalRegion1,
@@ -4297,7 +4286,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
         }
 
 
-        // 2. Combine faceZoneNames allocated on different processors
+        // 2.Combine faceZoneNames allocated on different processors
 
         Pstream::mapCombineGather(zonesToFaceZone, eqOp<word>());
         Pstream::mapCombineScatter(zonesToFaceZone);
@@ -4438,7 +4427,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
     //      - do a consistent orientation
     //      - check number of faces with consistent orientation
     //      - if <0 flip the whole patch
-    boolList meshFlipMap(mesh_.nFaces(), false);
+    PackedBoolList meshFlipMap(mesh_.nFaces(), false);
     {
         // Collect all data on zone faces without cellZones on either side.
         const indirectPrimitivePatch patch

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C

@@ -742,7 +742,7 @@ Foam::label Foam::meshRefinement::markInternalDistanceToFeatureRefinement
         }
     }
 
-    // Do test to see whether cells is inside/outside shell with higher level
+    // Do test to see whether cells are inside/outside shell with higher level
     labelList maxLevel;
     features_.findHigherLevel(testCc, testLevels, maxLevel);
 
@@ -820,7 +820,7 @@ Foam::label Foam::meshRefinement::markInternalRefinement
         }
     }
 
-    // Do test to see whether cells is inside/outside shell with higher level
+    // Do test to see whether cells are inside/outside shell with higher level
     labelList maxLevel;
     shells_.findHigherLevel(testCc, testLevels, maxLevel);
 
@@ -869,6 +869,56 @@ Foam::label Foam::meshRefinement::markInternalRefinement
 }
 
 
+Foam::label Foam::meshRefinement::unmarkInternalRefinement
+(
+    labelList& refineCell,
+    label& nRefine
+) const
+{
+    const labelList& cellLevel = meshCutter_.cellLevel();
+    const pointField& cellCentres = mesh_.cellCentres();
+
+    label oldNRefine = nRefine;
+
+    // Collect cells to test
+    pointField testCc(nRefine);
+    labelList testLevels(nRefine);
+    label testI = 0;
+
+    forAll(cellLevel, cellI)
+    {
+        if (refineCell[cellI] >= 0)
+        {
+            testCc[testI] = cellCentres[cellI];
+            testLevels[testI] = cellLevel[cellI];
+            testI++;
+        }
+    }
+
+    // Do test to see whether cells are inside/outside shell with higher level
+    labelList levelShell;
+    limitShells_.findLevel(testCc, testLevels, levelShell);
+
+    // Mark for refinement. Note that we didn't store the original cellID so
+    // now just reloop in same order.
+    testI = 0;
+    forAll(cellLevel, cellI)
+    {
+        if (refineCell[cellI] >= 0)
+        {
+            if (levelShell[testI] != -1)
+            {
+                refineCell[cellI] = -1;
+                nRefine--;
+            }
+            testI++;
+        }
+    }
+
+    return returnReduce(oldNRefine-nRefine, sumOp<label>());
+}
+
+
 // Collect faces that are intersected and whose neighbours aren't yet marked
 // for refinement.
 Foam::labelList Foam::meshRefinement::getRefineCandidateFaces
@@ -922,7 +972,6 @@ Foam::label Foam::meshRefinement::markSurfaceRefinement
 ) const
 {
     const labelList& cellLevel = meshCutter_.cellLevel();
-    const pointField& cellCentres = mesh_.cellCentres();
 
     label oldNRefine = nRefine;
 
@@ -941,36 +990,15 @@ Foam::label Foam::meshRefinement::markSurfaceRefinement
     pointField end(testFaces.size());
     labelList minLevel(testFaces.size());
 
-    forAll(testFaces, i)
+    calcCellCellRays
-    {
+    (
-        label faceI = testFaces[i];
+        neiCc,
-
+        neiLevel,
-        label own = mesh_.faceOwner()[faceI];
+        testFaces,
-
+        start,
-        if (mesh_.isInternalFace(faceI))
+        end,
-        {
+        minLevel
-            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(ROOTSMALL*(end-start));
-        start -= smallVec;
-        end += smallVec;
-    }
 
 
     // Do test for higher intersections
@@ -1130,6 +1158,7 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
     pointField end(testFaces.size());
     labelList minLevel(testFaces.size());
 
+    // Note: uses isMasterFace otherwise could be call to calcCellCellRays
     forAll(testFaces, i)
     {
         label faceI = testFaces[i];
@@ -1719,7 +1748,6 @@ Foam::label Foam::meshRefinement::markProximityRefinement
 ) const
 {
     const labelList& cellLevel = meshCutter_.cellLevel();
-    const pointField& cellCentres = mesh_.cellCentres();
 
     label oldNRefine = nRefine;
 
@@ -1739,36 +1767,15 @@ Foam::label Foam::meshRefinement::markProximityRefinement
     pointField end(testFaces.size());
     labelList minLevel(testFaces.size());
 
-    forAll(testFaces, i)
+    calcCellCellRays
-    {
+    (
-        label faceI = testFaces[i];
+        neiCc,
-
+        neiLevel,
-        label own = mesh_.faceOwner()[faceI];
+        testFaces,
-
+        start,
-        if (mesh_.isInternalFace(faceI))
+        end,
-        {
+        minLevel
-            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(ROOTSMALL*(end-start));
-        start -= smallVec;
-        end += smallVec;
-    }
 
 
     // Test for all intersections (with surfaces of higher gap level than
@@ -2045,7 +2052,9 @@ Foam::labelList Foam::meshRefinement::refineCandidates
     const bool internalRefinement,
     const bool surfaceRefinement,
     const bool curvatureRefinement,
+    const bool smallFeatureRefinement,
     const bool gapRefinement,
+    const bool bigGapRefinement,
     const label maxGlobalCells,
     const label maxLocalCells
 ) const
@@ -2096,6 +2105,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
         calcNeighbourData(neiLevel, neiCc);
 
 
+        const scalar planarCos = Foam::cos(degToRad(planarAngle));
+
 
         // Cells pierced by feature lines
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -2163,6 +2174,30 @@ Foam::labelList Foam::meshRefinement::refineCandidates
             );
             Info<< "Marked for refinement due to surface intersection          "
                 << ": " << nSurf << " cells."  << endl;
+
+            // Refine intersected-cells only inside gaps. See
+            // markInternalGapRefinement to refine all cells inside gaps.
+            if
+            (
+                planarCos >= -1
+             && planarCos <= 1
+             && max(shells_.maxGapLevel()) > 0
+            )
+            {
+                label nGapSurf = markSurfaceGapRefinement
+                (
+                    planarCos,
+                    nAllowRefine,
+                    neiLevel,
+                    neiCc,
+
+                    refineCell,
+                    nRefine
+                );
+                Info<< "Marked for refinement due to surface intersection"
+                    << " (at gaps)"
+                    << ": " << nGapSurf << " cells."  << endl;
+            }
         }
 
         // Refinement based on curvature of surface
@@ -2190,7 +2225,33 @@ Foam::labelList Foam::meshRefinement::refineCandidates
         }
 
 
-        const scalar planarCos = Foam::cos(degToRad(planarAngle));
+        // Refinement based on features smaller than cell size
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        if
+        (
+            smallFeatureRefinement
+         && (planarCos >= -1 && planarCos <= 1)
+         && max(shells_.maxGapLevel()) > 0
+        )
+        {
+            label nGap = markSmallFeatureRefinement
+            (
+                planarCos,
+                nAllowRefine,
+                neiLevel,
+                neiCc,
+
+                refineCell,
+                nRefine
+            );
+            Info<< "Marked for refinement due to close opposite surfaces       "
+                << ": " << nGap << " cells."  << endl;
+        }
+
+
+        // Refinement based on gap (only neighbouring cells)
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
         if
         (
@@ -2217,6 +2278,44 @@ Foam::labelList Foam::meshRefinement::refineCandidates
         }
 
 
+        // Refinement based on gaps larger than cell size
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        if
+        (
+            bigGapRefinement
+         && (planarCos >= -1 && planarCos <= 1)
+         && max(shells_.maxGapLevel()) > 0
+        )
+        {
+            // Refine based on gap information provided by shell and nearest
+            // surface
+            label nGap = markInternalGapRefinement
+            (
+                planarCos,
+                nAllowRefine,
+                refineCell,
+                nRefine
+            );
+            Info<< "Marked for refinement due to opposite surfaces"
+                << "             "
+                << ": " << nGap << " cells."  << endl;
+        }
+
+
+        // Limit refinement
+        // ~~~~~~~~~~~~~~~~
+
+        {
+            label nUnmarked = unmarkInternalRefinement(refineCell, nRefine);
+            if (nUnmarked > 0)
+            {
+                Info<< "Unmarked for refinement due to limit shells"
+                    << "                : " << nUnmarked << " cells."  << endl;
+            }
+        }
+
+
 
         // Pack cells-to-refine
         // ~~~~~~~~~~~~~~~~~~~~

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -157,11 +157,22 @@ Foam::refinementSurfaces::refinementSurfaces
     labelList globalMinLevel(surfI, 0);
     labelList globalMaxLevel(surfI, 0);
     labelList globalLevelIncr(surfI, 0);
+
+    FixedList<label, 3> nullGapLevel;
+    nullGapLevel[0] = 0;
+    nullGapLevel[1] = 0;
+    nullGapLevel[2] = 0;
+
+    List<FixedList<label, 3> > globalGapLevel(surfI);
+    List<volumeType> globalGapMode(surfI);
+
     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<FixedList<label, 3> > > regionGapLevel(surfI);
+    List<Map<volumeType> > regionGapMode(surfI);
     List<Map<scalar> > regionAngle(surfI);
     List<Map<autoPtr<dictionary> > > regionPatchInfo(surfI);
 
@@ -212,6 +223,44 @@ Foam::refinementSurfaces::refinementSurfaces
                     << exit(FatalIOError);
             }
 
+
+            // Optional gapLevel specification
+
+            globalGapLevel[surfI] = dict.lookupOrDefault
+            (
+                "gapLevel",
+                nullGapLevel
+            );
+            globalGapMode[surfI] = volumeType::names
+            [
+                dict.lookupOrDefault<word>
+                (
+                    "gapMode",
+                    volumeType::names[volumeType::MIXED]
+                )
+            ];
+            if
+            (
+                globalGapMode[surfI] == volumeType::UNKNOWN
+             || globalGapLevel[surfI][0] < 0
+             || globalGapLevel[surfI][1] < 0
+             || globalGapLevel[surfI][2] < 0
+             || globalGapLevel[surfI][1] > globalGapLevel[surfI][2]
+            )
+            {
+                FatalIOErrorIn
+                (
+                    "refinementSurfaces::refinementSurfaces"
+                    "(const searchableSurfaces&, const dictionary>&",
+                    dict
+                )   << "Illegal gapLevel specification for surface "
+                    << names_[surfI]
+                    << " : gapLevel:" << globalGapLevel[surfI]
+                    << " gapMode:" << volumeType::names[globalGapMode[surfI]]
+                    << exit(FatalIOError);
+            }
+
+
             const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
 
             // Surface zones
@@ -275,6 +324,54 @@ Foam::refinementSurfaces::refinementSurfaces
                                 << exit(FatalIOError);
                         }
 
+
+
+                        // Optional gapLevel specification
+
+                        FixedList<label, 3> gapSpec
+                        (
+                            regionDict.lookupOrDefault
+                            (
+                                "gapLevel",
+                                nullGapLevel
+                            )
+                        );
+                        regionGapLevel[surfI].insert(regionI, gapSpec);
+                        volumeType gapModeSpec
+                        (
+                            volumeType::names
+                            [
+                                regionDict.lookupOrDefault<word>
+                                (
+                                    "gapMode",
+                                    volumeType::names[volumeType::MIXED]
+                                )
+                            ]
+                        );
+                        regionGapMode[surfI].insert(regionI, gapModeSpec);
+                        if
+                        (
+                            gapModeSpec == volumeType::UNKNOWN
+                         || gapSpec[0] < 0
+                         || gapSpec[1] < 0
+                         || gapSpec[2] < 0
+                         || gapSpec[1] > gapSpec[2]
+                        )
+                        {
+                            FatalIOErrorIn
+                            (
+                                "refinementSurfaces::refinementSurfaces"
+                                "(const searchableSurfaces&,"
+                                " const dictionary>&",
+                                dict
+                            )   << "Illegal gapLevel specification for surface "
+                                << names_[surfI]
+                                << " : gapLevel:" << gapSpec
+                                << " gapMode:" << volumeType::names[gapModeSpec]
+                                << exit(FatalIOError);
+                        }
+
+
                         if (regionDict.found("perpendicularAngle"))
                         {
                             regionAngle[surfI].insert
@@ -331,6 +428,10 @@ Foam::refinementSurfaces::refinementSurfaces
     maxLevel_ = 0;
     gapLevel_.setSize(nRegions);
     gapLevel_ = -1;
+    extendedGapLevel_.setSize(nRegions);
+    extendedGapLevel_ = nullGapLevel;
+    extendedGapMode_.setSize(nRegions);
+    extendedGapMode_ = volumeType::UNKNOWN;
     perpendicularAngle_.setSize(nRegions);
     perpendicularAngle_ = -GREAT;
     patchInfo_.setSize(nRegions);
@@ -349,7 +450,8 @@ Foam::refinementSurfaces::refinementSurfaces
             gapLevel_[globalRegionI] =
                 maxLevel_[globalRegionI]
               + globalLevelIncr[surfI];
-
+            extendedGapLevel_[globalRegionI] = globalGapLevel[surfI];
+            extendedGapMode_[globalRegionI] = globalGapMode[surfI];
             perpendicularAngle_[globalRegionI] = globalAngle[surfI];
             if (globalPatchInfo.set(surfI))
             {
@@ -371,6 +473,10 @@ Foam::refinementSurfaces::refinementSurfaces
             gapLevel_[globalRegionI] =
                 maxLevel_[globalRegionI]
               + regionLevelIncr[surfI][iter.key()];
+            extendedGapLevel_[globalRegionI] =
+                regionGapLevel[surfI][iter.key()];
+            extendedGapMode_[globalRegionI] =
+                regionGapMode[surfI][iter.key()];
         }
         forAllConstIter(Map<scalar>, regionAngle[surfI], iter)
         {
@@ -482,12 +588,28 @@ Foam::refinementSurfaces::refinementSurfaces
 // }
 
 
+Foam::labelList Foam::refinementSurfaces::maxGapLevel() const
+{
+    labelList surfaceMax(surfaces_.size(), 0);
+
+    forAll(surfaces_, surfI)
+    {
+        const wordList& regionNames = allGeometry_[surfaces_[surfI]].regions();
+
+        forAll(regionNames, regionI)
+        {
+            label globalI = globalRegion(surfI, regionI);
+            const FixedList<label, 3>& gapInfo = extendedGapLevel_[globalI];
+            surfaceMax[surfI] = max(surfaceMax[surfI], gapInfo[2]);
+        }
+    }
+    return surfaceMax;
+}
+
+
 // Precalculate the refinement level for every element of the searchable
 // surface.
-void Foam::refinementSurfaces::setMinLevelFields
+void Foam::refinementSurfaces::setMinLevelFields(const shellSurfaces& shells)
-(
-    const shellSurfaces& shells
-)
 {
     forAll(surfaces_, surfI)
     {
@@ -1208,6 +1330,49 @@ void Foam::refinementSurfaces::findNearestIntersection
 }
 
 
+void Foam::refinementSurfaces::findNearestIntersection
+(
+    const pointField& start,
+    const pointField& end,
+
+    labelList& surface1,
+    vectorField& normal1
+) const
+{
+    // Initialize arguments
+    surface1.setSize(start.size());
+    surface1 = -1;
+    normal1.setSize(start.size());
+    normal1 = vector::zero;
+
+    // Current end of segment to test.
+    pointField nearest(end);
+    // Work array
+    List<pointIndexHit> nearestInfo(start.size());
+    labelList region;
+    vectorField normal;
+
+    forAll(surfaces_, surfI)
+    {
+        const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
+
+        // See if any intersection between start and current nearest
+        geom.findLine(start, nearest, nearestInfo);
+        geom.getNormal(nearestInfo, normal);
+
+        forAll(nearestInfo, pointI)
+        {
+            if (nearestInfo[pointI].hit())
+            {
+                surface1[pointI] = surfI;
+                normal1[pointI] = normal[pointI];
+                nearest[pointI] = nearestInfo[pointI].hitPoint();
+            }
+        }
+    }
+}
+
+
 void Foam::refinementSurfaces::findAnyIntersection
 (
     const pointField& start,

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H

@@ -42,6 +42,7 @@ SourceFiles
 #include "vectorList.H"
 #include "pointIndexHit.H"
 #include "surfaceZonesInfo.H"
+#include "volumeType.H"
 #include "pointList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -85,6 +86,12 @@ class refinementSurfaces
         //- From global region number to small-gap level
         labelList gapLevel_;
 
+        //- From global region number to small-gap level specification
+        List<FixedList<label, 3> > extendedGapLevel_;
+
+        //- From global region number to side of surface to detect
+        List<volumeType> extendedGapMode_;
+
         //- From global region number to perpendicular angle
         scalarField perpendicularAngle_;
 
@@ -188,6 +195,23 @@ public:
                 return gapLevel_;
             }
 
+            //- From global region number to specification of gap and its
+            //  refinement: 3 labels specifying
+            //  - minimum wanted number of cells in the gap
+            //  - minimum cell level when to start trying to detect gaps
+            //  - maximum cell level to refine to (so do not detect gaps if
+            //    cell >= maximum level)
+            const List<FixedList<label, 3> >& extendedGapLevel() const
+            {
+                return extendedGapLevel_;
+            }
+
+            //- From global region number to side of surface to detect
+            const List<volumeType>& extendedGapMode() const
+            {
+                return extendedGapMode_;
+            }
+
             //- From global region number to perpendicular angle
             const scalarField& perpendicularAngle() const
             {
@@ -226,6 +250,9 @@ public:
                 return minLevel_.size();
             }
 
+            //- Per surface the maximum extendedGapLevel over all its regions
+            labelList maxGapLevel() const;
+
             //- Calculate minLevelFields
             void setMinLevelFields
             (
@@ -314,6 +341,15 @@ public:
                 vectorField& normal2
             ) const;
 
+            //- Find nearest (to start only) intersection of edge
+            void findNearestIntersection
+            (
+                const pointField& start,
+                const pointField& end,
+                labelList& surfaces,
+                vectorField& normal
+            ) const;
+
             //- Used for debugging only: find intersection of edge.
             void findAnyIntersection
             (

src/mesh/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -124,7 +124,7 @@ void Foam::shellSurfaces::setAndCheckLevels
     }
     else
     {
-        if (!allGeometry_[shells_[shellI]].hasVolumeType())
+        if (!shell.hasVolumeType())
         {
             FatalErrorIn
             (
@@ -152,6 +152,46 @@ void Foam::shellSurfaces::setAndCheckLevels
 }
 
 
+void Foam::shellSurfaces::checkGapLevels
+(
+    const dictionary& shellDict,
+    const label shellI,
+    const List<FixedList<label, 3> >& levels
+)
+{
+    const searchableSurface& shell = allGeometry_[shells_[shellI]];
+
+    forAll(levels, regionI)
+    {
+        const FixedList<label, 3>& info = levels[regionI];
+
+        if (info[2] > 0)
+        {
+            if (modes_[shellI] == DISTANCE)
+            {
+                FatalIOErrorIn
+                (
+                    "shellSurfaces::shellSurfaces(..)",
+                    shellDict
+                )   << "'gapLevel' specification cannot be used with mode "
+                    << refineModeNames_[DISTANCE]
+                    << " for shell " << shell.name()
+                    << exit(FatalIOError);
+            }
+        }
+    }
+
+    // Hardcode for region 0
+    if (levels[0][0] > 0)
+    {
+        Info<< "Refinement level up to " << levels[0][2]
+            << " for all cells in gaps for shell "
+            << shell.name() << endl;
+    }
+}
+
+
+
 // Specifically orient triSurfaces using a calculated point outside.
 // Done since quite often triSurfaces not of consistent orientation which
 // is (currently) necessary for sideness calculation
@@ -288,18 +328,18 @@ void Foam::shellSurfaces::findHigherLevel
         );
 
         // Update maxLevel
-        forAll(nearInfo, candidateI)
+        forAll(nearInfo, i)
         {
-            if (nearInfo[candidateI].hit())
+            if (nearInfo[i].hit())
             {
                 // Check which level it actually is in.
                 label minDistI = findLower
                 (
                     distances,
-                    mag(nearInfo[candidateI].hitPoint()-candidates[candidateI])
+                    mag(nearInfo[i].hitPoint()-candidates[i])
                 );
 
-                label pointI = candidateMap[candidateI];
+                label pointI = candidateMap[i];
 
                 // pt is inbetween shell[minDistI] and shell[minDistI+1]
                 maxLevel[pointI] = levels[minDistI+1];
@@ -356,6 +396,197 @@ void Foam::shellSurfaces::findHigherLevel
 }
 
 
+void Foam::shellSurfaces::findHigherGapLevel
+(
+    const pointField& pt,
+    const labelList& ptLevel,
+    const label shellI,
+    labelList& gapShell,
+    List<FixedList<label, 3> >& gapInfo,
+    List<volumeType>& gapMode
+) const
+{
+    //TBD: hardcoded for region 0 information
+    const FixedList<label, 3>& info = extendedGapLevel_[shellI][0];
+    volumeType mode = extendedGapMode_[shellI][0];
+
+    if (info[2] == 0)
+    {
+        return;
+    }
+
+
+    // Collect all those points that have a current maxLevel less than the
+    // shell.
+
+    labelList candidateMap(pt.size());
+    label candidateI = 0;
+
+    forAll(ptLevel, pointI)
+    {
+        if (ptLevel[pointI] >= info[1] && ptLevel[pointI] < info[2])
+        {
+            candidateMap[candidateI++] = pointI;
+        }
+    }
+    candidateMap.setSize(candidateI);
+
+    // Do the expensive nearest test only for the candidate points.
+    List<volumeType> volType;
+    allGeometry_[shells_[shellI]].getVolumeType
+    (
+        pointField(pt, candidateMap),
+        volType
+    );
+
+    forAll(volType, i)
+    {
+        label pointI = candidateMap[i];
+
+        bool isInside = (volType[i] == volumeType::INSIDE);
+
+        if
+        (
+            (
+                (modes_[shellI] == INSIDE && isInside)
+             || (modes_[shellI] == OUTSIDE && !isInside)
+            )
+         && info[2] > gapInfo[pointI][2]
+        )
+        {
+            gapShell[pointI] = shellI;
+            gapInfo[pointI] = info;
+            gapMode[pointI] = mode;
+        }
+    }
+}
+
+
+void Foam::shellSurfaces::findLevel
+(
+    const pointField& pt,
+    const label shellI,
+    labelList& minLevel,
+    labelList& shell
+) const
+{
+    const labelList& levels = levels_[shellI];
+
+    if (modes_[shellI] == DISTANCE)
+    {
+        // Distance mode.
+
+        const scalarField& distances = distances_[shellI];
+
+        // Collect all those points that have a current level equal/greater
+        // (any of) the shell. Also collect the furthest distance allowable
+        // to any shell with a higher level.
+
+        pointField candidates(pt.size());
+        labelList candidateMap(pt.size());
+        scalarField candidateDistSqr(pt.size());
+        label candidateI = 0;
+
+        forAll(shell, pointI)
+        {
+            if (shell[pointI] == -1)
+            {
+                forAllReverse(levels, levelI)
+                {
+                    if (levels[levelI] <= minLevel[pointI])
+                    {
+                        candidates[candidateI] = pt[pointI];
+                        candidateMap[candidateI] = pointI;
+                        candidateDistSqr[candidateI] = sqr(distances[levelI]);
+                        candidateI++;
+                        break;
+                    }
+                }
+            }
+        }
+        candidates.setSize(candidateI);
+        candidateMap.setSize(candidateI);
+        candidateDistSqr.setSize(candidateI);
+
+        // Do the expensive nearest test only for the candidate points.
+        List<pointIndexHit> nearInfo;
+        allGeometry_[shells_[shellI]].findNearest
+        (
+            candidates,
+            candidateDistSqr,
+            nearInfo
+        );
+
+        // Update maxLevel
+        forAll(nearInfo, i)
+        {
+            if (nearInfo[i].hit())
+            {
+                // Check which level it actually is in.
+                label minDistI = findLower
+                (
+                    distances,
+                    mag(nearInfo[i].hitPoint()-candidates[i])
+                );
+
+                label pointI = candidateMap[i];
+
+                // pt is inbetween shell[minDistI] and shell[minDistI+1]
+                shell[pointI] = shellI;
+                minLevel[pointI] = levels[minDistI+1];
+            }
+        }
+    }
+    else
+    {
+        // Inside/outside mode
+
+        // Collect all those points that have a current maxLevel less than the
+        // shell.
+
+        pointField candidates(pt.size());
+        labelList candidateMap(pt.size());
+        label candidateI = 0;
+
+        forAll(shell, pointI)
+        {
+            if (shell[pointI] == -1 && levels[0] <= minLevel[pointI])
+            {
+                candidates[candidateI] = pt[pointI];
+                candidateMap[candidateI] = pointI;
+                candidateI++;
+            }
+        }
+        candidates.setSize(candidateI);
+        candidateMap.setSize(candidateI);
+
+        // Do the expensive nearest test only for the candidate points.
+        List<volumeType> volType;
+        allGeometry_[shells_[shellI]].getVolumeType(candidates, volType);
+
+        forAll(volType, i)
+        {
+            if
+            (
+                (
+                    modes_[shellI] == INSIDE
+                 && volType[i] == volumeType::INSIDE
+                )
+             || (
+                    modes_[shellI] == OUTSIDE
+                 && volType[i] == volumeType::OUTSIDE
+                )
+            )
+            {
+                label pointI = candidateMap[i];
+                shell[pointI] = shellI;
+                minLevel[pointI] = levels[0];
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::shellSurfaces::shellSurfaces
@@ -387,6 +618,15 @@ Foam::shellSurfaces::shellSurfaces
     distances_.setSize(shellI);
     levels_.setSize(shellI);
 
+    extendedGapLevel_.setSize(shellI);
+    extendedGapMode_.setSize(shellI);
+
+    FixedList<label, 3> nullGapLevel;
+    nullGapLevel[0] = 0;
+    nullGapLevel[1] = 0;
+    nullGapLevel[2] = 0;
+
+
     HashSet<word> unmatchedKeys(shellsDict.toc());
     shellI = 0;
 
@@ -407,6 +647,85 @@ Foam::shellSurfaces::shellSurfaces
             // Read pairs of distance+level
             setAndCheckLevels(shellI, dict.lookup("levels"));
 
+
+
+            // Gap specification
+            // ~~~~~~~~~~~~~~~~~
+
+
+            // Shell-wide gap level specification
+            const searchableSurface& surface = allGeometry_[geomI];
+            const wordList& regionNames = surface.regions();
+
+            FixedList<label, 3> gapSpec
+            (
+                dict.lookupOrDefault
+                (
+                    "gapLevel",
+                    nullGapLevel
+                )
+            );
+            extendedGapLevel_[shellI].setSize(regionNames.size());
+            extendedGapLevel_[shellI] = gapSpec;
+
+            volumeType gapModeSpec
+            (
+                volumeType::names
+                [
+                    dict.lookupOrDefault<word>
+                    (
+                        "gapMode",
+                        volumeType::names[volumeType::MIXED]
+                    )
+                ]
+            );
+            extendedGapMode_[shellI].setSize(regionNames.size());
+            extendedGapMode_[shellI] = gapModeSpec;
+
+
+            // Override on a per-region basis?
+
+            if (dict.found("regions"))
+            {
+                const dictionary& regionsDict = dict.subDict("regions");
+                forAll(regionNames, regionI)
+                {
+                    if (regionsDict.found(regionNames[regionI]))
+                    {
+                        // Get the dictionary for region
+                        const dictionary& regionDict = regionsDict.subDict
+                        (
+                            regionNames[regionI]
+                        );
+                        FixedList<label, 3> gapSpec
+                        (
+                            regionDict.lookupOrDefault
+                            (
+                                "gapLevel",
+                                nullGapLevel
+                            )
+                        );
+                        extendedGapLevel_[shellI][regionI] = gapSpec;
+
+                        volumeType gapModeSpec
+                        (
+                            volumeType::names
+                            [
+                                regionDict.lookupOrDefault<word>
+                                (
+                                    "gapMode",
+                                    volumeType::names[volumeType::MIXED]
+                                )
+                            ]
+                        );
+                        extendedGapMode_[shellI][regionI] = gapModeSpec;
+                    }
+                }
+            }
+
+
+            checkGapLevels(dict, shellI, extendedGapLevel_[shellI]);
+
             shellI++;
         }
     }
@@ -444,6 +763,22 @@ Foam::label Foam::shellSurfaces::maxLevel() const
 }
 
 
+Foam::labelList Foam::shellSurfaces::maxGapLevel() const
+{
+    labelList surfaceMax(extendedGapLevel_.size(), 0);
+
+    forAll(extendedGapLevel_, shellI)
+    {
+        const List<FixedList<label, 3> >& levels = extendedGapLevel_[shellI];
+        forAll(levels, i)
+        {
+            surfaceMax[shellI] = max(surfaceMax[shellI], levels[i][2]);
+        }
+    }
+    return surfaceMax;
+}
+
+
 void Foam::shellSurfaces::findHigherLevel
 (
     const pointField& pt,
@@ -461,4 +796,66 @@ void Foam::shellSurfaces::findHigherLevel
 }
 
 
+void Foam::shellSurfaces::findHigherGapLevel
+(
+    const pointField& pt,
+    const labelList& ptLevel,
+    labelList& gapShell,
+    List<FixedList<label, 3> >& gapInfo,
+    List<volumeType>& gapMode
+) const
+{
+    gapShell.setSize(pt.size());
+    gapShell = -1;
+
+    FixedList<label, 3> nullGapLevel;
+    nullGapLevel[0] = 0;
+    nullGapLevel[1] = 0;
+    nullGapLevel[2] = 0;
+
+    gapInfo.setSize(pt.size());
+    gapInfo = nullGapLevel;
+
+    gapMode.setSize(pt.size());
+    gapMode = volumeType::MIXED;
+
+    forAll(shells_, shellI)
+    {
+        findHigherGapLevel(pt, ptLevel, shellI, gapShell, gapInfo, gapMode);
+    }
+}
+
+
+void Foam::shellSurfaces::findHigherGapLevel
+(
+    const pointField& pt,
+    const labelList& ptLevel,
+    List<FixedList<label, 3> >& gapInfo,
+    List<volumeType>& gapMode
+) const
+{
+    labelList gapShell;
+    findHigherGapLevel(pt, ptLevel, gapShell, gapInfo, gapMode);
+}
+
+
+void Foam::shellSurfaces::findLevel
+(
+    const pointField& pt,
+    const labelList& ptLevel,
+    labelList& shell
+) const
+{
+    shell.setSize(pt.size());
+    shell = -1;
+
+    labelList minLevel(ptLevel);
+
+    forAll(shells_, shellI)
+    {
+        findLevel(pt, shellI, minLevel, shell);
+    }
+}
+
+
 // ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.H

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -85,6 +85,15 @@ private:
         labelListList levels_;
 
 
+        // Gap level refinement
+
+            //- Per shell, per region the small-gap level specification
+            List<List<FixedList<label, 3> > > extendedGapLevel_;
+
+            //- Per shell, per region the small-gap level specification
+            List<List<volumeType> > extendedGapMode_;
+
+
     // Private data
 
         //- refineMode names
@@ -93,15 +102,24 @@ private:
 
     // Private Member Functions
 
-        //- Helper function for initialisation.
+        //- Helper function for initialisation of levels
         void setAndCheckLevels
         (
             const label shellI,
             const List<Tuple2<scalar, label> >&
         );
 
+        //- Helper function for checking of gap information
+        void checkGapLevels
+        (
+            const dictionary&,
+            const label shellI,
+            const List<FixedList<label, 3> >& levels
+        );
+
         void orient();
 
+        //- Find first shell with a level higher than maxLevel
         void findHigherLevel
         (
             const pointField& pt,
@@ -109,6 +127,27 @@ private:
             labelList& maxLevel
         ) const;
 
+        //- Update highest min gap level
+        void findHigherGapLevel
+        (
+            const pointField& pt,
+            const labelList& ptLevel,
+            const label shellI,
+            labelList& gapShell,
+            List<FixedList<label, 3> >& gapInfo,
+            List<volumeType>& gapMode
+        ) const;
+
+        //- Find first shell with a level lower or equal to minLevel. Update
+        //  minLevel and shell.
+        void findLevel
+        (
+            const pointField& pt,
+            const label shellI,
+            labelList& minLevel,
+            labelList& shell
+        ) const;
+
 public:
 
     // Constructors
@@ -125,16 +164,11 @@ public:
 
         // Access
 
-            //const List<scalarField>& distances() const
+            //- Indices of surfaces that are shells
-            //{
+            const labelList& shells() const
-            //    return distances_;
+            {
-            //}
+                return shells_;
-            //
+            }
-            ////- Per shell per distance the refinement level
-            //const labelListList& levels() const
-            //{
-            //    return levels_;
-            //}
 
 
         // Query
@@ -142,6 +176,9 @@ public:
             //- Highest shell level
             label maxLevel() const;
 
+            //- Highest shell gap level
+            labelList maxGapLevel() const;
+
             //- Find shell level higher than ptLevel
             void findHigherLevel
             (
@@ -150,6 +187,33 @@ public:
                 labelList& maxLevel
             ) const;
 
+            //- Find a shell whose minimum gap level is >= ptLevel
+            void findHigherGapLevel
+            (
+                const pointField& pt,
+                const labelList& ptLevel,
+                labelList& gapShell,
+                List<FixedList<label, 3> >& gapInfo,
+                List<volumeType>& gapMode
+            ) const;
+
+            //- Find a shell whose minimum gap level is >= ptLevel. gapInfo
+            //  is (0 0 0) if no shell found
+            void findHigherGapLevel
+            (
+                const pointField& pt,
+                const labelList& ptLevel,
+                List<FixedList<label, 3> >& gapInfo,
+                List<volumeType>& gapMode
+            ) const;
+
+            //- Find first shell (or -1) with level equal or lower than ptLevel.
+            void findLevel
+            (
+                const pointField& pt,
+                const labelList& ptLevel,
+                labelList& shell
+            ) const;
 };
 
 

src/meshTools/Make/files

@@ -63,6 +63,7 @@ indexedOctree/treeDataTriSurface.C
 
 searchableSurface = searchableSurface
 $(searchableSurface)/searchableBox.C
+$(searchableSurface)/searchableRotatedBox.C
 $(searchableSurface)/searchableCylinder.C
 $(searchableSurface)/searchableDisk.C
 $(searchableSurface)/searchablePlane.C
@@ -75,6 +76,7 @@ $(searchableSurface)/searchableSurfacesQueries.C
 $(searchableSurface)/searchableSurfaceWithGaps.C
 $(searchableSurface)/triSurfaceMesh.C
 $(searchableSurface)/closedTriSurfaceMesh.C
+$(searchableSurface)/subTriSurfaceMesh.C
 
 topoSets = sets/topoSets
 $(topoSets)/cellSet.C

src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -408,7 +408,14 @@ void Foam::hierarchGeomDecomp::sortComponent
         {
             // No need for binary searching of bin size
             localSize = label(current.size()/n_[compI]);
-            rightCoord = sortedCoord[leftIndex+localSize];
+            if (leftIndex+localSize < sortedCoord.size())
+            {
+                rightCoord = sortedCoord[leftIndex+localSize];
+            }
+            else
+            {
+                rightCoord = maxCoord;
+            }
         }
         else
         {