diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index 67a9c3a944..7b4d8577c4 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/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;
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index 9b759b5ca4..6b40b2ab5e 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/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
-    // few cells. This setting will cause refinement to stop if <= minimumRefine
-    // are selected for refinement. Note: it will at least do one iteration
-    // (unless the number of cells to refine is 0)
+    // The surface refinement loop might spend lots of iterations refining just
+    // a few cells. This setting will cause refinement to stop if
+    // <= minimumRefine cells are selected for refinement. Note: it will
+    // at least do one iteration (unless the number of cells to refine is 0)
     minRefinementCells 0;
 
     // Allow a certain level of imbalance during refining
@@ -244,11 +244,27 @@ castellatedMeshControls
         }
         //sphere.stl
         //{
-        //    mode distance;
-        //    levels ((1.0 5) (2.0 3));
+        //    mode inside;
+        //    levels ((1.0 4));
+        //    // Optional override of uniform refinement level such
+        //    //  that in small gaps we're getting more cells.
+        //    //  The specification is
+        //    //  - numGapCells : minimum number of cells in the gap
+        //    //  - minLevel    : min refinement level at which to kick in
+        //    //  - maxLevel    : upper refinement level
+        //    // All three settings can be overridden on a surface by
+        //    // surface basis in the refinementSurfaces section.
+        //    gapLevel (<numGapCells> <minLevel> <maxlevel>);
         //}
     }
 
+
+    // Limit refinement in geometric region
+    limitRegions
+    {
+    }
+
+
     // Mesh selection
     // ~~~~~~~~~~~~~~
 
diff --git a/src/Allwmake b/src/Allwmake
index c0ea041c37..d70cb03fc2 100755
--- a/src/Allwmake
+++ b/src/Allwmake
@@ -59,6 +59,8 @@ parallel/Allwmake $targetType $*
 wmake $targetType ODE
 wmake $targetType randomProcesses
 
+wmake $targetType fvMotionSolver
+
 transportModels/Allwmake $targetType $*
 thermophysicalModels/Allwmake $targetType $*
 TurbulenceModels/Allwmake $targetType $*
@@ -69,7 +71,6 @@ mesh/Allwmake $targetType $*
 renumber/Allwmake $targetType $*
 fvAgglomerationMethods/Allwmake $targetType $*
 
-wmake $targetType fvMotionSolver
 wmake $targetType engine
 wmake $targetType fvOptions
 wmake $targetType regionCoupled
diff --git a/src/OpenFOAM/algorithms/indexedOctree/volumeType.H b/src/OpenFOAM/algorithms/indexedOctree/volumeType.H
index 7bdfc043e6..3d845c414b 100644
--- a/src/OpenFOAM/algorithms/indexedOctree/volumeType.H
+++ b/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
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     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
 
diff --git a/src/mesh/autoMesh/Make/files b/src/mesh/autoMesh/Make/files
index 081d70afbd..29001bafb7 100644
--- a/src/mesh/autoMesh/Make/files
+++ b/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
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index be908b76bb..bda68a1708 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -105,7 +105,10 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
                     false,              // internalRefinement
                     false,              // surfaceRefinement
                     false,              // curvatureRefinement
+                    false,              // smallFeatureRefinement
                     false,              // gapRefinement
+                    false,              // bigGapRefinement
+                    false,              // spreadGapSize
                     refineParams.maxGlobalCells(),
                     refineParams.maxLocalCells()
                 )
@@ -179,6 +182,128 @@ 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
+                false,              // spreadGapSize
+                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 +343,10 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
                 false,              // internalRefinement
                 true,               // surfaceRefinement
                 true,               // curvatureRefinement
+                false,              // smallFeatureRefinement
                 false,              // gapRefinement
+                false,              // bigGapRefinement
+                false,              // spreadGapSize
                 refineParams.maxGlobalCells(),
                 refineParams.maxLocalCells()
             )
@@ -352,7 +480,10 @@ Foam::label Foam::autoRefineDriver::gapOnlyRefine
                 false,              // internalRefinement
                 false,              // surfaceRefinement
                 false,              // curvatureRefinement
+                false,              // smallFeatureRefinement
                 true,               // gapRefinement
+                false,              // bigGapRefinement
+                false,              // spreadGapSize
                 refineParams.maxGlobalCells(),
                 refineParams.maxLocalCells()
             )
@@ -527,6 +658,148 @@ Foam::label Foam::autoRefineDriver::gapOnlyRefine
 }
 
 
+Foam::label Foam::autoRefineDriver::bigGapOnlyRefine
+(
+    const refinementParameters& refineParams,
+    const bool spreadGapSize,
+    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
+                spreadGapSize,      // spreadGapSize
+                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 +1381,10 @@ Foam::label Foam::autoRefineDriver::shellRefine
                 true,               // internalRefinement
                 false,              // surfaceRefinement
                 false,              // curvatureRefinement
+                false,              // smallFeatureRefinement
                 false,              // gapRefinement
+                false,              // bigGapRefinement
+                false,              // spreadGapSize
                 refineParams.maxGlobalCells(),
                 refineParams.maxLocalCells()
             )
@@ -1630,6 +1906,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 +1933,14 @@ void Foam::autoRefineDriver::doRefine
         1       // nBufferLayers
     );
 
+    // Refine consistently across narrow gaps (a form of shell refinement)
+    bigGapOnlyRefine
+    (
+        refineParams,
+        true,   // spreadGapSize
+        100     // maxIter
+    );
+
     // Internal mesh refinement
     shellRefine
     (
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
index 11ed1d9454..19f8a03fbd 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
+++ b/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,14 @@ class autoRefineDriver
             const label maxIter
         );
 
+        //- Refine all cells in large gaps
+        label bigGapOnlyRefine
+        (
+            const refinementParameters& refineParams,
+            const bool spreadGapSize,
+            const label maxIter
+        );
+
         //- Refine cells with almost all sides refined
         label danglingCellRefine
         (
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index e4bfa8840c..93818eb8d6 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/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,
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 31ef9f2e1a..2e71d97579 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/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,8 @@ SourceFiles
 #include "pointIndexHit.H"
 #include "wordPairHashTable.H"
 #include "surfaceZonesInfo.H"
+#include "volumeType.H"
+#include "DynamicField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -176,6 +179,9 @@ private:
         //- All shell-refinement interaction
         const shellSurfaces& shells_;
 
+        //- All limit-refinement interaction
+        const shellSurfaces& limitShells_;
+
         //- Refinement engine
         hexRef8 meshCutter_;
 
@@ -226,10 +232,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,
-            pointField& neiCc
+            const pointField& neiCc,
+            const labelList& neiLevel,
+            const labelList& testFaces,
+            pointField& start,
+            pointField& end,
+            labelList& minLevel
         ) const;
 
         //- Remove cells. Put exposedFaces into exposedPatchIDs.
@@ -241,12 +255,12 @@ private:
             removeCells& cellRemover
         );
 
-        // Get cells which are inside any closed surface. Note that
-        // all closed surfaces
-        // will have already been oriented to have keepPoint outside.
+        //- Get cells which are inside any closed surface. Note that
+        //  all closed surfaces
+        //  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 +316,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 +341,104 @@ 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 from nearPoint in direction of nearNormal
+            label generateRays
+            (
+                const point& nearPoint,
+                const vector& nearNormal,
+                const FixedList<label, 3>& gapInfo,
+                const volumeType& mode,
+
+                const label cLevel,
+
+                DynamicField<point>& start,
+                DynamicField<point>& end
+            ) const;
+
+            //- Generate pairs of rays through cell centre
+            //  Each ray pair has start, end, and expected gap size
+            label generateRays
+            (
+                const bool useSurfaceNormal,
+
+                const point& nearPoint,
+                const vector& nearNormal,
+                const FixedList<label, 3>& gapInfo,
+                const volumeType& mode,
+
+                const point& cc,
+                const label cLevel,
+
+                DynamicField<point>& start,
+                DynamicField<point>& end,
+                DynamicField<scalar>& gapSize,
+
+                DynamicField<point>& start2,
+                DynamicField<point>& end2,
+                DynamicField<scalar>& gapSize2
+            ) 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 bool spreadGapSize,
+                const label nAllowRefine,
+                labelList& refineCell,
+                label& nRefine,
+                labelList& numGapCells,
+                scalarField& gapSize
+            ) 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 +507,7 @@ private:
             void getIntersections
             (
                 const labelList& surfacesToTest,
+                const labelList& neiLevel,
                 const pointField& neiCc,
                 const labelList& testFaces,
 
@@ -567,7 +688,7 @@ private:
                 const labelList& cellToZone,
                 const labelList& neiCellZone,
                 const labelList& faceToZone,
-                const boolList& meshFlipMap,
+                const PackedBoolList& meshFlipMap,
                 polyTopoChange& meshMod
             ) const;
 
@@ -631,7 +752,7 @@ private:
                 const labelList& nMasterFaces,
                 const labelList& faceToZone,
                 const Map<label>& zoneToOrientation,
-                boolList& meshFlipMap
+                PackedBoolList& meshFlipMap
             ) const;
 
 
@@ -657,6 +778,7 @@ public:
             const bool overwrite,
             const refinementSurfaces&,
             const refinementFeatures&,
+            const shellSurfaces&,
             const shellSurfaces&
         );
 
@@ -855,7 +977,10 @@ public:
                 const bool internalRefinement,
                 const bool surfaceRefinement,
                 const bool curvatureRefinement,
+                const bool smallFeatureRefinement,
                 const bool gapRefinement,
+                const bool bigGapRefinement,
+                const bool spreadGapSize,
                 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 6f500c117a..302882a3fb 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/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];
-
-        label own = mesh_.faceOwner()[faceI];
-
-        if (mesh_.isInternalFace(faceI))
-        {
-            start[i] = cellCentres[own];
-            end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
-        }
-        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;
+        labelList minLevel;
+        calcCellCellRays
+        (
+            neiCc,
+            neiLevel,
+            testFaces,
+            start,
+            end,
+            minLevel
+        );
     }
 
 
@@ -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
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementGapRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementGapRefine.C
new file mode 100644
index 0000000000..d73ea12a0b
--- /dev/null
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementGapRefine.C
@@ -0,0 +1,1676 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "meshRefinement.H"
+#include "Time.H"
+#include "refinementSurfaces.H"
+#include "refinementFeatures.H"
+#include "shellSurfaces.H"
+#include "triSurfaceMesh.H"
+#include "treeDataCell.H"
+#include "searchableSurfaces.H"
+#include "DynamicField.H"
+#include "transportData.H"
+#include "FaceCellWave.H"
+#include "volFields.H"
+#include "zeroGradientFvPatchFields.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::label Foam::meshRefinement::markSurfaceGapRefinement
+(
+    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();
+
+    // Get the gap level for the shells
+    const labelList maxLevel(shells_.maxGapLevel());
+
+    label oldNRefine = nRefine;
+
+    if (max(maxLevel) > 0)
+    {
+        // Use cached surfaceIndex_ to detect if any intersection. If so
+        // re-intersect to determine level wanted.
+
+        // Collect candidate faces
+        // ~~~~~~~~~~~~~~~~~~~~~~~
+
+        labelList testFaces(getRefineCandidateFaces(refineCell));
+
+        // Collect segments
+        // ~~~~~~~~~~~~~~~~
+
+        pointField start(testFaces.size());
+        pointField end(testFaces.size());
+        {
+            labelList minLevel(testFaces.size());
+            calcCellCellRays
+            (
+                neiCc,
+                neiLevel,
+                testFaces,
+                start,
+                end,
+                minLevel
+            );
+        }
+
+
+        // Collect cells to test for inside/outside in shell
+        labelList cellToCompact(mesh_.nCells(), -1);
+        labelList bFaceToCompact(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
+        List<FixedList<label, 3> > shellGapInfo;
+        List<volumeType> shellGapMode;
+        {
+            DynamicField<point> compactToCc(mesh_.nCells()/10);
+            DynamicList<label> compactToLevel(compactToCc.capacity());
+            forAll(testFaces, i)
+            {
+                label faceI = testFaces[i];
+                label own = mesh_.faceOwner()[faceI];
+                if (cellToCompact[own] == -1)
+                {
+                    cellToCompact[own] = compactToCc.size();
+                    compactToCc.append(cellCentres[own]);
+                    compactToLevel.append(cellLevel[own]);
+                }
+                if (mesh_.isInternalFace(faceI))
+                {
+                    label nei = mesh_.faceNeighbour()[faceI];
+                    if (cellToCompact[nei] == -1)
+                    {
+                        cellToCompact[nei] = compactToCc.size();
+                        compactToCc.append(cellCentres[nei]);
+                        compactToLevel.append(cellLevel[nei]);
+                    }
+                }
+                else
+                {
+                    label bFaceI = faceI - mesh_.nInternalFaces();
+                    if (bFaceToCompact[bFaceI] == -1)
+                    {
+                        bFaceToCompact[bFaceI] = compactToCc.size();
+                        compactToCc.append(neiCc[bFaceI]);
+                        compactToLevel.append(neiLevel[bFaceI]);
+                    }
+                }
+            }
+
+            shells_.findHigherGapLevel
+            (
+                compactToCc,
+                compactToLevel,
+                shellGapInfo,
+                shellGapMode
+            );
+        }
+
+
+        const List<FixedList<label, 3> >& extendedGapLevel =
+            surfaces_.extendedGapLevel();
+        const List<volumeType>& extendedGapMode =
+            surfaces_.extendedGapMode();
+
+        labelList ccSurface1;
+        List<pointIndexHit> ccHit1;
+        labelList ccRegion1;
+        vectorField ccNormal1;
+
+        {
+            labelList ccSurface2;
+            List<pointIndexHit> ccHit2;
+            labelList ccRegion2;
+            vectorField ccNormal2;
+
+            surfaces_.findNearestIntersection
+            (
+                identity(surfaces_.surfaces().size()),
+                start,
+                end,
+
+                ccSurface1,
+                ccHit1,
+                ccRegion1,
+                ccNormal1,
+
+                ccSurface2,
+                ccHit2,
+                ccRegion2,
+                ccNormal2
+            );
+        }
+
+        start.clear();
+        end.clear();
+
+        DynamicField<point> rayStart(2*ccSurface1.size());
+        DynamicField<point> rayEnd(2*ccSurface1.size());
+        DynamicField<scalar> gapSize(2*ccSurface1.size());
+
+        DynamicField<point> rayStart2(2*ccSurface1.size());
+        DynamicField<point> rayEnd2(2*ccSurface1.size());
+        DynamicField<scalar> gapSize2(2*ccSurface1.size());
+
+        DynamicList<label> cellMap(2*ccSurface1.size());
+        DynamicList<label> compactMap(2*ccSurface1.size());
+
+        forAll(ccSurface1, i)
+        {
+            label surfI = ccSurface1[i];
+
+            if (surfI != -1)
+            {
+                label globalRegionI =
+                    surfaces_.globalRegion(surfI, ccRegion1[i]);
+
+                label faceI = testFaces[i];
+                const point& surfPt = ccHit1[i].hitPoint();
+
+                label own = mesh_.faceOwner()[faceI];
+                if (cellToCompact[own] != -1)
+                {
+                    // Combine info from shell and surface
+                    label compactI = cellToCompact[own];
+                    FixedList<label, 3> gapInfo;
+                    volumeType gapMode;
+                    mergeGapInfo
+                    (
+                        shellGapInfo[compactI],
+                        shellGapMode[compactI],
+                        extendedGapLevel[globalRegionI],
+                        extendedGapMode[globalRegionI],
+                        gapInfo,
+                        gapMode
+                    );
+
+                    const point& cc = cellCentres[own];
+                    label nRays = generateRays
+                    (
+                        false,
+                        surfPt,
+                        ccNormal1[i],
+                        gapInfo,
+                        gapMode,
+                        surfPt+((cc-surfPt)&ccNormal1[i])*ccNormal1[i],
+                        cellLevel[own],
+
+                        rayStart,
+                        rayEnd,
+                        gapSize,
+
+                        rayStart2,
+                        rayEnd2,
+                        gapSize2
+                    );
+                    for (label j = 0; j < nRays; j++)
+                    {
+                        cellMap.append(own);
+                        compactMap.append(i);
+                    }
+                }
+                if (mesh_.isInternalFace(faceI))
+                {
+                    label nei = mesh_.faceNeighbour()[faceI];
+                    if (cellToCompact[nei] != -1)
+                    {
+                        // Combine info from shell and surface
+                        label compactI = cellToCompact[nei];
+                        FixedList<label, 3> gapInfo;
+                        volumeType gapMode;
+                        mergeGapInfo
+                        (
+                            shellGapInfo[compactI],
+                            shellGapMode[compactI],
+                            extendedGapLevel[globalRegionI],
+                            extendedGapMode[globalRegionI],
+                            gapInfo,
+                            gapMode
+                        );
+
+                        const point& cc = cellCentres[nei];
+                        label nRays = generateRays
+                        (
+                            false,
+                            surfPt,
+                            ccNormal1[i],
+                            gapInfo,
+                            gapMode,
+                            surfPt+((cc-surfPt)&ccNormal1[i])*ccNormal1[i],
+                            cellLevel[nei],
+
+                            rayStart,
+                            rayEnd,
+                            gapSize,
+
+                            rayStart2,
+                            rayEnd2,
+                            gapSize2
+                        );
+                        for (label j = 0; j < nRays; j++)
+                        {
+                            cellMap.append(nei);
+                            compactMap.append(i);
+                        }
+                    }
+                }
+                else
+                {
+                    // Note: on coupled face. What cell are we going to
+                    // refine? We've got the neighbouring cell centre
+                    // and level but we cannot mark it for refinement on
+                    // this side...
+                    label bFaceI = faceI - mesh_.nInternalFaces();
+
+                    if (bFaceToCompact[bFaceI] != -1)
+                    {
+                        // Combine info from shell and surface
+                        label compactI = bFaceToCompact[bFaceI];
+                        FixedList<label, 3> gapInfo;
+                        volumeType gapMode;
+                        mergeGapInfo
+                        (
+                            shellGapInfo[compactI],
+                            shellGapMode[compactI],
+                            extendedGapLevel[globalRegionI],
+                            extendedGapMode[globalRegionI],
+                            gapInfo,
+                            gapMode
+                        );
+
+                        const point& cc = neiCc[bFaceI];
+                        label nRays = generateRays
+                        (
+                            false,
+                            surfPt,
+                            ccNormal1[i],
+                            gapInfo,
+                            gapMode,
+                            surfPt+((cc-surfPt)&ccNormal1[i])*ccNormal1[i],
+                            neiLevel[bFaceI],
+
+                            rayStart,
+                            rayEnd,
+                            gapSize,
+
+                            rayStart2,
+                            rayEnd2,
+                            gapSize2
+                        );
+                        for (label j = 0; j < nRays; j++)
+                        {
+                            cellMap.append(-1); // See above.
+                            compactMap.append(i);
+                        }
+                    }
+                }
+            }
+        }
+
+        Info<< "Shooting " << returnReduce(rayStart.size(), sumOp<label>())
+            << " rays from " << returnReduce(testFaces.size(), sumOp<label>())
+            << " intersected faces" << endl;
+
+        rayStart.shrink();
+        rayEnd.shrink();
+        gapSize.shrink();
+
+        rayStart2.shrink();
+        rayEnd2.shrink();
+        gapSize2.shrink();
+
+        cellMap.shrink();
+        compactMap.shrink();
+
+        testFaces.clear();
+        ccSurface1.clear();
+        ccHit1.clear();
+        ccRegion1.clear();
+        ccNormal1 = UIndirectList<vector>(ccNormal1, compactMap)();
+
+
+        // Do intersections in pairs
+        labelList surf1;
+        List<pointIndexHit> hit1;
+        vectorField normal1;
+        surfaces_.findNearestIntersection
+        (
+            rayStart,
+            rayEnd,
+            surf1,
+            hit1,
+            normal1
+        );
+
+        labelList surf2;
+        List<pointIndexHit> hit2;
+        vectorField normal2;
+        surfaces_.findNearestIntersection
+        (
+            rayStart2,
+            rayEnd2,
+            surf2,
+            hit2,
+            normal2
+        );
+
+        forAll(surf1, i)
+        {
+            if (surf1[i] != -1 && surf2[i] != -1)
+            {
+                // Found intersection with surface. Check opposite normal.
+                label cellI = cellMap[i];
+
+                if
+                (
+                    cellI != -1
+                 && (mag(normal1[i]&normal2[i]) > planarCos)
+                 && (
+                        magSqr(hit1[i].hitPoint()-hit2[i].hitPoint())
+                      < Foam::sqr(gapSize[i])
+                    )
+                )
+                {
+                    if
+                    (
+                       !markForRefine
+                        (
+                            surf1[i],
+                            nAllowRefine,
+                            refineCell[cellI],
+                            nRefine
+                        )
+                    )
+                    {
+                        break;
+                    }
+                }
+            }
+        }
+
+        if
+        (
+            returnReduce(nRefine, sumOp<label>())
+          > returnReduce(nAllowRefine, sumOp<label>())
+        )
+        {
+            Info<< "Reached refinement limit." << endl;
+        }
+    }
+
+    return returnReduce(nRefine-oldNRefine, sumOp<label>());
+}
+
+
+//Foam::meshRefinement::findNearestOppositeOp::findNearestOppositeOp
+//(
+//    const indexedOctree<treeDataTriSurface>& tree,
+//    const point& oppositePoint,
+//    const vector& oppositeNormal,
+//    const scalar minCos
+//)
+//:
+//    tree_(tree),
+//    oppositePoint_(oppositePoint),
+//    oppositeNormal_(oppositeNormal),
+//    minCos_(minCos)
+//{}
+//
+//
+//void Foam::meshRefinement::findNearestOppositeOp::operator()
+//(
+//    const labelUList& indices,
+//    const point& sample,
+//    scalar& nearestDistSqr,
+//    label& minIndex,
+//    point& nearestPoint
+//) const
+//{
+//    const treeDataTriSurface& shape = tree_.shapes();
+//    const triSurface& patch = shape.patch();
+//    const pointField& points = patch.points();
+//
+//    forAll(indices, i)
+//    {
+//        const label index = indices[i];
+//        const labelledTri& f = patch[index];
+//
+//        pointHit nearHit = f.nearestPoint(sample, points);
+//        scalar distSqr = sqr(nearHit.distance());
+//
+//        if (distSqr < nearestDistSqr)
+//        {
+//            // Nearer. Check if
+//            // - a bit way from other hit
+//            // - in correct search cone
+//            vector d(nearHit.rawPoint()-oppositePoint_);
+//            scalar normalDist(d&oppositeNormal_);
+//
+//            if (normalDist > Foam::sqr(SMALL) && normalDist/mag(d) > minCos_)
+//            {
+//                nearestDistSqr = distSqr;
+//                minIndex = index;
+//                nearestPoint = nearHit.rawPoint();
+//            }
+//        }
+//    }
+//}
+//
+//
+//void Foam::meshRefinement::searchCone
+//(
+//    const label surfI,
+//    labelList& nearMap,                 // cells
+//    scalarField& nearGap,               // gap size
+//    List<pointIndexHit>& nearInfo,      // nearest point on surface
+//    List<pointIndexHit>& oppositeInfo   // detected point on gap (or miss)
+//) const
+//{
+//    const labelList& cellLevel = meshCutter_.cellLevel();
+//    const pointField& cellCentres = mesh_.cellCentres();
+//    const scalar edge0Len = meshCutter_.level0EdgeLength();
+//
+//    const labelList& surfaceIndices = surfaces_.surfaces();
+//    const List<FixedList<label, 3> >& extendedGapLevel =
+//        surfaces_.extendedGapLevel();
+//    const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
+//
+//
+//    label geomI = surfaceIndices[surfI];
+//    const searchableSurface& geom = surfaces_.geometry()[geomI];
+//
+//    const triSurfaceMesh& s = refCast<const triSurfaceMesh>(geom);
+//    const indexedOctree<treeDataTriSurface>& tree = s.tree();
+//
+//
+//    const scalar searchCos(Foam::cos(degToRad(30)));
+//
+//    // Normals for ray shooting and inside/outside detection
+//    vectorField nearNormal;
+//    geom.getNormal(nearInfo, nearNormal);
+//    // Regions
+//    labelList nearRegion;
+//    geom.getRegion(nearInfo, nearRegion);
+//
+//
+//    // Now loop over all near points and search in the half cone
+//    labelList map(nearInfo.size());
+//    label compactI = 0;
+//
+//    oppositeInfo.setSize(nearInfo.size());
+//
+//    forAll(nearInfo, i)
+//    {
+//        label globalRegionI =
+//            surfaces_.globalRegion(surfI, nearRegion[i]);
+//
+//        // Get updated gap information now we have the region
+//        label nGapCells = extendedGapLevel[globalRegionI][0];
+//        label minLevel = extendedGapLevel[globalRegionI][1];
+//        label maxLevel = extendedGapLevel[globalRegionI][2];
+//        volumeType mode = extendedGapMode[globalRegionI];
+//
+//        label cellI = nearMap[i];
+//        label cLevel = cellLevel[cellI];
+//
+//        if (cLevel >= minLevel && cLevel < maxLevel)
+//        {
+//            scalar cellSize = edge0Len/pow(2, cLevel);
+//
+//            // Update gap size
+//            nearGap[i] = nGapCells*cellSize;
+//
+//            const point& nearPt = nearInfo[i].hitPoint();
+//            vector v(cellCentres[cellI]-nearPt);
+//            scalar magV = mag(v);
+//
+//            // Like with ray shooting we want to
+//            // - find triangles up to nearGap away on the wanted side of the
+//            //   surface
+//            // - find triangles up to 0.5*cellSize away on the unwanted side
+//            //   of the surface. This is for cells straddling the surface
+//            //   where
+//            //   the cell centre might be on the wrong side of the surface
+//
+//            // Tbd: check that cell centre is inbetween the gap hits
+//            // (only if the cell is far enough away)
+//
+//            scalar posNormalSize = 0.0;
+//            scalar negNormalSize = 0.0;
+//
+//            if (mode == volumeType::OUTSIDE)
+//            {
+//                posNormalSize = nearGap[i];
+//                if (magV < 0.5*cellSize)
+//                {
+//                    negNormalSize = 0.5*cellSize;
+//                }
+//            }
+//            else if (mode == volumeType::INSIDE)
+//            {
+//                if (magV < 0.5*cellSize)
+//                {
+//                    posNormalSize = 0.5*cellSize;
+//                }
+//                negNormalSize = nearGap[i];
+//            }
+//            else
+//            {
+//                posNormalSize = nearGap[i];
+//                negNormalSize = nearGap[i];
+//            }
+//
+//            // Test with positive normal
+//            oppositeInfo[compactI] = tree.findNearest
+//            (
+//                nearPt,
+//                sqr(posNormalSize),
+//                findNearestOppositeOp
+//                (
+//                    tree,
+//                    nearPt,
+//                    nearNormal[i],
+//                    searchCos
+//                )
+//            );
+//
+//            if (oppositeInfo[compactI].hit())
+//            {
+//                map[compactI++] = i;
+//            }
+//            else
+//            {
+//                // Test with negative normal
+//                oppositeInfo[compactI] = tree.findNearest
+//                (
+//                    nearPt,
+//                    sqr(negNormalSize),
+//                    findNearestOppositeOp
+//                    (
+//                        tree,
+//                        nearPt,
+//                        -nearNormal[i],
+//                        searchCos
+//                    )
+//                );
+//
+//                if (oppositeInfo[compactI].hit())
+//                {
+//                    map[compactI++] = i;
+//                }
+//            }
+//        }
+//    }
+//
+//    Info<< "Selected " << returnReduce(compactI, sumOp<label>())
+//        << " hits on the correct side out of "
+//        << returnReduce(map.size(), sumOp<label>()) << endl;
+//    map.setSize(compactI);
+//    oppositeInfo.setSize(compactI);
+//
+//    nearMap = UIndirectList<label>(nearMap, map)();
+//    nearGap = UIndirectList<scalar>(nearGap, map)();
+//    nearInfo = UIndirectList<pointIndexHit>(nearInfo, map)();
+//    nearNormal = UIndirectList<vector>(nearNormal, map)();
+//
+//    // Exclude hits which aren't opposite enough. E.g. you might find
+//    // a point on a perpendicular wall - but this does not consistute a gap.
+//    vectorField oppositeNormal;
+//    geom.getNormal(oppositeInfo, oppositeNormal);
+//
+//    compactI = 0;
+//    forAll(oppositeInfo, i)
+//    {
+//        if ((nearNormal[i] & oppositeNormal[i]) < -0.707)
+//        {
+//            map[compactI++] = i;
+//        }
+//    }
+//
+//    Info<< "Selected " << returnReduce(compactI, sumOp<label>())
+//        << " hits opposite the nearest out of "
+//        << returnReduce(map.size(), sumOp<label>()) << endl;
+//    map.setSize(compactI);
+//
+//    nearMap = UIndirectList<label>(nearMap, map)();
+//    nearGap = UIndirectList<scalar>(nearGap, map)();
+//    nearInfo = UIndirectList<pointIndexHit>(nearInfo, map)();
+//    oppositeInfo = UIndirectList<pointIndexHit>(oppositeInfo, map)();
+//}
+
+
+Foam::label Foam::meshRefinement::generateRays
+(
+    const point& nearPoint,
+    const vector& nearNormal,
+    const FixedList<label, 3>& gapInfo,
+    const volumeType& mode,
+
+    const label cLevel,
+
+    DynamicField<point>& start,
+    DynamicField<point>& end
+) const
+{
+    label nOldRays = start.size();
+
+    if (cLevel >= gapInfo[1] && cLevel < gapInfo[2])
+    {
+        scalar cellSize = meshCutter_.level0EdgeLength()/pow(2, cLevel);
+
+        // Calculate gap size
+        scalar nearGap = gapInfo[0]*cellSize;
+
+        const vector& n = nearNormal;
+
+        // Situation 'C' above: cell too close. Use surface
+        // -normal and -point to shoot rays
+
+        if (mode == volumeType::OUTSIDE)
+        {
+            start.append(nearPoint+1e-6*n);
+            end.append(nearPoint+nearGap*n);
+        }
+        else if (mode == volumeType::INSIDE)
+        {
+            start.append(nearPoint-1e-6*n);
+            end.append(nearPoint-nearGap*n);
+        }
+        else if (mode == volumeType::MIXED)
+        {
+            start.append(nearPoint+1e-6*n);
+            end.append(nearPoint+nearGap*n);
+
+            start.append(nearPoint-1e-6*n);
+            end.append(nearPoint-nearGap*n);
+        }
+    }
+
+    return start.size()-nOldRays;
+}
+
+
+Foam::label Foam::meshRefinement::generateRays
+(
+    const bool useSurfaceNormal,
+
+    const point& nearPoint,
+    const vector& nearNormal,
+    const FixedList<label, 3>& gapInfo,
+    const volumeType& mode,
+
+    const point& cc,
+    const label cLevel,
+
+    DynamicField<point>& start,
+    DynamicField<point>& end,
+    DynamicField<scalar>& gapSize,
+
+    DynamicField<point>& start2,
+    DynamicField<point>& end2,
+    DynamicField<scalar>& gapSize2
+) const
+{
+    // We want to handle the following cases:
+    // - surface: small gap (marked with 'surface'). gap might be
+    //            on inside or outside of surface.
+    // - A: cell well inside the gap.
+    // - B: cell well outside the gap.
+    // - C: cell straddling the gap. cell centre might be inside
+    //      or outside
+    //
+    //       +---+
+    //       | B |
+    //       +---+
+    //
+    //            +------+
+    //            |      |
+    //            |   C  |
+    //    --------|------|----surface
+    //            +------+
+    //
+    //        +---+
+    //        | A |
+    //        +---+
+    //
+    //
+    //    --------------------surface
+    //
+    // So:
+    // - find nearest point on surface
+    // - in situation A,B decide if on wanted side of surface
+    // - detect if locally a gap (and the cell inside the gap) by
+    //   shooting a ray from the point on the surface in the direction
+    //   of
+    //   - A,B: the cell centre
+    //   - C: the surface normal and/or negative surface normal
+    //   and see we hit anything
+    //
+    // Variations of this scheme:
+    // - always shoot in the direction of the surface normal. This needs
+    //   then an additional check to make sure the cell centre is
+    //   somewhere inside the gap
+    // - instead of ray shooting use a 'constrained' nearest search
+    //   by e.g. looking inside a search cone (implemented in searchCone).
+    //   The problem with this constrained nearest is that it still uses
+    //   the absolute nearest point on each triangle and only afterwards
+    //   checks if it is inside the search cone.
+
+
+    // Decide which near points are good:
+    // - with updated minLevel and maxLevel and nearGap make sure
+    //   the cell is still a candidate
+    //   NOTE: inside the gap the nearest point on the surface will
+    //         be HALF the gap size - otherwise we would have found
+    //         a point on the opposite side
+    // - if the mode is both sides
+    // - or if the hit is inside the current cell (situation 'C',
+    //   magV < 0.5cellSize)
+    // - or otherwise if on the correct side
+
+    label nOldRays = start.size();
+
+    if (cLevel >= gapInfo[1] && cLevel < gapInfo[2])
+    {
+        scalar cellSize = meshCutter_.level0EdgeLength()/pow(2, cLevel);
+
+        // Calculate gap size
+        scalar nearGap = gapInfo[0]*cellSize;
+
+        // Distance to nearest
+        vector v(cc-nearPoint);
+        scalar magV = mag(v);
+
+        if (useSurfaceNormal || magV < 0.5*cellSize)
+        {
+            const vector& n = nearNormal;
+
+            // Situation 'C' above: cell too close. Use surface
+            // -normal and -point to shoot rays
+
+            if (mode == volumeType::OUTSIDE)
+            {
+                start.append(nearPoint+1e-6*n);
+                end.append(nearPoint+nearGap*n);
+                gapSize.append(nearGap);
+                // Second vector so we get pairs of intersections
+                start2.append(nearPoint+1e-6*n);
+                end2.append(nearPoint-1e-6*n);
+                gapSize2.append(gapSize.last());
+            }
+            else if (mode == volumeType::INSIDE)
+            {
+                start.append(nearPoint-1e-6*n);
+                end.append(nearPoint-nearGap*n);
+                gapSize.append(nearGap);
+                // Second vector so we get pairs of intersections
+                start2.append(nearPoint-1e-6*n);
+                end2.append(nearPoint+1e-6*n);
+                gapSize2.append(gapSize.last());
+            }
+            else if (mode == volumeType::MIXED)
+            {
+                // Do both rays:
+                // Outside
+                {
+                    start.append(nearPoint+1e-6*n);
+                    end.append(nearPoint+nearGap*n);
+                    gapSize.append(nearGap);
+                    // Second vector so we get pairs of intersections
+                    start2.append(nearPoint+1e-6*n);
+                    end2.append(nearPoint-1e-6*n);
+                    gapSize2.append(gapSize.last());
+                }
+                // Inside
+                {
+                    start.append(nearPoint-1e-6*n);
+                    end.append(nearPoint-nearGap*n);
+                    gapSize.append(nearGap);
+                    // Second vector so we get pairs of intersections
+                    start2.append(nearPoint-1e-6*n);
+                    end2.append(nearPoint+1e-6*n);
+                    gapSize2.append(gapSize.last());
+                }
+            }
+        }
+        else
+        {
+            // Situation 'A' or 'B' above: cell well away. Test if
+            // cell on correct side of surface and shoot ray through
+            // cell centre. Note: no need to shoot ray in other
+            // direction since we're trying to detect cell inside
+            // the gap.
+
+            scalar s = (v&nearNormal);
+
+            if
+            (
+                (mode == volumeType::MIXED)
+             || (mode == volumeType::OUTSIDE && s > SMALL)
+             || (mode == volumeType::INSIDE && s < -SMALL)
+            )
+            {
+                //// Use single vector through cell centre
+                //vector n(v/(magV+ROOTVSMALL));
+                //
+                //start.append(cc);
+                //end.append(cc+nearGap*n);
+                //gapSize.append(nearGap);
+                //
+                //start2.append(cc);
+                //end2.append(cc-nearGap*n);
+                //gapSize2.append(nearGap);
+
+
+                //// Shoot some rays through the cell centre
+                //// X-direction:
+                //start.append(cc);
+                //end.append(cc+nearGap*vector(1, 0, 0));
+                //gapSize.append(nearGap);
+                //
+                //start2.append(cc);
+                //end2.append(cc-nearGap*vector(1, 0, 0));
+                //gapSize2.append(nearGap);
+                //
+                //// Y-direction:
+                //start.append(cc);
+                //end.append(cc+nearGap*vector(0, 1, 0));
+                //gapSize.append(nearGap);
+                //
+                //start2.append(cc);
+                //end2.append(cc-nearGap*vector(0, 1, 0));
+                //gapSize2.append(nearGap);
+                //
+                //// Z-direction:
+                //start.append(cc);
+                //end.append(cc+nearGap*vector(0, 0, 1));
+                //gapSize.append(nearGap);
+                //
+                //start2.append(cc);
+                //end2.append(cc-nearGap*vector(0, 0, 1));
+                //gapSize2.append(nearGap);
+
+
+                // 3 axes aligned with normal
+
+                // Use vector through cell centre
+                vector n(v/(magV+ROOTVSMALL));
+
+                // Get second vector. Make sure it is sufficiently perpendicular
+                vector e2(1, 0, 0);
+                scalar s = (e2 & n);
+                if (mag(s) < 0.9)
+                {
+                    e2 -= s*n;
+                }
+                else
+                {
+                    e2 = vector(0, 1, 0);
+                    e2 -= (e2 & n)*n;
+                }
+                e2 /= mag(e2);
+
+                // Third vector
+                vector e3 = n ^ e2;
+
+
+                // Rays in first direction
+                start.append(cc);
+                end.append(cc+nearGap*n);
+                gapSize.append(nearGap);
+
+                start2.append(cc);
+                end2.append(cc-nearGap*n);
+                gapSize2.append(nearGap);
+
+                // Rays in second direction
+                start.append(cc);
+                end.append(cc+nearGap*e2);
+                gapSize.append(nearGap);
+
+                start2.append(cc);
+                end2.append(cc-nearGap*e2);
+                gapSize2.append(nearGap);
+
+                // Rays in third direction
+                start.append(cc);
+                end.append(cc+nearGap*e3);
+                gapSize.append(nearGap);
+
+                start2.append(cc);
+                end2.append(cc-nearGap*e3);
+                gapSize2.append(nearGap);
+            }
+        }
+    }
+
+    return start.size()-nOldRays;
+}
+
+
+void Foam::meshRefinement::selectGapCandidates
+(
+    const labelList& refineCell,
+    const label nRefine,
+
+    labelList& cellMap,
+    List<FixedList<label, 3> >& shellGapInfo,
+    List<volumeType>& shellGapMode
+) const
+{
+    const labelList& cellLevel = meshCutter_.cellLevel();
+    const pointField& cellCentres = mesh_.cellCentres();
+
+    // Collect cells to test
+    cellMap.setSize(cellLevel.size()-nRefine);
+    label compactI = 0;
+
+    forAll(cellLevel, cellI)
+    {
+        if (refineCell[cellI] == -1)
+        {
+            cellMap[compactI++] = cellI;
+        }
+    }
+    Info<< "Selected " << returnReduce(compactI, sumOp<label>())
+        << " unmarked cells out of "
+        << mesh_.globalData().nTotalCells() << endl;
+    cellMap.setSize(compactI);
+
+    // Do test to see whether cells are inside/outside shell with
+    // applicable specification (minLevel <= celllevel < maxLevel)
+    shells_.findHigherGapLevel
+    (
+        pointField(cellCentres, cellMap),
+        UIndirectList<label>(cellLevel, cellMap)(),
+        shellGapInfo,
+        shellGapMode
+    );
+
+    // Compact out hits
+
+    labelList map(shellGapInfo.size());
+    compactI = 0;
+    forAll(shellGapInfo, i)
+    {
+        if (shellGapInfo[i][2] > 0)
+        {
+            map[compactI++] = i;
+        }
+    }
+
+    Info<< "Selected " << returnReduce(compactI, sumOp<label>())
+        << " cells inside gap shells out of "
+        << mesh_.globalData().nTotalCells() << endl;
+
+    map.setSize(compactI);
+    cellMap = UIndirectList<label>(cellMap, map)();
+    shellGapInfo = UIndirectList<FixedList<label, 3> >(shellGapInfo, map)();
+    shellGapMode = UIndirectList<volumeType>(shellGapMode, map)();
+}
+
+
+void Foam::meshRefinement::mergeGapInfo
+(
+    const FixedList<label, 3>& shellGapInfo,
+    const volumeType shellGapMode,
+    const FixedList<label, 3>& surfGapInfo,
+    const volumeType surfGapMode,
+    FixedList<label, 3>& gapInfo,
+    volumeType& gapMode
+) const
+{
+    if (surfGapInfo[0] == 0)
+    {
+        gapInfo = shellGapInfo;
+        gapMode = shellGapMode;
+    }
+    else if (shellGapInfo[0] == 0)
+    {
+        gapInfo = surfGapInfo;
+        gapMode = surfGapMode;
+    }
+    else
+    {
+        // Both specify a level. Does surface level win? Or does information
+        // need to be merged?
+
+        //gapInfo[0] = max(surfGapInfo[0], shellGapInfo[0]);
+        //gapInfo[1] = min(surfGapInfo[1], shellGapInfo[1]);
+        //gapInfo[2] = max(surfGapInfo[2], shellGapInfo[2]);
+        gapInfo = surfGapInfo;
+        gapMode = surfGapMode;
+    }
+}
+
+
+Foam::label Foam::meshRefinement::markInternalGapRefinement
+(
+    const scalar planarCos,
+    const bool spreadGapSize,
+    const label nAllowRefine,
+
+    labelList& refineCell,
+    label& nRefine,
+    labelList& numGapCells,
+    scalarField& detectedGapSize
+) const
+{
+    detectedGapSize.setSize(mesh_.nCells());
+    detectedGapSize = GREAT;
+    numGapCells.setSize(mesh_.nCells());
+    numGapCells = -1;
+
+    const labelList& cellLevel = meshCutter_.cellLevel();
+    const pointField& cellCentres = mesh_.cellCentres();
+    const scalar edge0Len = meshCutter_.level0EdgeLength();
+
+    const List<FixedList<label, 3> >& extendedGapLevel =
+        surfaces_.extendedGapLevel();
+    const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
+
+    // Get the gap level for the shells
+    const labelList maxLevel(shells_.maxGapLevel());
+
+    label oldNRefine = nRefine;
+
+    if (max(maxLevel) > 0)
+    {
+        // Collect cells to test
+        labelList cellMap;
+        List<FixedList<label, 3> > shellGapInfo;
+        List<volumeType> shellGapMode;
+        selectGapCandidates
+        (
+            refineCell,
+            nRefine,
+
+            cellMap,
+            shellGapInfo,
+            shellGapMode
+        );
+
+        // Find nearest point and normal on the surfaces
+        List<pointIndexHit> nearInfo;
+        vectorField nearNormal;
+        labelList nearSurface;
+        labelList nearRegion;
+        {
+            // Now we have both the cell-level and the gap size information. Use
+            // this to calculate the gap size
+            scalarField gapSize(cellMap.size());
+            forAll(cellMap, i)
+            {
+                label cellI = cellMap[i];
+                scalar cellSize = edge0Len/pow(2, cellLevel[cellI]);
+                gapSize[i] = shellGapInfo[i][0]*cellSize;
+            }
+
+            surfaces_.findNearestRegion
+            (
+                identity(surfaces_.surfaces().size()),
+                pointField(cellCentres, cellMap),
+                sqr(gapSize),
+                nearSurface,
+                nearInfo,
+                nearRegion,
+                nearNormal
+            );
+        }
+
+
+
+        DynamicList<label> map(nearInfo.size());
+        DynamicField<point> rayStart(nearInfo.size());
+        DynamicField<point> rayEnd(nearInfo.size());
+        DynamicField<scalar> gapSize(nearInfo.size());
+
+        DynamicField<point> rayStart2(nearInfo.size());
+        DynamicField<point> rayEnd2(nearInfo.size());
+        DynamicField<scalar> gapSize2(nearInfo.size());
+
+        label nTestCells = 0;
+
+        forAll(nearInfo, i)
+        {
+            if (nearInfo[i].hit())
+            {
+                label globalRegionI = surfaces_.globalRegion
+                (
+                    nearSurface[i],
+                    nearRegion[i]
+                );
+
+                // Combine info from shell and surface
+                FixedList<label, 3> gapInfo;
+                volumeType gapMode;
+                mergeGapInfo
+                (
+                    shellGapInfo[i],
+                    shellGapMode[i],
+                    extendedGapLevel[globalRegionI],
+                    extendedGapMode[globalRegionI],
+                    gapInfo,
+                    gapMode
+                );
+
+                // Store wanted number of cells in gap
+                label cellI = cellMap[i];
+                label cLevel = cellLevel[cellI];
+                if (cLevel >= gapInfo[1] && cLevel < gapInfo[2])
+                {
+                    numGapCells[cellI] = max(numGapCells[cellI], gapInfo[0]);
+                }
+
+                // Construct one or more rays to test for oppositeness
+                label nRays = generateRays
+                (
+                    false,
+                    nearInfo[i].hitPoint(),
+                    nearNormal[i],
+                    gapInfo,
+                    gapMode,
+
+                    cellCentres[cellI],
+                    cLevel,
+
+                    rayStart,
+                    rayEnd,
+                    gapSize,
+
+                    rayStart2,
+                    rayEnd2,
+                    gapSize2
+                );
+                if (nRays > 0)
+                {
+                    nTestCells++;
+                    for (label j = 0; j < nRays; j++)
+                    {
+                        map.append(i);
+                    }
+                }
+            }
+        }
+
+        Info<< "Selected " << returnReduce(nTestCells, sumOp<label>())
+            << " cells for testing out of "
+            << mesh_.globalData().nTotalCells() << endl;
+        map.shrink();
+        rayStart.shrink();
+        rayEnd.shrink();
+        gapSize.shrink();
+
+        rayStart2.shrink();
+        rayEnd2.shrink();
+        gapSize2.shrink();
+
+        cellMap = UIndirectList<label>(cellMap, map)();
+        nearNormal = UIndirectList<vector>(nearNormal, map)();
+        shellGapInfo.clear();
+        shellGapMode.clear();
+        nearInfo.clear();
+        nearSurface.clear();
+        nearRegion.clear();
+
+
+        // Do intersections in pairs
+        labelList surf1;
+        List<pointIndexHit> hit1;
+        vectorField normal1;
+        surfaces_.findNearestIntersection
+        (
+            rayStart,
+            rayEnd,
+            surf1,
+            hit1,
+            normal1
+        );
+
+        labelList surf2;
+        List<pointIndexHit> hit2;
+        vectorField normal2;
+        surfaces_.findNearestIntersection
+        (
+            rayStart2,
+            rayEnd2,
+            surf2,
+            hit2,
+            normal2
+        );
+
+        // Extract cell based gap size
+        forAll(surf1, i)
+        {
+            if (surf1[i] != -1 && surf2[i] != -1)
+            {
+                // Found intersections with surface. Check for
+                // - small gap
+                // - coplanar normals
+
+                label cellI = cellMap[i];
+
+                scalar d2 = magSqr(hit1[i].hitPoint()-hit2[i].hitPoint());
+
+                if
+                (
+                    cellI != -1
+                 && (mag(normal1[i]&normal2[i]) > planarCos)
+                 && (d2 < Foam::sqr(gapSize[i]))
+                )
+                {
+                    detectedGapSize[cellI] = min
+                    (
+                        detectedGapSize[cellI],
+                        Foam::sqrt(d2)
+                    );
+                }
+            }
+        }
+
+        // Spread it
+        if (spreadGapSize)
+        {
+            // Field on cells and faces
+            List<transportData> cellData(mesh_.nCells());
+            List<transportData> faceData(mesh_.nFaces());
+
+            // Start of walk
+            const pointField& faceCentres = mesh_.faceCentres();
+
+            DynamicList<label> frontFaces(mesh_.nFaces());
+            DynamicList<transportData> frontData(mesh_.nFaces());
+            for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+            {
+                label own = mesh_.faceOwner()[faceI];
+                label nei = mesh_.faceNeighbour()[faceI];
+
+                scalar minSize = min
+                (
+                    detectedGapSize[own],
+                    detectedGapSize[nei]
+                );
+
+                if (minSize < GREAT)
+                {
+                    frontFaces.append(faceI);
+                    frontData.append
+                    (
+                        transportData
+                        (
+                            faceCentres[faceI],
+                            minSize,
+                            0.0
+                        )
+                    );
+                }
+            }
+            for
+            (
+                label faceI = mesh_.nInternalFaces();
+                faceI < mesh_.nFaces();
+                faceI++
+            )
+            {
+                label own = mesh_.faceOwner()[faceI];
+
+                if (detectedGapSize[own] < GREAT)
+                {
+                    frontFaces.append(faceI);
+                    frontData.append
+                    (
+                        transportData
+                        (
+                            faceCentres[faceI],
+                            detectedGapSize[own],
+                            0.0
+                        )
+                    );
+                }
+            }
+
+            Info<< "Selected "
+                << returnReduce(frontFaces.size(), sumOp<label>())
+                << " faces for spreading gap size out of "
+                << mesh_.globalData().nTotalFaces() << endl;
+
+
+            transportData::trackData td(surfaceIndex_);
+
+            FaceCellWave<transportData, transportData::trackData> deltaCalc
+            (
+                mesh_,
+                frontFaces,
+                frontData,
+                faceData,
+                cellData,
+                mesh_.globalData().nTotalCells()+1,
+                td
+            );
+
+
+            forAll(cellMap, i)
+            {
+                label cellI = cellMap[i];
+                if
+                (
+                    cellI != -1
+                 && cellData[cellI].valid(deltaCalc.data())
+                 && numGapCells[cellI] != -1
+                )
+                {
+                    // Update transported gap size
+                    detectedGapSize[cellI] = min
+                    (
+                        detectedGapSize[cellI],
+                        cellData[cellI].data()
+                    );
+                }
+            }
+        }
+
+
+        // Use it
+        forAll(cellMap, i)
+        {
+            label cellI = cellMap[i];
+
+            if (cellI != -1 && numGapCells[cellI] != -1)
+            {
+                // Needed gap size
+                label cLevel = cellLevel[cellI];
+                scalar cellSize = meshCutter_.level0EdgeLength()/pow(2, cLevel);
+                scalar neededGapSize = numGapCells[cellI]*cellSize;
+
+                if (neededGapSize > detectedGapSize[cellI])
+                {
+                    if
+                    (
+                       !markForRefine
+                        (
+                            123,
+                            nAllowRefine,
+                            refineCell[cellI],
+                            nRefine
+                        )
+                    )
+                    {
+                        break;
+                    }
+                }
+            }
+        }
+
+
+        if
+        (
+            returnReduce(nRefine, sumOp<label>())
+          > returnReduce(nAllowRefine, sumOp<label>())
+        )
+        {
+            Info<< "Reached refinement limit." << endl;
+        }
+    }
+
+    return returnReduce(nRefine-oldNRefine, sumOp<label>());
+}
+
+
+Foam::label Foam::meshRefinement::markSmallFeatureRefinement
+(
+    const scalar planarCos,
+    const label nAllowRefine,
+    const labelList& neiLevel,
+    const pointField& neiCc,
+
+    labelList& refineCell,
+    label& nRefine
+) const
+{
+    const labelList& cellLevel = meshCutter_.cellLevel();
+    const labelList& surfaceIndices = surfaces_.surfaces();
+    const List<FixedList<label, 3> >& extendedGapLevel =
+        surfaces_.extendedGapLevel();
+    const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
+
+    label oldNRefine = nRefine;
+
+    // Check that we're using any gap refinement
+    labelList shellMaxLevel(shells_.maxGapLevel());
+
+    if (max(shellMaxLevel) == 0)
+    {
+        return 0;
+    }
+
+    //- Force calculation of tetBasePt
+    (void)mesh_.tetBasePtIs();
+    (void)mesh_.cellTree();
+
+
+    forAll(surfaceIndices, surfI)
+    {
+        label geomI = surfaceIndices[surfI];
+        const searchableSurface& geom = surfaces_.geometry()[geomI];
+
+
+        // Get the element index in a roundabout way. Problem is e.g.
+        // distributed surface where local indices differ from global
+        // ones (needed for getRegion call)
+
+        pointField ctrs;
+        labelList region;
+        vectorField normal;
+        {
+            // Representative local coordinates and bounding sphere
+            scalarField radiusSqr;
+            geom.boundingSpheres(ctrs, radiusSqr);
+
+            List<pointIndexHit> info;
+            geom.findNearest(ctrs, radiusSqr, info);
+
+            forAll(info, i)
+            {
+                if (!info[i].hit())
+                {
+                    FatalErrorIn("meshRefinement::markSmallFeatureRefinement")
+                        << "fc:" << ctrs[i]
+                        << " radius:" << radiusSqr[i]
+                        << exit(FatalError);
+                }
+            }
+
+            geom.getRegion(info, region);
+            geom.getNormal(info, normal);
+        }
+
+        // Do test to see whether triangles are inside/outside shell with
+        // applicable specification (minLevel <= celllevel < maxLevel)
+        List<FixedList<label, 3> > shellGapInfo;
+        List<volumeType> shellGapMode;
+        shells_.findHigherGapLevel
+        (
+            ctrs,
+            labelList(ctrs.size(), 0),
+            shellGapInfo,
+            shellGapMode
+        );
+
+
+        DynamicList<label> map(ctrs.size());
+        DynamicList<label> cellMap(ctrs.size());
+
+        DynamicField<point> rayStart(ctrs.size());
+        DynamicField<point> rayEnd(ctrs.size());
+        DynamicField<scalar> gapSize(ctrs.size());
+
+        label nTestCells = 0;
+
+        forAll(ctrs, i)
+        {
+            if (shellGapInfo[i][2] > 0)
+            {
+                label globalRegionI = surfaces_.globalRegion(surfI, region[i]);
+
+                // Combine info from shell and surface
+                FixedList<label, 3> gapInfo;
+                volumeType gapMode;
+                mergeGapInfo
+                (
+                    shellGapInfo[i],
+                    shellGapMode[i],
+                    extendedGapLevel[globalRegionI],
+                    extendedGapMode[globalRegionI],
+                    gapInfo,
+                    gapMode
+                );
+
+                //- Option 1: use octree nearest searching inside polyMesh
+                //label cellI = mesh_.findCell(pt);
+
+                //- Option 2: use octree 'inside' searching inside polyMesh. Is
+                //            much faster.
+                label cellI = -1;
+                const indexedOctree<treeDataCell>& tree = mesh_.cellTree();
+                if (tree.nodes().size() && tree.bb().contains(ctrs[i]))
+                {
+                    cellI = tree.findInside(ctrs[i]);
+                }
+
+                if (cellI != -1 && refineCell[cellI] == -1)
+                {
+                    // Construct one or two rays to test for oppositeness
+                    // Note that we always want to use the surface normal
+                    // and not the vector from cell centre to surface point
+
+                    label nRays = generateRays
+                    (
+                        ctrs[i],
+                        normal[i],
+                        gapInfo,
+                        gapMode,
+
+                        cellLevel[cellI],
+
+                        rayStart,
+                        rayEnd
+                    );
+
+                    if (nRays > 0)
+                    {
+                        nTestCells++;
+                        for (label j = 0; j < nRays; j++)
+                        {
+                            cellMap.append(cellI);
+                            map.append(i);
+                        }
+                    }
+                }
+            }
+        }
+
+        Info<< "Selected " << returnReduce(nTestCells, sumOp<label>())
+            << " cells containing triangle centres out of "
+            << mesh_.globalData().nTotalCells() << endl;
+        map.shrink();
+        cellMap.shrink();
+        rayStart.shrink();
+        rayEnd.shrink();
+
+        ctrs.clear();
+        region.clear();
+        shellGapInfo.clear();
+        shellGapMode.clear();
+        normal = UIndirectList<vector>(normal, map)();
+
+        // Do intersections.
+        labelList surfaceHit;
+        vectorField surfaceNormal;
+        surfaces_.findNearestIntersection
+        (
+            rayStart,
+            rayEnd,
+            surfaceHit,
+            surfaceNormal
+        );
+
+
+        label nOldRefine = 0;
+
+
+        forAll(surfaceHit, i)
+        {
+            if (surfaceHit[i] != -1) // && surf2[i] != -1)
+            {
+                // Found intersection with surface. Check coplanar normals.
+                label cellI = cellMap[i];
+
+                if (mag(normal[i]&surfaceNormal[i]) > planarCos)
+                {
+                    if
+                    (
+                       !markForRefine
+                        (
+                            surfaceHit[i],
+                            nAllowRefine,
+                            refineCell[cellI],
+                            nRefine
+                        )
+                    )
+                    {
+                        break;
+                    }
+                }
+            }
+        }
+
+        Info<< "For surface " << geom.name() << " found "
+            << returnReduce(nRefine-nOldRefine, sumOp<label>())
+            << " cells in small gaps" << endl;
+
+        if
+        (
+            returnReduce(nRefine, sumOp<label>())
+          > returnReduce(nAllowRefine, sumOp<label>())
+        )
+        {
+            Info<< "Reached refinement limit." << endl;
+        }
+    }
+
+    return returnReduce(nRefine-oldNRefine, sumOp<label>());
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index b27c4fae3f..4ffa1248f8 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/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)
-    {
-        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;
-    }
+    calcCellCellRays
+    (
+        neiCc,
+        neiLevel,
+        testFaces,
+        start,
+        end,
+        minLevel
+    );
 
 
     // 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)
-    {
-        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;
-    }
+    calcCellCellRays
+    (
+        neiCc,
+        neiLevel,
+        testFaces,
+        start,
+        end,
+        minLevel
+    );
 
 
     // Test for all intersections (with surfaces of higher gap level than
@@ -2045,7 +2052,10 @@ Foam::labelList Foam::meshRefinement::refineCandidates
     const bool internalRefinement,
     const bool surfaceRefinement,
     const bool curvatureRefinement,
+    const bool smallFeatureRefinement,
     const bool gapRefinement,
+    const bool bigGapRefinement,
+    const bool spreadGapSize,
     const label maxGlobalCells,
     const label maxLocalCells
 ) const
@@ -2096,6 +2106,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
         calcNeighbourData(neiLevel, neiCc);
 
 
+        const scalar planarCos = Foam::cos(degToRad(planarAngle));
+
 
         // Cells pierced by feature lines
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -2163,6 +2175,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 +2226,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 +2279,50 @@ 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
+            labelList numGapCells;
+            scalarField gapSize;
+            label nGap = markInternalGapRefinement
+            (
+                planarCos,
+                spreadGapSize,
+                nAllowRefine,
+
+                refineCell,
+                nRefine,
+                numGapCells,
+                gapSize
+            );
+            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
         // ~~~~~~~~~~~~~~~~~~~~
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/transportData.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/transportData.H
new file mode 100644
index 0000000000..eec088e6d2
--- /dev/null
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/transportData.H
@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::transportData
+
+Description
+    Holds information (coordinate and distance). Walks out 0.5*distance.
+
+SourceFiles
+    transportDataI.H
+    transportData.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef transportData_H
+#define transportData_H
+
+#include "wallPointData.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class transportData Declaration
+\*---------------------------------------------------------------------------*/
+
+class transportData
+:
+    public wallPointData<scalar>
+{
+public:
+
+    //- Class used to pass additional data in
+    class trackData
+    {
+    public:
+
+        //- Per face the index of the surface hit
+        const labelList& surfaceIndex_;
+
+        trackData(const labelList& surfaceIndex)
+        :
+            surfaceIndex_(surfaceIndex)
+        {}
+    };
+
+
+private:
+
+    // Private Member Functions
+
+        //- Evaluate distance to point. Update distSqr, origin from whomever
+        //  is nearer pt. Return true if w2 is closer to point,
+        //  false otherwise.
+        template<class TrackingData>
+        inline bool update
+        (
+            const point&,
+            const transportData& w2,
+            const scalar tol,
+            TrackingData& td
+        );
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline transportData();
+
+        //- Construct from origin, gapSize, distance
+        inline transportData
+        (
+            const point& origin,
+            const scalar gapSize,
+            const scalar distSqr
+        );
+
+
+    // Member Functions
+
+        // Needed by FaceCellWave
+
+            //- Influence of neighbouring face.
+            //  Calls update(...) with cellCentre of cellI
+            template<class TrackingData>
+            inline bool updateCell
+            (
+                const polyMesh& mesh,
+                const label thisCellI,
+                const label neighbourFaceI,
+                const transportData& neighbourWallInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of neighbouring cell.
+            //  Calls update(...) with faceCentre of faceI
+            template<class TrackingData>
+            inline bool updateFace
+            (
+                const polyMesh& mesh,
+                const label thisFaceI,
+                const label neighbourCellI,
+                const transportData& neighbourWallInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of different value on same face.
+            //  Merge new and old info.
+            //  Calls update(...) with faceCentre of faceI
+            template<class TrackingData>
+            inline bool updateFace
+            (
+                const polyMesh& mesh,
+                const label thisFaceI,
+                const transportData& neighbourWallInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+};
+
+
+//- Data associated with transportData type is same as underlying
+template<>
+inline bool contiguous<transportData>()
+{
+    return contiguous<wallPointData<scalar> >();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "transportDataI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/transportDataI.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/transportDataI.H
new file mode 100644
index 0000000000..e0bbeaeaf3
--- /dev/null
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/transportDataI.H
@@ -0,0 +1,180 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "polyMesh.H"
+#include "transform.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Update this with w2 if w2 nearer to pt.
+template<class TrackingData>
+inline bool Foam::transportData::update
+(
+    const point& pt,
+    const transportData& w2,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    scalar dist2 = magSqr(pt - w2.origin());
+
+    if (valid(td))
+    {
+        scalar diff = distSqr() - dist2;
+
+        if (diff < 0)
+        {
+            // Already nearer to pt
+            return false;
+        }
+
+        if ((diff < SMALL) || ((distSqr() > SMALL) && (diff/distSqr() < tol)))
+        {
+            // Don't propagate small changes
+            return false;
+        }
+    }
+
+    // Either *this is not yet valid or w2 is closer
+    {
+        // current not yet set so use any value
+        distSqr() = dist2;
+        origin() = w2.origin();
+        data() = w2.data();
+
+        if (distSqr() > sqr(0.25*data()))
+        {
+            // No need to transport gap data since too far away
+            return false;
+        }
+
+        return true;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::transportData::transportData()
+:
+    wallPointData<scalar>()
+{}
+
+
+inline Foam::transportData::transportData
+(
+    const point& origin,
+    const scalar gapSize,
+    const scalar distSqr
+)
+:
+    wallPointData<scalar>(origin, gapSize, distSqr)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class TrackingData>
+inline bool Foam::transportData::updateCell
+(
+    const polyMesh& mesh,
+    const label cellI,
+    const label faceI,
+    const transportData& neighbourWallInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    if (td.surfaceIndex_[faceI] != -1)
+    {
+        return false;
+    }
+
+    const vectorField& cellCentres = mesh.primitiveMesh::cellCentres();
+
+    bool updated = update
+    (
+        cellCentres[cellI],
+        neighbourWallInfo,
+        tol,
+        td
+    );
+
+    return updated;
+}
+
+
+template<class TrackingData>
+inline bool Foam::transportData::updateFace
+(
+    const polyMesh& mesh,
+    const label thisFaceI,
+    const label neighbourCellI,
+    const transportData& neighbourWallInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    if (td.surfaceIndex_[thisFaceI] != -1)
+    {
+        return false;
+    }
+
+    return update
+    (
+        mesh.faceCentres()[thisFaceI],
+        neighbourWallInfo,
+        tol,
+        td
+    );
+}
+
+
+template<class TrackingData>
+inline bool Foam::transportData::updateFace
+(
+    const polyMesh& mesh,
+    const label thisFaceI,
+    const transportData& neighbourWallInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    if (td.surfaceIndex_[thisFaceI] != -1)
+    {
+        return false;
+    }
+
+    return update
+    (
+        mesh.faceCentres()[thisFaceI],
+        neighbourWallInfo,
+        tol,
+        td
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index 7d47a612c9..b6fb43d98a 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/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
-(
-    const shellSurfaces& shells
-)
+void Foam::refinementSurfaces::setMinLevelFields(const shellSurfaces& shells)
 {
     forAll(surfaces_, surfI)
     {
@@ -1208,6 +1330,99 @@ 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::findNearestIntersection
+(
+    const pointField& start,
+    const pointField& end,
+
+    labelList& surface1,
+    List<pointIndexHit>& hitInfo1,
+    vectorField& normal1
+) const
+{
+    // 1. intersection from start to end
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Initialize arguments
+    surface1.setSize(start.size());
+    surface1 = -1;
+    hitInfo1.setSize(start.size());
+    hitInfo1 = pointIndexHit();
+    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;
+                hitInfo1[pointI] = nearestInfo[pointI];
+                normal1[pointI] = normal[pointI];
+                nearest[pointI] = nearestInfo[pointI].hitPoint();
+            }
+        }
+    }
+}
+
+
 void Foam::refinementSurfaces::findAnyIntersection
 (
     const pointField& start,
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
index fadb9aaa26..7e2bcdf946 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/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,25 @@ 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;
+
+            //- Find nearest (to start only) intersection of edge
+            void findNearestIntersection
+            (
+                const pointField& start,
+                const pointField& end,
+                labelList& surfaces,
+                List<pointIndexHit>&,
+                vectorField& normal
+            ) const;
+
             //- Used for debugging only: find intersection of edge.
             void findAnyIntersection
             (
diff --git a/src/mesh/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C b/src/mesh/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C
index cd39802553..128d2892ec 100644
--- a/src/mesh/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C
+++ b/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);
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.H b/src/mesh/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.H
index 9b0e4c9e05..13817de16f 100644
--- a/src/mesh/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.H
+++ b/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
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     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
-            //{
-            //    return distances_;
-            //}
-            //
-            ////- Per shell per distance the refinement level
-            //const labelListList& levels() const
-            //{
-            //    return levels_;
-            //}
+            //- Indices of surfaces that are shells
+            const labelList& shells() const
+            {
+                return shells_;
+            }
 
 
         // 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;
 };
 
 
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index 0c1b794f7c..23b11884e2 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -63,7 +63,9 @@ indexedOctree/treeDataTriSurface.C
 
 searchableSurface = searchableSurface
 $(searchableSurface)/searchableBox.C
+$(searchableSurface)/searchableRotatedBox.C
 $(searchableSurface)/searchableCylinder.C
+$(searchableSurface)/searchableCone.C
 $(searchableSurface)/searchableDisk.C
 $(searchableSurface)/searchablePlane.C
 $(searchableSurface)/searchablePlate.C
@@ -75,6 +77,7 @@ $(searchableSurface)/searchableSurfacesQueries.C
 $(searchableSurface)/searchableSurfaceWithGaps.C
 $(searchableSurface)/triSurfaceMesh.C
 $(searchableSurface)/closedTriSurfaceMesh.C
+$(searchableSurface)/subTriSurfaceMesh.C
 
 topoSets = sets/topoSets
 $(topoSets)/cellSet.C
diff --git a/src/meshTools/searchableSurface/searchableCone.C b/src/meshTools/searchableSurface/searchableCone.C
new file mode 100644
index 0000000000..b4c212a92b
--- /dev/null
+++ b/src/meshTools/searchableSurface/searchableCone.C
@@ -0,0 +1,1146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "searchableCone.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(searchableCone, 0);
+addToRunTimeSelectionTable(searchableSurface, searchableCone, dict);
+
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::pointField> Foam::searchableCone::coordinates() const
+{
+    tmp<pointField> tCtrs(new pointField(1, 0.5*(point1_ + point2_)));
+
+    return tCtrs;
+}
+
+
+void Foam::searchableCone::boundingSpheres
+(
+    pointField& centres,
+    scalarField& radiusSqr
+) const
+{
+    centres.setSize(1);
+    centres[0] = 0.5*(point1_ + point2_);
+
+    radiusSqr.setSize(1);
+    if (radius1_ > radius2_)
+    {
+        radiusSqr[0] = Foam::magSqr(point1_-centres[0]) + Foam::sqr(radius1_);
+    }
+    else
+    {
+        radiusSqr[0] = Foam::magSqr(point2_-centres[0]) + Foam::sqr(radius2_);
+    }
+
+    // Add a bit to make sure all points are tested inside
+    radiusSqr += Foam::sqr(SMALL);
+}
+
+
+Foam::tmp<Foam::pointField> Foam::searchableCone::points() const
+{
+    tmp<pointField> tPts(new pointField(2));
+    pointField& pts = tPts();
+
+    pts[0] = point1_;
+    pts[1] = point2_;
+
+    return tPts;
+}
+
+
+void Foam::searchableCone::findNearestAndNormal
+(
+    const point& sample,
+    const scalar nearestDistSqr,
+    pointIndexHit& info,
+    vector& nearNormal
+) const
+{
+    vector v(sample - point1_);
+
+    // Decompose sample-point1 into normal and parallel component
+    scalar parallel = (v & unitDir_);
+
+    // Remove the parallel component and normalise
+    v -= parallel*unitDir_;
+
+    scalar magV = mag(v);
+    if (magV < ROOTVSMALL)
+    {
+       v = vector::zero;
+    }
+    else
+    {
+       v /= magV;
+    }
+
+    // Nearest and normal on disk at point1
+    point disk1Point(point1_ + min(max(magV, innerRadius1_), radius1_)*v);
+    vector disk1Normal(-unitDir_);
+
+    // Nearest and normal on disk at point2
+    point disk2Point(point2_ + min(max(magV, innerRadius2_), radius2_)*v);
+    vector disk2Normal(unitDir_);
+
+    // Nearest and normal on cone. Initialise to far-away point so if not
+    // set it picks one of the disk points
+    point nearCone(GREAT, GREAT, GREAT);
+    vector normalCone(1, 0, 0);
+
+    // Nearest and normal on inner cone. Initialise to far-away point
+    point iCnearCone(GREAT, GREAT, GREAT);
+    vector iCnormalCone(1, 0, 0);
+
+    point projPt1 = point1_+ radius1_*v;
+    point projPt2 = point2_+ radius2_*v;
+
+    vector p1 = (projPt1 - point1_);
+    if (mag(p1) > ROOTVSMALL)
+    {
+        p1 /= mag(p1);
+
+        // Find vector along the two end of cone
+        vector b(projPt2 - projPt1);
+        scalar magS = mag(b);
+        b /= magS;
+
+        // Find the vector along sample pt and pt at one end of conde
+        vector a(sample - projPt1);
+
+        if (mag(a) <= ROOTVSMALL)
+        {
+            // Exception: sample on disk1. Redo with projPt2.
+            vector a(sample - projPt2);
+            // Find normal unitvector
+            nearCone = (a & b)*b+projPt2;
+            vector b1 = (p1 & b)*b;
+            normalCone = p1 - b1;
+            normalCone /= mag(normalCone);
+        }
+        else
+        {
+            // Find neartest point on cone surface
+            nearCone = (a & b)*b+projPt1;
+            // Find projection along surface of cone
+            vector b1 = (p1 & b)*b;
+            normalCone = p1 - b1;
+            normalCone /= mag(normalCone);
+        }
+
+        if (innerRadius1_ > 0 || innerRadius2_ > 0)
+        {
+            // Same for inner radius
+            point iCprojPt1 = point1_+ innerRadius1_*v;
+            point iCprojPt2 = point2_+ innerRadius2_*v;
+
+            vector iCp1 = (iCprojPt1 - point1_);
+            iCp1 /= mag(iCp1);
+
+            // Find vector along the two end of cone
+            vector iCb(iCprojPt2 - iCprojPt1);
+            magS = mag(iCb);
+            iCb /= magS;
+
+
+            // Find the vector along sample pt and pt at one end of conde
+            vector iCa(sample - iCprojPt1);
+
+            if (mag(iCa) <= ROOTVSMALL)
+            {
+                vector iCa(sample - iCprojPt2);
+
+                // Find normal unitvector
+                iCnearCone = (iCa & iCb)*iCb+iCprojPt2;
+                vector b1 = (iCp1 & iCb)*iCb;
+                iCnormalCone = iCp1 - b1;
+                iCnormalCone /= mag(iCnormalCone);
+            }
+            else
+            {
+                // Find nearest point on cone surface
+                iCnearCone = (iCa & iCb)*iCb+iCprojPt1;
+                // Find projection along surface of cone
+                vector b1 = (iCp1 & iCb)*iCb;
+                iCnormalCone = iCp1 - b1;
+                iCnormalCone /= mag(iCnormalCone);
+            }
+        }
+    }
+
+
+    // Select nearest out of the 4 points (outer cone, disk1, disk2, inner
+    // cone)
+
+    FixedList<scalar, 4> dist;
+    dist[0] = magSqr(nearCone-sample);
+    dist[1] = magSqr(disk1Point-sample);
+    dist[2] = magSqr(disk2Point-sample);
+    dist[3] = magSqr(iCnearCone-sample);
+
+    label minI = findMin(dist);
+
+
+    // Snap the point to the corresponding surface
+
+    if (minI == 0)  // Outer cone
+    {
+        // Closest to (infinite) outer cone. See if needs clipping to end disks
+
+        {
+            vector v1(nearCone-point1_);
+            scalar para = (v1 & unitDir_);
+            // Remove the parallel component and normalise
+            v1 -= para*unitDir_;
+            scalar magV1 = mag(v1);
+            v1 = v1/magV1;
+            if (para < 0.0 && magV1 >= radius1_)
+            {
+                // Near point 1. Set point to intersection of disk and cone.
+                // Keep normal from cone.
+                nearCone = disk1Point;
+            }
+            else if (para < 0.0 && magV1 < radius1_)
+            {
+                // On disk1
+                nearCone = disk1Point;
+                normalCone = disk1Normal;
+            }
+            else if (para > magDir_ && magV1 >= radius2_)
+            {
+                // Near point 2. Set point to intersection of disk and cone.
+                // Keep normal from cone.
+                nearCone = disk2Point;
+            }
+            else if (para > magDir_ && magV1 < radius2_)
+            {
+                // On disk2
+                nearCone = disk2Point;
+                normalCone = disk2Normal;
+            }
+        }
+        info.setPoint(nearCone);
+        nearNormal = normalCone;
+    }
+    else if (minI == 1) // Near to disk1
+    {
+        info.setPoint(disk1Point);
+        nearNormal = disk1Normal;
+    }
+    else if (minI == 2) // Near to disk2
+    {
+        info.setPoint(disk2Point);
+        nearNormal = disk2Normal;
+    }
+    else                // Near to inner cone
+    {
+        {
+            vector v1(iCnearCone-point1_);
+            scalar para = (v1 & unitDir_);
+            // Remove the parallel component and normalise
+            v1 -= para*unitDir_;
+            scalar magV1 = mag(v1);
+            v1 = v1/magV1;
+
+            if (para < 0.0 && magV1 >= innerRadius1_)
+            {
+               iCnearCone = disk1Point;
+            }
+            else if (para < 0.0 && magV1 < innerRadius1_)
+            {
+               iCnearCone = disk1Point;
+               iCnormalCone = disk1Normal;
+            }
+            else if (para > magDir_ && magV1 >= innerRadius2_)
+            {
+               iCnearCone = disk2Point;
+            }
+            else if (para > magDir_ && magV1 < innerRadius2_)
+            {
+               iCnearCone = disk2Point;
+               iCnormalCone = disk2Normal;
+            }
+        }
+        info.setPoint(iCnearCone);
+        nearNormal = iCnormalCone;
+    }
+
+
+    if (magSqr(sample - info.rawPoint()) < nearestDistSqr)
+    {
+        info.setHit();
+        info.setIndex(0);
+    }
+}
+
+
+Foam::scalar Foam::searchableCone::radius2
+(
+    const searchableCone& cone,
+    const point& pt
+)
+{
+    const vector x = (pt-cone.point1_) ^ cone.unitDir_;
+    return x&x;
+}
+
+
+// From mrl.nyu.edu/~dzorin/rend05/lecture2.pdf,
+// Ray Tracing II, Infinite cone ray intersection.
+void Foam::searchableCone::findLineAll
+(
+    const searchableCone& cone,
+    const scalar innerRadius1,
+    const scalar innerRadius2,
+    const point& start,
+    const point& end,
+    pointIndexHit& near,
+    pointIndexHit& far
+) const
+{
+    near.setMiss();
+    far.setMiss();
+
+    vector point1Start(start-cone.point1_);
+    vector point2Start(start-cone.point2_);
+    vector point1End(end-cone.point1_);
+
+    // Quick rejection of complete vector outside endcaps
+    scalar s1 = point1Start&(cone.unitDir_);
+    scalar s2 = point1End&(cone.unitDir_);
+
+    if ((s1 < 0.0 && s2 < 0.0) || (s1 > cone.magDir_ && s2 > cone.magDir_))
+    {
+        return;
+    }
+
+    // Line as P = start+t*V where V is unit vector and t=[0..mag(end-start)]
+    vector V(end-start);
+    scalar magV = mag(V);
+    if (magV < ROOTVSMALL)
+    {
+        return;
+    }
+    V /= magV;
+
+
+    // We now get the nearest intersections to start. This can either be
+    // the intersection with the end plane or with the cylinder side.
+
+    // Get the two points (expressed in t) on the end planes. This is to
+    // clip any cylinder intersection against.
+    scalar tPoint1;
+    scalar tPoint2;
+
+    // Maintain the two intersections with the endcaps
+    scalar tNear = VGREAT;
+    scalar tFar = VGREAT;
+
+    scalar radius_sec = cone.radius1_;
+
+    {
+        // Find dot product: mag(s)>VSMALL suggest that it is greater
+        scalar s = (V&unitDir_);
+        if (mag(s) > VSMALL)
+        {
+            tPoint1 = -s1/s;
+            tPoint2 = -(point2Start&(cone.unitDir_))/s;
+
+            if (tPoint2 < tPoint1)
+            {
+                Swap(tPoint1, tPoint2);
+            }
+            if (tPoint1 > magV || tPoint2 < 0)
+            {
+                return;
+            }
+            // See if the points on the endcaps are actually inside the cylinder
+            if (tPoint1 >= 0 && tPoint1 <= magV)
+            {
+                scalar r2 = radius2(cone, start+tPoint1*V);
+                vector p1 = (start+tPoint1*V-point1_);
+                vector p2 = (start+tPoint1*V-point2_);
+                radius_sec = cone.radius1_;
+                scalar inC_radius_sec = innerRadius1_;
+
+                if (mag(p2&(cone.unitDir_)) < mag(p1&(cone.unitDir_)))
+                {
+                    radius_sec = cone.radius2_;
+                    inC_radius_sec = innerRadius2_;
+                }
+
+                if (r2 <= sqr(radius_sec) && r2 >= sqr(inC_radius_sec))
+                {
+                    tNear = tPoint1;
+                }
+            }
+            if (tPoint2 >= 0 && tPoint2 <= magV)
+            {
+                // Decompose sample-point1 into normal and parallel component
+                vector p1 = (start+tPoint2*V-cone.point1_);
+                vector p2 = (start+tPoint2*V-cone.point2_);
+                radius_sec = cone.radius1_;
+                scalar inC_radius_sec = innerRadius1_;
+                if (mag(p2&(cone.unitDir_)) < mag(p1&(cone.unitDir_)))
+                {
+                    radius_sec = cone.radius2_;
+                    inC_radius_sec = innerRadius2_;
+                }
+                scalar r2 = radius2(cone, start+tPoint2*V);
+                if (r2 <= sqr(radius_sec) && r2 >= sqr(inC_radius_sec))
+                {
+                    // Check if already have a near hit from point1
+                    if (tNear <= 1)
+                    {
+                        tFar = tPoint2;
+                    }
+                    else
+                    {
+                        tNear = tPoint2;
+                    }
+                }
+            }
+        }
+        else
+        {
+            // Vector perpendicular to cylinder. Check for outside already done
+            // above so just set tpoint to allow all.
+            tPoint1 = -VGREAT;
+            tPoint2 = VGREAT;
+        }
+    }
+
+
+    // Second order equation of the form a*t^2 + b*t + c
+    scalar a, b, c;
+
+    scalar deltaRadius = cone.radius2_-cone.radius1_;
+    if (mag(deltaRadius) <= ROOTVSMALL)
+    {
+        vector point1Start(start-cone.point1_);
+        const vector x = point1Start ^ cone.unitDir_;
+        const vector y = V ^ cone.unitDir_;
+        const scalar d = sqr(0.5*(cone.radius1_ + cone.radius2_));
+
+        a = (y&y);
+        b = 2*(x&y);
+        c = (x&x)-d;
+    }
+    else
+    {
+        vector va = cone.unitDir_;
+        vector v1(end-start);
+        v1 = v1/mag(v1);
+        scalar p  = (va&v1);
+        vector a1 = (v1-p*va);
+
+        // Determine the end point of the cone
+        point pa =
+            cone.unitDir_*cone.radius1_*mag(cone.point2_-cone.point1_)
+           /(-deltaRadius)
+           +cone.point1_;
+
+        scalar l2 = sqr(deltaRadius)+sqr(cone.magDir_);
+        scalar sqrCosAlpha = sqr(cone.magDir_)/l2;
+        scalar sqrSinAlpha = sqr(deltaRadius)/l2;
+
+
+        vector delP(start-pa);
+        vector p1 = (delP-(delP&va)*va);
+
+        a = sqrCosAlpha*((v1-p*va)&(v1-p*va))-sqrSinAlpha*p*p;
+        b =
+            2.0*sqrCosAlpha*(a1&p1)
+           -2.0*sqrSinAlpha*(v1&va)*(delP&va);
+        c =
+            sqrCosAlpha
+           *((delP-(delP&va)*va)&(delP-(delP&va)*va))
+           -sqrSinAlpha*sqr(delP&va);
+    }
+
+    const scalar disc = b*b-4.0*a*c;
+
+    scalar t1 = -VGREAT;
+    scalar t2 = VGREAT;
+
+    if (disc < 0)
+    {
+        // Fully outside
+        return;
+    }
+    else if (disc < ROOTVSMALL)
+    {
+        // Single solution
+        if (mag(a) > ROOTVSMALL)
+        {
+            t1 = -b/(2.0*a);
+
+            if (t1 >= 0.0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
+            {
+                // Valid. Insert sorted.
+                if (t1 < tNear)
+                {
+                    tFar = tNear;
+                    tNear = t1;
+                }
+                else if (t1 < tFar)
+                {
+                    tFar = t1;
+                }
+            }
+            else
+            {
+                return;
+            }
+        }
+        else
+        {
+            // Aligned with axis. Check if outside radius
+            if (c > 0.0)
+            {
+                return;
+            }
+        }
+    }
+    else
+    {
+        if (mag(a) > ROOTVSMALL)
+        {
+            scalar sqrtDisc = sqrt(disc);
+
+            t1 = (-b - sqrtDisc)/(2.0*a);
+            t2 = (-b + sqrtDisc)/(2.0*a);
+            if (t2 < t1)
+            {
+                Swap(t1, t2);
+            }
+
+            if (t1 >= 0.0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
+            {
+                // Valid. Insert sorted.
+                if (t1 < tNear)
+                {
+                    tFar = tNear;
+                    tNear = t1;
+                }
+                else if (t1 < tFar)
+                {
+                    tFar = t1;
+                }
+            }
+            if (t2>=0 && t2 <= magV && t2 >= tPoint1 && t2 <= tPoint2)
+            {
+                // Valid. Insert sorted.
+                if (t2 < tNear)
+                {
+                    tFar = tNear;
+                    tNear = t2;
+                }
+                else if (t2 < tFar)
+                {
+                    tFar = t2;
+                }
+            }
+        }
+        else
+        {
+            // Aligned with axis. Check if outside radius
+            if (c > 0.0)
+            {
+                return;
+            }
+        }
+    }
+
+    // Check tNear, tFar
+    if (tNear>=0 && tNear <= magV)
+    {
+        near.setPoint(start+tNear*V);
+        near.setHit();
+        near.setIndex(0);
+        if (tFar <= magV)
+        {
+            far.setPoint(start+tFar*V);
+            far.setHit();
+            far.setIndex(0);
+        }
+    }
+    else if (tFar>=0 && tFar <= magV)
+    {
+        near.setPoint(start+tFar*V);
+        near.setHit();
+        near.setIndex(0);
+    }
+}
+
+
+void Foam::searchableCone::insertHit
+(
+    const point& start,
+    const point& end,
+    List<pointIndexHit>& info,
+    const pointIndexHit& hit
+) const
+{
+    scalar smallDistSqr = SMALL*magSqr(end-start);
+
+    scalar hitMagSqr = magSqr(hit.hitPoint()-start);
+
+    forAll(info, i)
+    {
+        scalar d2 = magSqr(info[i].hitPoint()-start);
+
+        if (d2 > hitMagSqr+smallDistSqr)
+        {
+            // Insert at i.
+            label sz = info.size();
+            info.setSize(sz+1);
+            for (label j = sz; j > i; --j)
+            {
+                info[j] = info[j-1];
+            }
+            info[i] = hit;
+            return;
+        }
+        else if (d2 > hitMagSqr-smallDistSqr)
+        {
+            // hit is same point as info[i].
+            return;
+        }
+    }
+    // Append
+    label sz = info.size();
+    info.setSize(sz+1);
+    info[sz] = hit;
+}
+
+
+Foam::boundBox Foam::searchableCone::calcBounds() const
+{
+    // Adapted from
+    // http://www.gamedev.net/community/forums
+    //       /topic.asp?topic_id=338522&forum_id=20&gforum_id=0
+
+    // Let cylinder have end points A,B and radius r,
+
+    // Bounds in direction X (same for Y and Z) can be found as:
+    // Let A.X<B.X (otherwise swap points)
+    // Good approximate lowest bound is A.X-r and highest is B.X+r (precise for
+    // capsule). At worst, in one direction it can be larger than needed by 2*r.
+
+    // Accurate bounds for cylinder is
+    // A.X-kx*r, B.X+kx*r
+    // where
+    // kx=sqrt(((A.Y-B.Y)^2+(A.Z-B.Z)^2)/((A.X-B.X)^2+(A.Y-B.Y)^2+(A.Z-B.Z)^2))
+
+    // similar thing for Y and Z
+    // (i.e.
+    // ky=sqrt(((A.X-B.X)^2+(A.Z-B.Z)^2)/((A.X-B.X)^2+(A.Y-B.Y)^2+(A.Z-B.Z)^2))
+    // kz=sqrt(((A.X-B.X)^2+(A.Y-B.Y)^2)/((A.X-B.X)^2+(A.Y-B.Y)^2+(A.Z-B.Z)^2))
+    // )
+
+    // How derived: geometric reasoning. Bounds of cylinder is same as for 2
+    // circles centered on A and B. This sqrt thingy gives sine of angle between
+    // axis and direction, used to find projection of radius.
+
+    vector kr
+    (
+        sqrt(sqr(unitDir_.y()) + sqr(unitDir_.z())),
+        sqrt(sqr(unitDir_.x()) + sqr(unitDir_.z())),
+        sqrt(sqr(unitDir_.x()) + sqr(unitDir_.y()))
+    );
+
+    if (radius2_ >= radius1_)
+    {
+        kr *= radius2_;
+    }
+    else
+    {
+        kr *= radius1_;
+    }
+
+    point min = point1_ - kr;
+    point max = point1_ + kr;
+
+    min = ::Foam::min(min, point2_ - kr);
+    max = ::Foam::max(max, point2_ + kr);
+
+    return boundBox(min, max);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::searchableCone::searchableCone
+(
+    const IOobject& io,
+    const point& point1,
+    const scalar radius1,
+    const scalar innerRadius1,
+    const point& point2,
+    const scalar radius2,
+    const scalar innerRadius2
+)
+:
+    searchableSurface(io),
+    point1_(point1),
+    radius1_(radius1),
+    innerRadius1_(innerRadius1),
+    point2_(point2),
+    radius2_(radius2),
+    innerRadius2_(innerRadius2),
+    magDir_(mag(point2_-point1_)),
+    unitDir_((point2_-point1_)/magDir_)
+{
+    bounds() = calcBounds();
+}
+
+
+Foam::searchableCone::searchableCone
+(
+    const IOobject& io,
+    const dictionary& dict
+)
+:
+    searchableSurface(io),
+    point1_(dict.lookup("point1")),
+    radius1_(readScalar(dict.lookup("radius1"))),
+    innerRadius1_(dict.lookupOrDefault("innerRadius1", 0.0)),
+    point2_(dict.lookup("point2")),
+    radius2_(readScalar(dict.lookup("radius2"))),
+    innerRadius2_(dict.lookupOrDefault("innerRadius2", 0.0)),
+    magDir_(mag(point2_-point1_)),
+    unitDir_((point2_-point1_)/magDir_)
+{
+    bounds() = calcBounds();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::searchableCone::~searchableCone()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::wordList& Foam::searchableCone::regions() const
+{
+    if (regions_.empty())
+    {
+        regions_.setSize(1);
+        regions_[0] = "region0";
+    }
+    return regions_;
+}
+
+
+void Foam::searchableCone::findNearest
+(
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+    List<pointIndexHit>& info
+) const
+{
+    info.setSize(samples.size());
+    forAll(samples, i)
+    {
+        vector normal;
+        findNearestAndNormal(samples[i], nearestDistSqr[i], info[i], normal);
+    }
+}
+
+
+void Foam::searchableCone::findLine
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    info.setSize(start.size());
+
+    forAll(start, i)
+    {
+        // Pick nearest intersection. If none intersected take second one.
+        pointIndexHit b;
+        findLineAll
+        (
+            *this,
+            innerRadius1_,
+            innerRadius2_,
+            start[i],
+            end[i],
+            info[i],
+            b
+        );
+        if (!info[i].hit() && b.hit())
+        {
+            info[i] = b;
+        }
+    }
+
+
+    if (innerRadius1_ > 0.0 || innerRadius2_ > 0.0)
+    {
+        IOobject io(*this);
+        io.rename(name()+"Inner");
+        searchableCone innerCone
+        (
+            io,
+            point1_,
+            innerRadius1_,
+            0.0,
+            point2_,
+            innerRadius2_,
+            0.0
+        );
+
+        forAll(info, i)
+        {
+            point newEnd;
+            if (info[i].hit())
+            {
+                newEnd = info[i].hitPoint();
+            }
+            else
+            {
+                newEnd = end[i];
+            }
+            pointIndexHit near;
+            pointIndexHit far;
+            findLineAll
+            (
+                innerCone,
+                innerRadius1_,
+                innerRadius2_,
+                start[i],
+                newEnd,
+                near,
+                far
+            );
+
+            if (near.hit())
+            {
+                info[i] = near;
+            }
+            else if (far.hit())
+            {
+                info[i] = far;
+            }
+        }
+    }
+}
+
+
+void Foam::searchableCone::findLineAny
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    info.setSize(start.size());
+
+    forAll(start, i)
+    {
+        // Discard far intersection
+        pointIndexHit b;
+        findLineAll
+        (
+            *this,
+            innerRadius1_,
+            innerRadius2_,
+            start[i],
+            end[i],
+            info[i],
+            b
+        );
+        if (!info[i].hit() && b.hit())
+        {
+            info[i] = b;
+        }
+    }
+    if (innerRadius1_ > 0.0 || innerRadius2_ > 0.0)
+    {
+        IOobject io(*this);
+        io.rename(name()+"Inner");
+        searchableCone cone
+        (
+            io,
+            point1_,
+            innerRadius1_,
+            0.0,
+            point2_,
+            innerRadius2_,
+            0.0
+        );
+
+        forAll(info, i)
+        {
+            if (!info[i].hit())
+            {
+                pointIndexHit b;
+                findLineAll
+                (
+                    cone,
+                    innerRadius1_,
+                    innerRadius2_,
+                    start[i],
+                    end[i],
+                    info[i],
+                    b
+                );
+                if (!info[i].hit() && b.hit())
+                {
+                    info[i] = b;
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::searchableCone::findLineAll
+(
+    const pointField& start,
+    const pointField& end,
+    List<List<pointIndexHit> >& info
+) const
+{
+    info.setSize(start.size());
+
+    forAll(start, i)
+    {
+        pointIndexHit near, far;
+        findLineAll
+        (
+            *this,
+            innerRadius1_,
+            innerRadius2_,
+            start[i],
+            end[i],
+            near,
+            far
+        );
+
+        if (near.hit())
+        {
+            if (far.hit())
+            {
+                info[i].setSize(2);
+                info[i][0] = near;
+                info[i][1] = far;
+            }
+            else
+            {
+                info[i].setSize(1);
+                info[i][0] = near;
+            }
+        }
+        else
+        {
+            if (far.hit())
+            {
+                info[i].setSize(1);
+                info[i][0] = far;
+            }
+            else
+            {
+                info[i].clear();
+            }
+        }
+    }
+
+    if (innerRadius1_ > 0.0 || innerRadius2_ > 0.0)
+    {
+        IOobject io(*this);
+        io.rename(name()+"Inner");
+        searchableCone cone
+        (
+            io,
+            point1_,
+            innerRadius1_,
+            0.0,
+            point2_,
+            innerRadius2_,
+            0.0
+        );
+
+        forAll(info, i)
+        {
+            pointIndexHit near;
+            pointIndexHit far;
+            findLineAll
+            (
+                cone,
+                innerRadius1_,
+                innerRadius2_,
+                start[i],
+                end[i],
+                near,
+                far
+            );
+
+            if (near.hit())
+            {
+                insertHit(start[i], end[i], info[i], near);
+            }
+            if (far.hit())
+            {
+                insertHit(start[i], end[i], info[i], far);
+            }
+        }
+    }
+}
+
+
+void Foam::searchableCone::getRegion
+(
+    const List<pointIndexHit>& info,
+    labelList& region
+) const
+{
+    region.setSize(info.size());
+    region = 0;
+}
+
+
+void Foam::searchableCone::getNormal
+(
+    const List<pointIndexHit>& info,
+    vectorField& normal
+) const
+{
+    normal.setSize(info.size());
+    normal = vector::zero;
+
+    forAll(info, i)
+    {
+        if (info[i].hit())
+        {
+            pointIndexHit nearInfo;
+            findNearestAndNormal
+            (
+                info[i].hitPoint(),
+                Foam::sqr(GREAT),
+                nearInfo,
+                normal[i]
+            );
+        }
+    }
+}
+
+
+void Foam::searchableCone::getVolumeType
+(
+    const pointField& points,
+    List<volumeType>& volType
+) const
+{
+    volType.setSize(points.size());
+    volType = volumeType::INSIDE;
+
+    scalar magDir = magDir_;
+
+    bool swap = false;
+    bool swapInner = false;
+
+    if (radius1_ > radius2_)
+    {
+         swap = true;
+    }
+
+    if (innerRadius1_ > innerRadius2_)
+    {
+         swapInner = true;
+    }
+
+    forAll(points, pointI)
+    {
+        const point& pt = points[pointI];
+
+        vector v(pt - point1_);
+
+        // Decompose sample-point1 into normal and parallel component
+        scalar parallel = v & unitDir_;
+        scalar comp = parallel;
+        scalar compInner = parallel;
+        if (swap)
+        {
+            comp = -(magDir-parallel);
+        }
+
+        if (swapInner)
+        {
+            compInner = -(magDir-parallel);
+        }
+
+        scalar radius_sec = radius1_+comp*(radius2_-radius1_)/magDir_;
+
+        scalar radius_sec_inner =
+            innerRadius1_
+           +compInner*(innerRadius2_-innerRadius1_)/magDir_;
+
+        if (parallel < 0)
+        {
+            // Left of point1 endcap
+            volType[pointI] = volumeType::OUTSIDE;
+        }
+        else if (parallel > magDir_)
+        {
+            // Right of point2 endcap
+            volType[pointI] = volumeType::OUTSIDE;
+        }
+        else
+        {
+            // Remove the parallel component
+            v -= parallel*unitDir_;
+            if (mag(v) >= radius_sec_inner && mag(v) <= radius_sec)
+            {
+                volType[pointI] = volumeType::INSIDE;
+            }
+            else
+            {
+                volType[pointI] = volumeType::OUTSIDE;
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/searchableCone.H b/src/meshTools/searchableSurface/searchableCone.H
new file mode 100644
index 0000000000..3b4d7c3d1f
--- /dev/null
+++ b/src/meshTools/searchableSurface/searchableCone.H
@@ -0,0 +1,296 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::searchableCone
+
+Description
+    Searching on (optionally hollow) cone.
+
+    \heading Function object usage
+    \table
+        Property     | Description             | Required    | Default value
+        point1       | coordinate of endpoint  | yes         |
+        radius1      | radius at point1        | yes         | yes
+        innerRadius1 | inner radius at point1  | no          |
+        point2       | coordinate of endpoint  | yes         |
+        radius2      | radius at point2        | yes         | yes
+        innerRadius2 | inner radius at point2  | no          |
+    \endtable
+
+Note
+    Initial implementation, might suffer from robustness (e.g. radius1==radius2)
+
+SourceFiles
+    searchableCone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchableCone_H
+#define searchableCone_H
+
+#include "searchableSurface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class searchableCone Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchableCone
+:
+    public searchableSurface
+{
+    // Private Member Data
+
+        //- 'Left' point
+        const point point1_;
+
+        //- Outer radius at point1
+        const scalar radius1_;
+
+        //- Inner radius at point1
+        const scalar innerRadius1_;
+
+
+        //- 'Right' point
+        const point point2_;
+
+        //- Outer radius at point2
+        const scalar radius2_;
+
+        //- Inner radius at point2
+        const scalar innerRadius2_;
+
+
+        //- Length of vector point2-point1
+        const scalar magDir_;
+
+        //- Normalised vector point2-point1
+        const vector unitDir_;
+
+
+        //- Names of regions
+        mutable wordList regions_;
+
+
+   // Private Member Functions
+
+        //- Find nearest point on cylinder.
+        void findNearestAndNormal
+        (
+            const point& sample,
+            const scalar nearestDistSqr,
+            pointIndexHit & nearInfo,
+            vector& normal
+        ) const;
+
+        //- Determine radial coordinate (squared)
+        static scalar radius2(const searchableCone& cone, const point& pt);
+
+        //- Find both intersections with cone. innerRadii supplied externally.
+        void findLineAll
+        (
+            const searchableCone& cone,
+            const scalar innerRadius1,
+            const scalar innerRadius2,
+            const point& start,
+            const point& end,
+            pointIndexHit& near,
+            pointIndexHit& far
+        ) const;
+
+        //- Insert a hit if it differs (by a tolerance) from the existing ones
+        void insertHit
+        (
+            const point& start,
+            const point& end,
+            List<pointIndexHit>& info,
+            const pointIndexHit& hit
+        ) const;
+
+        //- Return the boundBox of the cylinder
+        boundBox calcBounds() const;
+
+        //- Disallow default bitwise copy construct
+        searchableCone(const searchableCone&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const searchableCone&);
+
+
+public:
+
+    //- Runtime type information
+   TypeName("searchableCone");
+
+
+    // Constructors
+
+        //- Construct from components
+       searchableCone
+       (
+           const IOobject& io,
+           const point& point1,
+           const scalar radius1,
+           const scalar innerRadius1,
+           const point& point2,
+           const scalar radius2,
+           const scalar innerRadius2
+       );
+
+        //- Construct from dictionary (used by searchableSurface)
+        searchableCone
+        (
+            const IOobject& io,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+
+        virtual ~searchableCone();
+
+
+    // Member Functions
+
+        virtual const wordList& regions() const;
+
+        //- Whether supports volume type below
+        virtual bool hasVolumeType() const
+        {
+            return true;
+        }
+
+        //- Range of local indices that can be returned.
+        virtual label size() const
+        {
+            return 1;
+        }
+
+        //- Get representative set of element coordinates
+        //  Usually the element centres (should be of length size()).
+        virtual tmp<pointField> coordinates() const;
+
+        //- Get bounding spheres (centre and radius squared), one per element.
+        //  Any point on element is guaranteed to be inside.
+        virtual void boundingSpheres
+        (
+            pointField& centres,
+            scalarField& radiusSqr
+        ) const;
+
+        //- Get the points that define the surface.
+        virtual tmp<pointField> points() const;
+
+        //- Does any part of the surface overlap the supplied bound box?
+        virtual bool overlaps(const boundBox& bb) const
+        {
+            notImplemented
+            (
+                "searchableCone::overlaps(const boundBox&) const"
+            );
+
+            return false;
+        }
+
+
+        // Multiple point queries.
+
+            //- Find nearest point on cylinder
+            virtual void findNearest
+            (
+                const pointField& sample,
+                const scalarField& nearestDistSqr,
+                List<pointIndexHit>&
+            ) const;
+
+            //- Find nearest intersection on line from start to end
+            virtual void findLine
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            //- Find all intersections in order from start to end
+            virtual void findLineAll
+            (
+                const pointField& start,
+                const pointField& end,
+                List<List<pointIndexHit> >&
+            ) const;
+
+            //- Find any intersection on line from start to end
+            virtual void findLineAny
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            //- From a set of points and indices get the region
+            virtual void getRegion
+            (
+                const List<pointIndexHit>&,
+                labelList& region
+            ) const;
+
+            //- From a set of points and indices get the normal
+            virtual void getNormal
+            (
+                const List<pointIndexHit>&,
+                vectorField& normal
+            ) const;
+
+            //- Determine type (inside/outside/mixed) for point. unknown if
+            //  cannot be determined (e.g. non-manifold surface)
+            virtual void getVolumeType
+            (
+                const pointField&,
+                List<volumeType>&
+            ) const;
+
+
+        // regIOobject implementation
+
+            virtual bool writeData(Ostream&) const
+            {
+                notImplemented("searchableCone::writeData(Ostream&) const");
+                return false;
+            }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/searchableCylinder.H b/src/meshTools/searchableSurface/searchableCylinder.H
index fe00bc4258..60880aa0cf 100644
--- a/src/meshTools/searchableSurface/searchableCylinder.H
+++ b/src/meshTools/searchableSurface/searchableCylinder.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -43,8 +43,6 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
-
 /*---------------------------------------------------------------------------*\
                            Class searchableCylinder Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/meshTools/searchableSurface/searchableRotatedBox.C b/src/meshTools/searchableSurface/searchableRotatedBox.C
index cbedeeeeed..3e9660ca80 100644
--- a/src/meshTools/searchableSurface/searchableRotatedBox.C
+++ b/src/meshTools/searchableSurface/searchableRotatedBox.C
@@ -55,8 +55,7 @@ Foam::searchableRotatedBox::searchableRotatedBox
             io.db(),
             io.readOpt(),
             io.writeOpt(),
-            false,      //io.registerObject(),
-            io.globalObject()
+            false      //io.registerObject(),
         ),
         treeBoundBox(point::zero, dict.lookup("span"))
     ),
diff --git a/src/meshTools/searchableSurface/searchableRotatedBox.H b/src/meshTools/searchableSurface/searchableRotatedBox.H
index 1d1b9c452c..5c5c9206ad 100644
--- a/src/meshTools/searchableSurface/searchableRotatedBox.H
+++ b/src/meshTools/searchableSurface/searchableRotatedBox.H
@@ -56,10 +56,8 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
-
 /*---------------------------------------------------------------------------*\
-                           Class searchableRotatedBox Declaration
+                    Class searchableRotatedBox Declaration
 \*---------------------------------------------------------------------------*/
 
 class searchableRotatedBox
@@ -145,6 +143,9 @@ public:
 
         // Single point queries.
 
+            //- Inherit findNearest from searchableSurface
+            using searchableSurface::findNearest;
+
             //- Calculate nearest point on surface.
             //  Returns
             //  - bool : any point found nearer than nearestDistSqr
diff --git a/src/meshTools/searchableSurface/subTriSurfaceMesh.C b/src/meshTools/searchableSurface/subTriSurfaceMesh.C
new file mode 100644
index 0000000000..dfffbd793e
--- /dev/null
+++ b/src/meshTools/searchableSurface/subTriSurfaceMesh.C
@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "subTriSurfaceMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(subTriSurfaceMesh, 0);
+addToRunTimeSelectionTable(searchableSurface, subTriSurfaceMesh, dict);
+
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::wordList Foam::subTriSurfaceMesh::patchNames(const triSurface& s)
+{
+    const geometricSurfacePatchList& patches = s.patches();
+
+    wordList names(patches.size());
+    forAll(patches, patchI)
+    {
+        names[patchI] = patches[patchI].name();
+    }
+    return names;
+}
+
+
+Foam::labelList Foam::subTriSurfaceMesh::selectedRegions
+(
+    const triSurface& s,
+    const wordReList& regionNames
+)
+{
+    const wordList names(patchNames(s));
+
+    labelList regions(names.size());
+
+    label compactI = 0;
+
+    forAll(names, regionI)
+    {
+        const word& name = names[regionI];
+
+        forAll(regionNames, i)
+        {
+            if (regionNames[i].match(name))
+            {
+                regions[compactI++] = regionI;
+            }
+        }
+    }
+    regions.setSize(compactI);
+
+    return regions;
+}
+
+
+Foam::triSurface Foam::subTriSurfaceMesh::subset
+(
+    const IOobject& io,
+    const dictionary& dict
+)
+{
+    const word subGeomName(dict.lookup("surface"));
+
+    const triSurfaceMesh& s =
+        io.db().lookupObject<triSurfaceMesh>(subGeomName);
+
+    const wordReList regionNames(dict.lookup("patches"));
+
+    labelList regionMap(selectedRegions(s, regionNames));
+
+    if (regionMap.size() == 0)
+    {
+        FatalIOErrorIn
+        (
+            "subTriSurfaceMesh::subset"
+            "(\n"
+            "    const IOobject&,\n"
+            "    const dictionary&\n"
+            ")",
+            dict
+        )   << "Found no regions in triSurface matching " << regionNames
+            << ". Valid regions are " << patchNames(s)
+            << exit(FatalIOError);
+    }
+
+    labelList reverseRegionMap(s.patches().size(), -1);
+    forAll(regionMap, i)
+    {
+        reverseRegionMap[regionMap[i]] = i;
+    }
+
+    boolList isSelected(s.size(), false);
+    forAll(s, triI)
+    {
+        if (reverseRegionMap[s.triSurface::operator[](triI).region()] != -1)
+        {
+            isSelected[triI] = true;
+        }
+    }
+
+    labelList pointMap;
+    labelList faceMap;
+    return s.subsetMesh(isSelected, pointMap, faceMap);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::subTriSurfaceMesh::subTriSurfaceMesh
+(
+    const IOobject& io,
+    const dictionary& dict
+)
+:
+    triSurfaceMesh(io, subset(io, dict))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::subTriSurfaceMesh::~subTriSurfaceMesh()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/subTriSurfaceMesh.H b/src/meshTools/searchableSurface/subTriSurfaceMesh.H
new file mode 100644
index 0000000000..9d0f508e14
--- /dev/null
+++ b/src/meshTools/searchableSurface/subTriSurfaceMesh.H
@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::subTriSurfaceMesh
+
+Description
+    On-the-fly subset of triSurfaceMesh (based on selected patches)
+
+Note
+
+    Example usage:
+    \verbatim
+    spheres.stl
+    {
+        type subTriSurfaceMesh;
+
+        // Surface to operate on. Avoid duplicate loading
+        surface spheres.stl;
+
+        // Regions to operate on (regular expressions allowed)
+        patches (solid1);
+    }
+    \endverbatim
+
+
+SourceFiles
+    subTriSurfaceMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef subTriSurfaceMesh_H
+#define subTriSurfaceMesh_H
+
+#include "triSurfaceMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class subTriSurfaceMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class subTriSurfaceMesh
+:
+    public triSurfaceMesh
+{
+    // Private member functions
+
+        //- Extract patch names of triSurface
+        static wordList patchNames(const triSurface& s);
+
+        //- Select regions by name
+        static labelList selectedRegions
+        (
+            const triSurface& s,
+            const wordReList& regionNames
+        );
+
+        //- Subset triSurface based on regions
+        static triSurface subset(const IOobject&, const dictionary&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("subTriSurfaceMesh");
+
+
+    // Constructors
+
+        //- Construct from IO and dictionary.
+        subTriSurfaceMesh
+        (
+            const IOobject& io,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~subTriSurfaceMesh();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C b/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
index 9928f39f47..0d19394afd 100644
--- a/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
+++ b/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
         {
diff --git a/tutorials/mesh/snappyHexMesh/Allrun b/tutorials/mesh/snappyHexMesh/Allrun
index d37c93a33d..aeac877153 100755
--- a/tutorials/mesh/snappyHexMesh/Allrun
+++ b/tutorials/mesh/snappyHexMesh/Allrun
@@ -12,6 +12,11 @@ cd ${0%/*} || exit 1    # Run from this directory
     ./Allrun
 )
 
+(
+    cd gap_detection || exit
+    ./Allrun
+)
+
 exit 0
 
 # These cases are links to solver test cases and are run when the Allrun
diff --git a/tutorials/mesh/snappyHexMesh/gap_detection/Allclean b/tutorials/mesh/snappyHexMesh/gap_detection/Allclean
new file mode 100755
index 0000000000..8bfe3dc5d7
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/gap_detection/Allclean
@@ -0,0 +1,7 @@
+#!/bin/sh
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+\rm -f constant/polyMesh/boundary > /dev/null 2>&1
+cleanCase
diff --git a/tutorials/mesh/snappyHexMesh/gap_detection/Allrun b/tutorials/mesh/snappyHexMesh/gap_detection/Allrun
new file mode 100755
index 0000000000..c2975fc527
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/gap_detection/Allrun
@@ -0,0 +1,15 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+
+runApplication decomposePar
+
+runParallel snappyHexMesh 8 -overwrite
+
+runApplication reconstructParMesh -constant
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/mesh/snappyHexMesh/gap_detection/README.txt b/tutorials/mesh/snappyHexMesh/gap_detection/README.txt
new file mode 100644
index 0000000000..883972db8b
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/gap_detection/README.txt
@@ -0,0 +1,7 @@
+Testcase for automatic gap refinement. The geometry is two nested boxes
+with a single small pipe between them.
+
+- 1 cell initial mesh
+- consistent normal orientation of surface so
+  specify a 'gapMode' to limit refinement only to
+  gaps on the 'outside' of the surface
diff --git a/tutorials/mesh/snappyHexMesh/gap_detection/constant/polyMesh/blockMeshDict b/tutorials/mesh/snappyHexMesh/gap_detection/constant/polyMesh/blockMeshDict
new file mode 100644
index 0000000000..94bddf09d3
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/gap_detection/constant/polyMesh/blockMeshDict
@@ -0,0 +1,102 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+vertices
+(
+    (-400 -100  -400)
+    ( 400 -100  -400)
+    ( 400  700  -400)
+    (-400  700  -400)
+    (-400 -100   400)
+    ( 400 -100   400)
+    ( 400  700   400)
+    (-400  700   400)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (1 1 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    maxY
+    {
+        type patch;
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+
+    minX
+    {
+        type patch;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+
+    maxX
+    {
+        type patch;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+
+    minY
+    {
+        type patch;
+        faces
+        (
+            (1 5 4 0)
+        );
+    }
+
+    ground
+    {
+        type patch;
+        faces
+        (
+            (0 3 2 1)
+        );
+    }
+
+    maxZ
+    {
+        type patch;
+        faces
+        (
+            (4 5 6 7)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/gap_detection/constant/triSurface/mech_test.obj b/tutorials/mesh/snappyHexMesh/gap_detection/constant/triSurface/mech_test.obj
new file mode 100644
index 0000000000..d078954dc0
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/gap_detection/constant/triSurface/mech_test.obj
@@ -0,0 +1,124 @@
+# Wavefront OBJ file
+# Regions:
+#     0    patch0
+#
+# points    : 40
+# triangles : 76
+#
+v 80 20 -20
+v -80 20 -20
+v 100 0 0
+v -100 0 0
+v 80 20 -280
+v -80 20 -280
+v 80 280 -100
+v -80 280 -100
+v 100 0 -300
+v -100 0 -300
+v 100 300 -120
+v -100 300 -120
+v 80 280 -280
+v -80 280 -280
+v 100 300 -300
+v -100 300 -300
+v 0 445 0
+v 0 445 -20
+v 3.53538 446.464 0
+v -3.53538 446.464 0
+v -3.53538 446.464 -20
+v 3.53538 446.464 -20
+v 5 450 0
+v -5 450 0
+v -5 450 -20
+v 5 450 -20
+v -3.53538 453.536 0
+v 3.53538 453.536 0
+v -3.53538 453.536 -20
+v 3.53538 453.536 -20
+v 0 455 0
+v 0 455 -20
+v 80 580 -20
+v -80 580 -20
+v 80 580 -100
+v -80 580 -100
+v 100 600 0
+v -100 600 0
+v 100 600 -120
+v -100 600 -120
+g patch0
+f 39 37 3
+f 4 3 37
+f 11 39 3
+f 15 11 3
+f 9 15 3
+f 10 9 3
+f 10 3 4
+f 40 37 39
+f 19 4 37
+f 38 37 40
+f 27 37 38
+f 23 19 37
+f 28 23 37
+f 31 28 37
+f 27 31 37
+f 40 39 11
+f 16 11 15
+f 12 11 16
+f 40 11 12
+f 16 15 9
+f 16 9 10
+f 12 10 4
+f 38 12 4
+f 24 38 4
+f 20 24 4
+f 17 20 4
+f 19 17 4
+f 12 16 10
+f 38 40 12
+f 27 38 24
+f 25 24 20
+f 29 27 24
+f 29 24 25
+f 21 20 17
+f 21 25 20
+f 18 17 19
+f 21 17 18
+f 22 19 23
+f 18 19 22
+f 26 23 28
+f 22 23 26
+f 30 28 31
+f 26 28 30
+f 32 31 27
+f 30 31 32
+f 32 27 29
+f 7 5 1
+f 2 1 5
+f 33 7 1
+f 26 33 1
+f 34 1 2
+f 21 1 34
+f 22 26 1
+f 18 22 1
+f 21 18 1
+f 7 13 5
+f 6 5 13
+f 6 2 5
+f 14 13 7
+f 6 13 14
+f 33 35 7
+f 8 7 35
+f 14 7 8
+f 34 35 33
+f 36 35 34
+f 8 35 36
+f 30 34 33
+f 26 30 33
+f 36 34 2
+f 8 36 2
+f 14 8 2
+f 6 14 2
+f 29 25 34
+f 21 34 25
+f 32 29 34
+f 30 32 34
diff --git a/tutorials/mesh/snappyHexMesh/gap_detection/system/controlDict b/tutorials/mesh/snappyHexMesh/gap_detection/system/controlDict
new file mode 100644
index 0000000000..2d847ac651
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/gap_detection/system/controlDict
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev   s                               |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     snappyHexMesh;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         2000;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   100;
+
+purgeWrite      0;
+
+writeFormat     binary;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/gap_detection/system/decomposeParDict b/tutorials/mesh/snappyHexMesh/gap_detection/system/decomposeParDict
new file mode 100644
index 0000000000..773056e68b
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/gap_detection/system/decomposeParDict
@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 8;
+
+method          hierarchical;
+
+simpleCoeffs
+{
+    n               ( 4 1 1 );
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               ( 2 2 2 );
+    delta           0.001;
+    order           xyz;
+}
+
+manualCoeffs
+{
+    dataFile        "cellDecomposition";
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/gap_detection/system/fvSchemes b/tutorials/mesh/snappyHexMesh/gap_detection/system/fvSchemes
new file mode 100644
index 0000000000..bba4ac8e0c
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/gap_detection/system/fvSchemes
@@ -0,0 +1,63 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+
+    div(phi,U)      bounded Gauss upwind;
+    div(phi,T)      bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div(phi,R)      bounded Gauss upwind;
+    div(R)          Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear limited corrected 0.333;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         limited corrected 0.333;
+}
+
+fluxRequired
+{
+    default         no;
+    p_rgh           ;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/gap_detection/system/fvSolution b/tutorials/mesh/snappyHexMesh/gap_detection/system/fvSolution
new file mode 100644
index 0000000000..627615c534
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/gap_detection/system/fvSolution
@@ -0,0 +1,69 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p_rgh
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-08;
+        relTol          0.01;
+    }
+
+    "(U|T|k|epsilon)"
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-07;
+        relTol          0.1;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 2;
+    pRefCell        0;
+    pRefValue       0;
+
+    residualControl
+    {
+        p_rgh           1e-2;
+        U               1e-4;
+        T               1e-3;
+
+        // possibly check turbulence fields
+        "(k|epsilon|omega)" 1e-3;
+    }
+}
+
+relaxationFactors
+{
+    fields
+    {
+        p_rgh           0.7;
+    }
+    equations
+    {
+        U               0.2;
+        T               0.5;
+        "(k|epsilon)"   0.7;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/gap_detection/system/meshQualityDict b/tutorials/mesh/snappyHexMesh/gap_detection/system/meshQualityDict
new file mode 100644
index 0000000000..fd0701a01c
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/gap_detection/system/meshQualityDict
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      meshQualityDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Include defaults parameters from master dictionary
+#include "$WM_PROJECT_DIR/etc/caseDicts/meshQualityDict"
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/gap_detection/system/snappyHexMeshDict b/tutorials/mesh/snappyHexMesh/gap_detection/system/snappyHexMeshDict
new file mode 100644
index 0000000000..6da20c38ca
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/gap_detection/system/snappyHexMeshDict
@@ -0,0 +1,297 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      snappyHexMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Which of the steps to run
+castellatedMesh true;
+snap            false;
+addLayers       false;
+
+
+// Geometry. Definition of all surfaces. All surfaces are of class
+// searchableSurface.
+// Surfaces are used
+// - to specify refinement for any mesh cell intersecting it
+// - to specify refinement for any mesh cell inside/outside/near
+// - to 'snap' the mesh boundary to the surface
+geometry
+{
+    mech_test.obj
+    {
+        type triSurfaceMesh;
+    }
+    all
+    {
+        type searchableBox;
+        min (-1000 -1000 -1000);
+        max (1000 1000 1000);
+    }
+
+};
+
+
+
+// Settings for the castellatedMesh generation.
+castellatedMeshControls
+{
+
+    // Refinement parameters
+    // ~~~~~~~~~~~~~~~~~~~~~
+
+    // If local number of cells is >= maxLocalCells on any processor
+    // switches from from refinement followed by balancing
+    // (current method) to (weighted) balancing before refinement.
+    maxLocalCells 100000;
+
+    // Overall cell limit (approximately). Refinement will stop immediately
+    // upon reaching this number so a refinement level might not complete.
+    // Note that this is the number of cells before removing the part which
+    // is not 'visible' from the keepPoint. The final number of cells might
+    // actually be a lot less.
+    maxGlobalCells 2000000;
+
+    // The surface refinement loop might spend lots of iterations refining just a
+    // few cells. This setting will cause refinement to stop if <= minimumRefine
+    // are selected for refinement. Note: it will at least do one iteration
+    // (unless the number of cells to refine is 0)
+    minRefinementCells 0;
+
+    // Number of buffer layers between different levels.
+    // 1 means normal 2:1 refinement restriction, larger means slower
+    // refinement.
+    nCellsBetweenLevels 1;
+
+
+
+    // Explicit feature edge refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies a level for any cell intersected by its edges.
+    // This is a featureEdgeMesh, read from constant/triSurface for now.
+    features
+    (
+    );
+
+
+
+    // Surface based refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies two levels for every surface. The first is the minimum level,
+    // every cell intersecting a surface gets refined up to the minimum level.
+    // The second level is the maximum level. Cells that 'see' multiple
+    // intersections where the intersections make an
+    // angle > resolveFeatureAngle get refined up to the maximum level.
+
+    refinementSurfaces
+    {
+        mech_test.obj
+        {
+            // Surface-wise min and max refinement level
+            level (0 0);
+        }
+    }
+
+    // Resolve sharp angles
+    resolveFeatureAngle 60;
+
+
+    // Region-wise refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies refinement level for cells in relation to a surface. One of
+    // three modes
+    // - distance. 'levels' specifies per distance to the surface the
+    //   wanted refinement level. The distances need to be specified in
+    //   descending order.
+    // - inside. 'levels' is only one entry and only the level is used. All
+    //   cells inside the surface get refined up to the level. The surface
+    //   needs to be closed for this to be possible.
+    // - outside. Same but cells outside.
+
+    refinementRegions
+    {
+        all
+        {
+            mode        inside;
+
+            // Dummy base level
+            levels      ((10000 0));
+
+            // If cells
+            // - have level 0..9
+            // - and are in a gap < 3 cell sizes across
+            // - with the gap on the inside ('inside'), outside ('outside')
+            //   or both ('mixed') of the surface
+            // refine them
+            gapLevel    (4 0 10);
+            gapMode     outside;
+        }
+    }
+
+
+    // Mesh selection
+    // ~~~~~~~~~~~~~~
+
+    // After refinement patches get added for all refinementSurfaces and
+    // all cells intersecting the surfaces get put into these patches. The
+    // section reachable from the locationInMesh is kept.
+    // NOTE: This point should never be on a face, always inside a cell, even
+    // after refinement.
+    locationInMesh (-100 -5 -300);
+
+
+    // Whether any faceZones (as specified in the refinementSurfaces)
+    // are only on the boundary of corresponding cellZones or also allow
+    // free-standing zone faces. Not used if there are no faceZones.
+    allowFreeStandingZoneFaces false;
+}
+
+
+
+// Settings for the snapping.
+snapControls
+{
+    //- Number of patch smoothing iterations before finding correspondence
+    //  to surface
+    nSmoothPatch 3;
+
+    //- Relative distance for points to be attracted by surface feature point
+    //  or edge. True distance is this factor times local
+    //  maximum edge length.
+    tolerance 2.0;
+
+    //- Number of mesh displacement relaxation iterations.
+    nSolveIter 30;
+
+    //- Maximum number of snapping relaxation iterations. Should stop
+    //  before upon reaching a correct mesh.
+    nRelaxIter 5;
+
+
+    // Feature snapping
+
+        //- Number of feature edge snapping iterations.
+        //  Leave out altogether to disable.
+        nFeatureSnapIter 10;
+
+        //- Detect (geometric) features by sampling the surface (default=false)
+        implicitFeatureSnap true;
+
+        //- Use castellatedMeshControls::features (default = true)
+        explicitFeatureSnap false;
+}
+
+
+
+// Settings for the layer addition.
+addLayersControls
+{
+    // Are the thickness parameters below relative to the undistorted
+    // size of the refined cell outside layer (true) or absolute sizes (false).
+    relativeSizes true;
+
+    // Per final patch (so not geometry!) the layer information
+    layers
+    {
+    }
+
+    // Expansion factor for layer mesh
+    expansionRatio 1.0;
+
+    // Wanted thickness of final added cell layer. If multiple layers
+    // is the thickness of the layer furthest away from the wall.
+    // Relative to undistorted size of cell outside layer.
+    // is the thickness of the layer furthest away from the wall.
+    // See relativeSizes parameter.
+    finalLayerThickness 0.5;
+
+    // Minimum thickness of cell layer. If for any reason layer
+    // cannot be above minThickness do not add layer.
+    // Relative to undistorted size of cell outside layer.
+    // See relativeSizes parameter.
+    minThickness 0.25;
+
+    // If points get not extruded do nGrow layers of connected faces that are
+    // also not grown. This helps convergence of the layer addition process
+    // close to features.
+    // Note: changed(corrected) w.r.t 17x! (didn't do anything in 17x)
+    nGrow 0;
+
+
+    // Advanced settings
+
+    // When not to extrude surface. 0 is flat surface, 90 is when two faces
+    // are perpendicular
+    featureAngle 60;
+
+    // Maximum number of snapping relaxation iterations. Should stop
+    // before upon reaching a correct mesh.
+    nRelaxIter 5;
+
+    // Number of smoothing iterations of surface normals
+    nSmoothSurfaceNormals 1;
+
+    // Number of smoothing iterations of interior mesh movement direction
+    nSmoothNormals 3;
+
+    // Smooth layer thickness over surface patches
+    nSmoothThickness 10;
+
+    // Stop layer growth on highly warped cells
+    maxFaceThicknessRatio 0.5;
+
+    // Reduce layer growth where ratio thickness to medial
+    // distance is large
+    maxThicknessToMedialRatio 0.3;
+
+    // Angle used to pick up medial axis points
+    // Note: changed(corrected) w.r.t 16x! 90 degrees corresponds to 130 in 16x.
+    minMedianAxisAngle 90;
+
+    // Create buffer region for new layer terminations
+    nBufferCellsNoExtrude 0;
+
+
+    // Overall max number of layer addition iterations. The mesher will exit
+    // if it reaches this number of iterations; possibly with an illegal
+    // mesh.
+    nLayerIter 50;
+}
+
+
+
+// Generic mesh quality settings. At any undoable phase these determine
+// where to undo.
+meshQualityControls
+{
+    #include "meshQualityDict"
+
+    // Advanced
+
+    //- Number of error distribution iterations
+    nSmoothScale 4;
+    //- amount to scale back displacement at error points
+    errorReduction 0.75;
+}
+
+
+// Advanced
+
+// Merge tolerance. Is fraction of overall bounding box of initial mesh.
+// Note: the write tolerance needs to be higher than this.
+mergeTolerance 1e-6;
+
+// ************************************************************************* //