diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index ea073cd9f0..6b40b2ab5e 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -244,8 +244,17 @@ 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>);
         //}
     }
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
index e952af9fd6..19f8a03fbd 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
@@ -109,6 +109,7 @@ class autoRefineDriver
         label bigGapOnlyRefine
         (
             const refinementParameters& refineParams,
+            const bool spreadGapSize,
             const label maxIter
         );
 
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index b2cbbd983f..2e71d97579 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -55,6 +55,7 @@ SourceFiles
 #include "wordPairHashTable.H"
 #include "surfaceZonesInfo.H"
 #include "volumeType.H"
+#include "DynamicField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -353,8 +354,23 @@ private:
                 label& nRefine
             ) const;
 
-            //- Generate single ray for cell to detect gap
-            bool generateRay
+            //- 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,
 
@@ -366,10 +382,13 @@ private:
                 const point& cc,
                 const label cLevel,
 
-                point& start,
-                point& end,
-                point& start2,
-                point& end2
+                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
@@ -400,9 +419,12 @@ private:
             label markInternalGapRefinement
             (
                 const scalar planarCos,
+                const bool spreadGapSize,
                 const label nAllowRefine,
                 labelList& refineCell,
-                label& nRefine
+                label& nRefine,
+                labelList& numGapCells,
+                scalarField& gapSize
             ) const;
 
             //- Refine cells containing small gaps
@@ -958,6 +980,7 @@ public:
                 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/meshRefinementGapRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementGapRefine.C
index 62a285928c..cd9cd153db 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementGapRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementGapRefine.C
@@ -28,11 +28,14 @@ License
 #include "refinementSurfaces.H"
 #include "refinementFeatures.H"
 #include "shellSurfaces.H"
-#include "OBJstream.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  * * * * * * * * * * * //
 
@@ -71,17 +74,18 @@ Foam::label Foam::meshRefinement::markSurfaceGapRefinement
 
         pointField start(testFaces.size());
         pointField end(testFaces.size());
-        labelList minLevel(testFaces.size());
-
-        calcCellCellRays
-        (
-            neiCc,
-            neiLevel,
-            testFaces,
-            start,
-            end,
-            minLevel
-        );
+        {
+            labelList minLevel(testFaces.size());
+            calcCellCellRays
+            (
+                neiCc,
+                neiLevel,
+                testFaces,
+                start,
+                end,
+                minLevel
+            );
+        }
 
 
         // Collect cells to test for inside/outside in shell
@@ -134,24 +138,21 @@ Foam::label Foam::meshRefinement::markSurfaceGapRefinement
         }
 
 
-        //OBJstream str(mesh_.time().timePath()/"markSurfaceRefinement.obj");
-        //Info<< "Dumping rays to " << str.name( ) << endl;
-
         const List<FixedList<label, 3> >& extendedGapLevel =
             surfaces_.extendedGapLevel();
         const List<volumeType>& extendedGapMode =
             surfaces_.extendedGapMode();
 
-        labelList surface1;
-        List<pointIndexHit> hit1;
-        labelList region1;
-        vectorField normal1;
+        labelList ccSurface1;
+        List<pointIndexHit> ccHit1;
+        labelList ccRegion1;
+        vectorField ccNormal1;
 
         {
-            labelList surface2;
-            List<pointIndexHit> hit2;
-            labelList region2;
-            vectorField normal2;
+            labelList ccSurface2;
+            List<pointIndexHit> ccHit2;
+            labelList ccRegion2;
+            vectorField ccNormal2;
 
             surfaces_.findNearestIntersection
             (
@@ -159,38 +160,43 @@ Foam::label Foam::meshRefinement::markSurfaceGapRefinement
                 start,
                 end,
 
-                surface1,
-                hit1,
-                region1,
-                normal1,
+                ccSurface1,
+                ccHit1,
+                ccRegion1,
+                ccNormal1,
 
-                surface2,
-                hit2,
-                region2,
-                normal2
+                ccSurface2,
+                ccHit2,
+                ccRegion2,
+                ccNormal2
             );
         }
 
+        start.clear();
+        end.clear();
 
-        start.setSize(2*surface1.size());
-        end.setSize(2*surface1.size());
-        labelList map(2*surface1.size());
-        pointField start2(2*surface1.size());
-        pointField end2(2*surface1.size());
-        labelList cellMap(2*surface1.size());
+        DynamicField<point> rayStart(2*ccSurface1.size());
+        DynamicField<point> rayEnd(2*ccSurface1.size());
+        DynamicField<scalar> gapSize(2*ccSurface1.size());
 
-        label compactI = 0;
+        DynamicField<point> rayStart2(2*ccSurface1.size());
+        DynamicField<point> rayEnd2(2*ccSurface1.size());
+        DynamicField<scalar> gapSize2(2*ccSurface1.size());
 
-        forAll(surface1, i)
+        DynamicList<label> cellMap(2*ccSurface1.size());
+        DynamicList<label> compactMap(2*ccSurface1.size());
+
+        forAll(ccSurface1, i)
         {
-            label surfI = surface1[i];
+            label surfI = ccSurface1[i];
 
             if (surfI != -1)
             {
-                label globalRegionI = surfaces_.globalRegion(surfI, region1[i]);
+                label globalRegionI =
+                    surfaces_.globalRegion(surfI, ccRegion1[i]);
 
                 label faceI = testFaces[i];
-                const point& surfPt = hit1[i].hitPoint();
+                const point& surfPt = ccHit1[i].hitPoint();
 
                 label own = mesh_.faceOwner()[faceI];
                 if (cellToCompact[own] != -1)
@@ -210,26 +216,28 @@ Foam::label Foam::meshRefinement::markSurfaceGapRefinement
                     );
 
                     const point& cc = cellCentres[own];
-                    bool okRay = generateRay
+                    label nRays = generateRays
                     (
                         false,
                         surfPt,
-                        normal1[i],
+                        ccNormal1[i],
                         gapInfo,
                         gapMode,
-                        surfPt+((cc-surfPt)&normal1[i])*normal1[i],
+                        surfPt+((cc-surfPt)&ccNormal1[i])*ccNormal1[i],
                         cellLevel[own],
 
-                        start[compactI],
-                        end[compactI],
-                        start2[compactI],
-                        end2[compactI]
-                    );
-                    cellMap[compactI] = own;
+                        rayStart,
+                        rayEnd,
+                        gapSize,
 
-                    if (okRay)
+                        rayStart2,
+                        rayEnd2,
+                        gapSize2
+                    );
+                    for (label j = 0; j < nRays; j++)
                     {
-                        map[compactI++] = i;
+                        cellMap.append(own);
+                        compactMap.append(i);
                     }
                 }
                 if (mesh_.isInternalFace(faceI))
@@ -252,31 +260,37 @@ Foam::label Foam::meshRefinement::markSurfaceGapRefinement
                         );
 
                         const point& cc = cellCentres[nei];
-                        bool okRay = generateRay
+                        label nRays = generateRays
                         (
                             false,
                             surfPt,
-                            normal1[i],
+                            ccNormal1[i],
                             gapInfo,
                             gapMode,
-                            surfPt+((cc-surfPt)&normal1[i])*normal1[i],
+                            surfPt+((cc-surfPt)&ccNormal1[i])*ccNormal1[i],
                             cellLevel[nei],
 
-                            start[compactI],
-                            end[compactI],
-                            start2[compactI],
-                            end2[compactI]
-                        );
-                        cellMap[compactI] = nei;
+                            rayStart,
+                            rayEnd,
+                            gapSize,
 
-                        if (okRay)
+                            rayStart2,
+                            rayEnd2,
+                            gapSize2
+                        );
+                        for (label j = 0; j < nRays; j++)
                         {
-                            map[compactI++] = i;
+                            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)
@@ -296,99 +310,103 @@ Foam::label Foam::meshRefinement::markSurfaceGapRefinement
                         );
 
                         const point& cc = neiCc[bFaceI];
-                        bool okRay = generateRay
+                        label nRays = generateRays
                         (
                             false,
                             surfPt,
-                            normal1[i],
+                            ccNormal1[i],
                             gapInfo,
                             gapMode,
-                            surfPt+((cc-surfPt)&normal1[i])*normal1[i],
+                            surfPt+((cc-surfPt)&ccNormal1[i])*ccNormal1[i],
                             neiLevel[bFaceI],
 
-                            start[compactI],
-                            end[compactI],
-                            start2[compactI],
-                            end2[compactI]
-                        );
-                        cellMap[compactI] = -1;
+                            rayStart,
+                            rayEnd,
+                            gapSize,
 
-                        if (okRay)
+                            rayStart2,
+                            rayEnd2,
+                            gapSize2
+                        );
+                        for (label j = 0; j < nRays; j++)
                         {
-                            map[compactI++] = i;
+                            cellMap.append(-1); // See above.
+                            compactMap.append(i);
                         }
                     }
                 }
             }
         }
 
-        //Info<< "Retesting " << returnReduce(compactI, sumOp<label>())
-        //    << " out of " << returnReduce(start.size(), sumOp<label>())
-        //    << endl;
+        Info<< "Shooting " << returnReduce(rayStart.size(), sumOp<label>())
+            << " rays from " << returnReduce(testFaces.size(), sumOp<label>())
+            << " intersected faces" << endl;
 
-        start.setSize(compactI);
-        end.setSize(compactI);
-        start2.setSize(compactI);
-        end2.setSize(compactI);
-        map.setSize(compactI);
-        cellMap.setSize(compactI);
+        rayStart.shrink();
+        rayEnd.shrink();
+        gapSize.shrink();
 
-        testFaces = UIndirectList<label>(testFaces, map)();
-        minLevel = UIndirectList<label>(minLevel, map)();
-        surface1.clear();
-        hit1.clear();
-        region1.clear();
-        normal1 = UIndirectList<vector>(normal1, map)();
+        rayStart2.shrink();
+        rayEnd2.shrink();
+        gapSize2.shrink();
 
-        // Do intersections in first direction
-        labelList surfaceHit;
-        vectorField surfaceNormal;
+        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
         (
-            start,
-            end,
-            surfaceHit,
-            surfaceNormal
+            rayStart,
+            rayEnd,
+            surf1,
+            hit1,
+            normal1
         );
+
+        labelList surf2;
+        List<pointIndexHit> hit2;
+        vectorField normal2;
+        surfaces_.findNearestIntersection
+        (
+            rayStart2,
+            rayEnd2,
+            surf2,
+            hit2,
+            normal2
+        );
+
+        forAll(surf1, i)
         {
-            // Do intersections in second direction and merge
-            labelList surfaceHit2;
-            vectorField surfaceNormal2;
-            surfaces_.findNearestIntersection
-            (
-                start2,
-                end2,
-                surfaceHit2,
-                surfaceNormal2
-            );
-            forAll(surfaceHit, i)
-            {
-                if (surfaceHit[i] == -1 && surfaceHit2[i] != -1)
-                {
-                    surfaceHit[i] = surfaceHit2[i];
-                    surfaceNormal[i] = surfaceNormal2[i];
-                }
-            }
-        }
-
-
-        forAll(surfaceHit, i)
-        {
-            label surfI = surfaceHit[i];
-
-            if (surfI != -1)
+            if (surf1[i] != -1 && surf2[i] != -1)
             {
                 // Found intersection with surface. Check opposite normal.
-
                 label cellI = cellMap[i];
 
-                if (cellI != -1 && mag(normal1[i]&surfaceNormal[i]) > planarCos)
+                if
+                (
+                    cellI != -1
+                 && (mag(normal1[i]&normal2[i]) > planarCos)
+                 && (
+                        magSqr(hit1[i].hitPoint()-hit2[i].hitPoint())
+                      < Foam::sqr(gapSize[i])
+                    )
+                )
                 {
                     if
                     (
                        !markForRefine
                         (
-                            surfI,
+                            surf1[i],
                             nAllowRefine,
                             refineCell[cellI],
                             nRefine
@@ -415,84 +433,6 @@ Foam::label Foam::meshRefinement::markSurfaceGapRefinement
 }
 
 
-//Foam::labelList Foam::meshRefinement::extractHits
-//(
-//    const List<pointIndexHit>& info
-//)
-//{
-//    labelList compactMap(info.size());
-//    label compactI = 0;
-//    forAll(info, i)
-//    {
-//        if (info[i].hit())
-//        {
-//            compactMap[compactI++] = i;
-//        }
-//    }
-//    compactMap.setSize(compactI);
-//    return compactMap;
-//}
-//
-//
-//void Foam::meshRefinement::mergeHits
-//(
-//    const pointField& samples,
-//    const List<pointIndexHit>& extraInfo,
-//    List<pointIndexHit>& info
-//)
-//{
-//    if (extraInfo.size() != info.size())
-//    {
-//        FatalErrorIn
-//        (
-//            "meshRefinement::mergeHits(const List<pointIndexHit>&"
-//            ", List<pointIndexHit>& info)"
-//        )   << "Sizes differ " << extraInfo.size() << ' ' << info.size()
-//            << exit(FatalError);
-//    }
-//
-//    // Merge oppositeInfo2 into oppositeInfo
-//    forAll(extraInfo, i)
-//    {
-//        if (!info[i].hit())
-//        {
-//            if (extraInfo[i].hit())
-//            {
-//                info[i] = extraInfo[i];
-//            }
-//        }
-//        else if (extraInfo[i].hit())
-//        {
-//            const point& nearPt = samples[i];   //nearInfo[i].hitPoint();
-//
-//            scalar d = magSqr(info[i].hitPoint()-nearPt);
-//            scalar d2 = magSqr(extraInfo[i].hitPoint()-nearPt);
-//
-//            if (d2 < d)
-//            {
-//                info[i] = extraInfo[i];
-//            }
-//        }
-//    }
-//}
-//
-//
-//Foam::tmp<Foam::pointField> Foam::meshRefinement::extractPoints
-//(
-//    const List<pointIndexHit>& info
-//)
-//{
-//    tmp<pointField> tfld(new pointField(info.size()));
-//    pointField& fld = tfld();
-//
-//    forAll(info, i)
-//    {
-//        fld[i] = info[i].rawPoint();
-//    }
-//    return tfld;
-//}
-//
-//
 //Foam::meshRefinement::findNearestOppositeOp::findNearestOppositeOp
 //(
 //    const indexedOctree<treeDataTriSurface>& tree,
@@ -730,7 +670,58 @@ Foam::label Foam::meshRefinement::markSurfaceGapRefinement
 //}
 
 
-bool Foam::meshRefinement::generateRay
+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,
 
@@ -742,10 +733,13 @@ bool Foam::meshRefinement::generateRay
     const point& cc,
     const label cLevel,
 
-    point& start,
-    point& end,
-    point& start2,
-    point& end2
+    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:
@@ -805,7 +799,7 @@ bool Foam::meshRefinement::generateRay
     //   magV < 0.5cellSize)
     // - or otherwise if on the correct side
 
-    bool okRay = false;
+    label nOldRays = start.size();
 
     if (cLevel >= gapInfo[1] && cLevel < gapInfo[2])
     {
@@ -823,36 +817,52 @@ bool Foam::meshRefinement::generateRay
             const vector& n = nearNormal;
 
             // Situation 'C' above: cell too close. Use surface
-            // normal to shoot rays
+            // -normal and -point to shoot rays
 
             if (mode == volumeType::OUTSIDE)
             {
-                start = nearPoint+1e-6*n;
-                end = nearPoint+nearGap*n;
-                // Small second vector just to make sure to
-                // refine cell
-                start2 = nearPoint-1e-6*n;
-                end2 = nearPoint-0.5*cellSize*n;
+                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 = nearPoint-1e-6*n;
-                end = nearPoint-nearGap*n;
-                // Small second vector just to make sure to
-                // refine cell
-                start2 = nearPoint+1e-6*n;
-                end2 = nearPoint+0.5*cellSize*n;
+                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)
             {
-                start = nearPoint+1e-6*n;
-                end = nearPoint+nearGap*n;
-
-                start2 = nearPoint-1e-6*n;
-                end2 = nearPoint-nearGap*n;
+                // 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());
+                }
             }
-
-            okRay = true;
         }
         else
         {
@@ -871,117 +881,53 @@ bool Foam::meshRefinement::generateRay
              || (mode == volumeType::INSIDE && s < -SMALL)
             )
             {
-                // Use vector through cell centre
-                vector n(v/(magV+ROOTVSMALL));
+                //// Use 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);
 
-                start = nearPoint+1e-6*n;
-                end = nearPoint+nearGap*n;
-                // Dummy second vector
-                start2 = start;
-                end2 = start2;
 
-                okRay = true;
+                // 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);
             }
         }
     }
 
-    return okRay;
+    return start.size()-nOldRays;
 }
 
 
-//void Foam::meshRefinement::shootRays
-//(
-//    const label surfI,
-//    labelList& nearMap,
-//    scalarField& nearGap,
-//    List<pointIndexHit>& nearInfo,
-//    List<pointIndexHit>& oppositeInfo
-//) const
-//{
-//    const labelList& cellLevel = meshCutter_.cellLevel();
-//    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];
-//
-//
-//    // Normals for ray shooting and inside/outside detection
-//    vectorField nearNormal;
-//    geom.getNormal(nearInfo, nearNormal);
-//    // Regions
-//    labelList nearRegion;
-//    geom.getRegion(nearInfo, nearRegion);
-//
-//    labelList map(nearInfo.size());
-//
-//    pointField start(nearInfo.size());
-//    pointField end(nearInfo.size());
-//    pointField start2(nearInfo.size());
-//    pointField end2(nearInfo.size());
-//
-//    label compactI = 0;
-//    forAll(nearInfo, i)
-//    {
-//        label globalRegionI =
-//            surfaces_.globalRegion(surfI, nearRegion[i]);
-//
-//        label cellI = nearMap[i];
-//        label cLevel = cellLevel[cellI];
-//        scalar cellSize = edge0Len/pow(2, cLevel);
-//
-//        // Update gap size
-//        nearGap[i] = extendedGapLevel[globalRegionI][0]*cellSize;
-//
-//        // Construct one or two rays to test for oppositeness
-//        bool okRay = generateRay
-//        (
-//            false,
-//            nearInfo[i].hitPoint(),
-//            nearNormal[i],
-//            extendedGapLevel[globalRegionI],
-//            extendedGapMode[globalRegionI],
-//
-//            cellCentres[cellI],
-//            cellLevel[cellI],
-//
-//            start[compactI],
-//            end[compactI],
-//            start2[compactI],
-//            end2[compactI]
-//        );
-//        if (okRay)
-//        {
-//            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);
-//    start.setSize(compactI);
-//    end.setSize(compactI);
-//    start2.setSize(compactI);
-//    end2.setSize(compactI);
-//
-//    nearMap = UIndirectList<label>(nearMap, map)();
-//    nearGap = UIndirectList<scalar>(nearGap, map)();
-//    nearInfo = UIndirectList<pointIndexHit>(nearInfo, map)();
-//
-//    geom.findLineAny(start, end, oppositeInfo);
-//
-//    List<pointIndexHit> oppositeInfo2;
-//    geom.findLineAny(start2, end2, oppositeInfo2);
-//    mergeHits(extractPoints(nearInfo), oppositeInfo2, oppositeInfo);
-//}
-
-
 void Foam::meshRefinement::selectGapCandidates
 (
     const labelList& refineCell,
@@ -1081,12 +1027,20 @@ void Foam::meshRefinement::mergeGapInfo
 Foam::label Foam::meshRefinement::markInternalGapRefinement
 (
     const scalar planarCos,
+    const bool spreadGapSize,
     const label nAllowRefine,
 
     labelList& refineCell,
-    label& nRefine
+    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();
@@ -1116,7 +1070,6 @@ Foam::label Foam::meshRefinement::markInternalGapRefinement
             shellGapMode
         );
 
-
         // Find nearest point and normal on the surfaces
         List<pointIndexHit> nearInfo;
         vectorField nearNormal;
@@ -1147,13 +1100,16 @@ Foam::label Foam::meshRefinement::markInternalGapRefinement
 
 
 
-        labelList map(nearInfo.size());
-        pointField start(nearInfo.size());
-        pointField end(nearInfo.size());
-        pointField start2(nearInfo.size());
-        pointField end2(nearInfo.size());
+        DynamicList<label> map(nearInfo.size());
+        DynamicField<point> rayStart(nearInfo.size());
+        DynamicField<point> rayEnd(nearInfo.size());
+        DynamicField<scalar> gapSize(nearInfo.size());
 
-        label compactI = 0;
+        DynamicField<point> rayStart2(nearInfo.size());
+        DynamicField<point> rayEnd2(nearInfo.size());
+        DynamicField<scalar> gapSize2(nearInfo.size());
+
+        label nTestCells = 0;
 
         forAll(nearInfo, i)
         {
@@ -1178,9 +1134,16 @@ Foam::label Foam::meshRefinement::markInternalGapRefinement
                     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 two rays to test for oppositeness
-                bool okRay = generateRay
+                // Construct one or more rays to test for oppositeness
+                label nRays = generateRays
                 (
                     false,
                     nearInfo[i].hitPoint(),
@@ -1188,29 +1151,39 @@ Foam::label Foam::meshRefinement::markInternalGapRefinement
                     gapInfo,
                     gapMode,
 
-                    cellCentres[cellMap[i]],
-                    cellLevel[cellMap[i]],
+                    cellCentres[cellI],
+                    cLevel,
 
-                    start[compactI],
-                    end[compactI],
-                    start2[compactI],
-                    end2[compactI]
+                    rayStart,
+                    rayEnd,
+                    gapSize,
+
+                    rayStart2,
+                    rayEnd2,
+                    gapSize2
                 );
-                if (okRay)
+                if (nRays > 0)
                 {
-                    map[compactI++] = i;
+                    nTestCells++;
+                    for (label j = 0; j < nRays; j++)
+                    {
+                        map.append(i);
+                    }
                 }
             }
         }
 
-        Info<< "Selected " << returnReduce(compactI, sumOp<label>())
+        Info<< "Selected " << returnReduce(nTestCells, sumOp<label>())
             << " cells for testing out of "
             << mesh_.globalData().nTotalCells() << endl;
-        map.setSize(compactI);
-        start.setSize(compactI);
-        end.setSize(compactI);
-        start2.setSize(compactI);
-        end2.setSize(compactI);
+        map.shrink();
+        rayStart.shrink();
+        rayEnd.shrink();
+        gapSize.shrink();
+
+        rayStart2.shrink();
+        rayEnd2.shrink();
+        gapSize2.shrink();
 
         cellMap = UIndirectList<label>(cellMap, map)();
         nearNormal = UIndirectList<vector>(nearNormal, map)();
@@ -1221,64 +1194,191 @@ Foam::label Foam::meshRefinement::markInternalGapRefinement
         nearRegion.clear();
 
 
-        // Do intersections in first direction
-        labelList surfaceHit;
-        vectorField surfaceNormal;
+        // Do intersections in pairs
+        labelList surf1;
+        List<pointIndexHit> hit1;
+        vectorField normal1;
         surfaces_.findNearestIntersection
         (
-            start,
-            end,
-            surfaceHit,
-            surfaceNormal
+            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)
         {
-            // Do intersections in second direction and merge
-            labelList surfaceHit2;
-            vectorField surfaceNormal2;
-            surfaces_.findNearestIntersection
-            (
-                start2,
-                end2,
-                surfaceHit2,
-                surfaceNormal2
-            );
-            forAll(surfaceHit, i)
+            if (surf1[i] != -1 && surf2[i] != -1)
             {
-                if (surfaceHit[i] == -1 && surfaceHit2[i] != -1)
-                {
-                    surfaceHit[i] = surfaceHit2[i];
-                    surfaceNormal[i] = surfaceNormal2[i];
-                }
-            }
-        }
-
-
-        forAll(surfaceHit, i)
-        {
-            label surfI = surfaceHit[i];
-
-            if (surfI != -1 && mag(nearNormal[i]&surfaceNormal[i]) > planarCos)
-            {
-                // Found intersection with surface. Check opposite normal.
+                // Found intersections with surface. Check for
+                // - small gap
+                // - coplanar normals
 
                 label cellI = cellMap[i];
 
+                scalar d2 = magSqr(hit1[i].hitPoint()-hit2[i].hitPoint());
+
                 if
                 (
-                   !markForRefine
-                    (
-                        surfI,
-                        nAllowRefine,
-                        refineCell[cellI],
-                        nRefine
-                    )
+                    cellI != -1
+                 && (mag(normal1[i]&normal2[i]) > planarCos)
+                 && (d2 < Foam::sqr(gapSize[i]))
                 )
                 {
-                    break;
+                    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;
+
+
+            FaceCellWave<transportData> deltaCalc
+            (
+                mesh_,
+                frontFaces,
+                frontData,
+                faceData,
+                cellData,
+                mesh_.globalData().nTotalCells()+1
+            );
+
+
+            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>())
@@ -1305,8 +1405,6 @@ Foam::label Foam::meshRefinement::markSmallFeatureRefinement
 ) const
 {
     const labelList& cellLevel = meshCutter_.cellLevel();
-    const pointField& cellCentres = mesh_.cellCentres();
-
     const labelList& surfaceIndices = surfaces_.surfaces();
     const List<FixedList<label, 3> >& extendedGapLevel =
         surfaces_.extendedGapLevel();
@@ -1376,14 +1474,14 @@ Foam::label Foam::meshRefinement::markSmallFeatureRefinement
         );
 
 
-        labelList map(ctrs.size());
-        labelList cellMap(ctrs.size());
-        label compactI = 0;
+        DynamicList<label> map(ctrs.size());
+        DynamicList<label> cellMap(ctrs.size());
 
-        pointField start(ctrs.size());
-        pointField end(ctrs.size());
-        pointField start2(ctrs.size());
-        pointField end2(ctrs.size());
+        DynamicField<point> rayStart(ctrs.size());
+        DynamicField<point> rayEnd(ctrs.size());
+        DynamicField<scalar> gapSize(ctrs.size());
+
+        label nTestCells = 0;
 
         forAll(ctrs, i)
         {
@@ -1422,40 +1520,39 @@ Foam::label Foam::meshRefinement::markSmallFeatureRefinement
                     // Note that we always want to use the surface normal
                     // and not the vector from cell centre to surface point
 
-                    bool okRay = generateRay
+                    label nRays = generateRays
                     (
-                        true,               // always use surface normal
                         ctrs[i],
                         normal[i],
                         gapInfo,
                         gapMode,
 
-                        cellCentres[cellI],
                         cellLevel[cellI],
 
-                        start[compactI],
-                        end[compactI],
-                        start2[compactI],
-                        end2[compactI]
+                        rayStart,
+                        rayEnd
                     );
-                    if (okRay)
+
+                    if (nRays > 0)
                     {
-                        cellMap[compactI] = cellI;
-                        map[compactI++] = i;
+                        nTestCells++;
+                        for (label j = 0; j < nRays; j++)
+                        {
+                            cellMap.append(cellI);
+                            map.append(i);
+                        }
                     }
                 }
             }
         }
 
-        Info<< "Selected " << returnReduce(compactI, sumOp<label>())
+        Info<< "Selected " << returnReduce(nTestCells, sumOp<label>())
             << " cells containing triangle centres out of "
             << mesh_.globalData().nTotalCells() << endl;
-        map.setSize(compactI);
-        cellMap.setSize(compactI);
-        start.setSize(compactI);
-        end.setSize(compactI);
-        start2.setSize(compactI);
-        end2.setSize(compactI);
+        map.shrink();
+        cellMap.shrink();
+        rayStart.shrink();
+        rayEnd.shrink();
 
         ctrs.clear();
         region.clear();
@@ -1463,70 +1560,49 @@ Foam::label Foam::meshRefinement::markSmallFeatureRefinement
         shellGapMode.clear();
         normal = UIndirectList<vector>(normal, map)();
 
-
-        // Do intersections in first direction.
-        // Note passing in of -1 for cell level such that it does not
-        // discard surfaces with level 0. (since can be overridden with
-        // gap level)
+        // Do intersections.
         labelList surfaceHit;
         vectorField surfaceNormal;
         surfaces_.findNearestIntersection
         (
-            start,
-            end,
+            rayStart,
+            rayEnd,
             surfaceHit,
             surfaceNormal
         );
-        {
-            // Do intersections in second direction and merge
-            labelList surfaceHit2;
-            vectorField surfaceNormal2;
-            surfaces_.findNearestIntersection
-            (
-                start2,
-                end2,
-                surfaceHit2,
-                surfaceNormal2
-            );
-            forAll(surfaceHit, i)
-            {
-                if (surfaceHit[i] == -1 && surfaceHit2[i] != -1)
-                {
-                    surfaceHit[i] = surfaceHit2[i];
-                    surfaceNormal[i] = surfaceNormal2[i];
-                }
-            }
-        }
 
 
-        label nGaps = 0;
+        label nOldRefine = 0;
+
 
         forAll(surfaceHit, i)
         {
-            label surfI = surfaceHit[i];
-
-            if (surfI != -1 && mag(normal[i]&surfaceNormal[i]) > planarCos)
+            if (surfaceHit[i] != -1) // && surf2[i] != -1)
             {
-                nGaps++;
+                // Found intersection with surface. Check coplanar normals.
+                label cellI = cellMap[i];
 
-                if
-                (
-                   !markForRefine
-                    (
-                        surfI,
-                        nAllowRefine,
-                        refineCell[cellMap[i]],
-                        nRefine
-                    )
-                )
+                if (mag(normal[i]&surfaceNormal[i]) > planarCos)
                 {
-                    break;
+                    if
+                    (
+                       !markForRefine
+                        (
+                            surfaceHit[i],
+                            nAllowRefine,
+                            refineCell[cellI],
+                            nRefine
+                        )
+                    )
+                    {
+                        break;
+                    }
                 }
             }
         }
 
         Info<< "For surface " << geom.name() << " found "
-            << returnReduce(nGaps, sumOp<label>())
+            << returnReduce(nRefine-nOldRefine, sumOp<label>())
             << " cells in small gaps" << endl;
 
         if
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index 8429a3c117..4ffa1248f8 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -2055,6 +2055,7 @@ Foam::labelList Foam::meshRefinement::refineCandidates
     const bool smallFeatureRefinement,
     const bool gapRefinement,
     const bool bigGapRefinement,
+    const bool spreadGapSize,
     const label maxGlobalCells,
     const label maxLocalCells
 ) const
@@ -2290,12 +2291,18 @@ Foam::labelList Foam::meshRefinement::refineCandidates
         {
             // Refine based on gap information provided by shell and nearest
             // surface
+            labelList numGapCells;
+            scalarField gapSize;
             label nGap = markInternalGapRefinement
             (
                 planarCos,
+                spreadGapSize,
                 nAllowRefine,
+
                 refineCell,
-                nRefine
+                nRefine,
+                numGapCells,
+                gapSize
             );
             Info<< "Marked for refinement due to opposite surfaces"
                 << "             "
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/transportData.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/transportData.H
new file mode 100644
index 0000000000..f601d7e315
--- /dev/null
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/transportData.H
@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 hexRef8;
+
+/*---------------------------------------------------------------------------*\
+                           Class transportData Declaration
+\*---------------------------------------------------------------------------*/
+
+class transportData
+:
+    public wallPointData<scalar>
+{
+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..16b497f4c0
--- /dev/null
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/transportDataI.H
@@ -0,0 +1,165 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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,
+    const transportData& neighbourWallInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    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
+)
+{
+    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
+)
+{
+    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 6f36ce5904..b6fb43d98a 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -1373,6 +1373,56 @@ void Foam::refinementSurfaces::findNearestIntersection
 }
 
 
+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 72f20c1777..7e2bcdf946 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -350,6 +350,16 @@ public:
                 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
             (