diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
index 40607fbe18..7e3b43f9cb 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
@@ -555,8 +555,21 @@ private:
                 label& nRefine
             ) const;
 
-            //- Mark cells for surface proximity based refinement.
-            label markProximityRefinementWave
+            //- Mark faces for additional proximity based testing. Extends
+            //  testFaces
+            void markProximityCandidateFaces
+            (
+                const labelList& blockedSurfaces,
+                const List<scalarList>& regionToBlockSize,
+                const labelList& neiLevel,
+                const pointField& neiCc,
+
+                labelList& testFaces
+            ) const;
+
+            //- Mark cells for surface proximity based refinement. Uses
+            //  ray-shooting from markInternalGapRefinement.
+            label markFakeGapRefinement
             (
                 const scalar planarCos,
                 const labelList& blockedSurfaces,
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBlock.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBlock.C
index fc27ab08af..435cadd56d 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBlock.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBlock.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -455,16 +455,14 @@ void Foam::meshRefinement::growSet
 //}
 
 
-Foam::label Foam::meshRefinement::markProximityRefinementWave
+void Foam::meshRefinement::markProximityCandidateFaces
 (
-    const scalar planarCos,
     const labelList& blockedSurfaces,
-    const label nAllowRefine,
+    const List<scalarList>& regionToBlockSize,
     const labelList& neiLevel,
     const pointField& neiCc,
 
-    labelList& refineCell,
-    label& nRefine
+    labelList& testFaces
 ) const
 {
     labelListList faceZones(surfaces_.surfaces().size());
@@ -498,35 +496,6 @@ Foam::label Foam::meshRefinement::markProximityRefinementWave
     }
 
 
-    // Re-work the blockLevel of the blockedSurfaces into a length scale
-    // and a number of cells to cross
-    List<scalarList> regionToBlockSize(surfaces_.surfaces().size());
-    for (const label surfi : blockedSurfaces)
-    {
-        const label geomi = surfaces_.surfaces()[surfi];
-        const searchableSurface& s = surfaces_.geometry()[geomi];
-        const label nRegions = s.regions().size();
-        regionToBlockSize[surfi].setSize(nRegions);
-        for (label regioni = 0; regioni < nRegions; regioni++)
-        {
-            const label globalRegioni = surfaces_.globalRegion(surfi, regioni);
-            const label bLevel = surfaces_.blockLevel()[globalRegioni];
-            regionToBlockSize[surfi][regioni] =
-                meshCutter_.level0EdgeLength()/pow(2.0, bLevel);
-
-            //const label mLevel = surfaces_.maxLevel()[globalRegioni];
-            //// TBD: check for higher cached level of surface due to vol
-            ////      refinement. Problem: might still miss refinement bubble
-            ////      fully inside thin channel
-            //if (isA<triSurfaceMesh>(s) && !isA<distributedTriSurfaceMesh>(s))
-            //{
-            //    const triSurfaceMesh& surf = refCast<const triSurfaceMesh>(s);
-            //}
-
-            //nIters = max(nIters, (2<<(mLevel-bLevel)));
-        }
-    }
-
     // Clever limiting of the number of iterations (= max cells in the channel)
     // since it has too many problematic issues, e.g. with volume refinement
     // and the real check uses the proper distance anyway just disable.
@@ -535,7 +504,7 @@ Foam::label Foam::meshRefinement::markProximityRefinementWave
 
     // Collect candidate faces (i.e. intersecting any surface and
     // owner/neighbour not yet refined)
-    const labelList testFaces(getRefineCandidateFaces(refineCell));
+    //const labelList testFaces(getRefineCandidateFaces(refineCell));
 
     // Collect segments
     pointField start(testFaces.size());
@@ -729,208 +698,309 @@ Foam::label Foam::meshRefinement::markProximityRefinementWave
     wallDistCalc.iterate(nIters);
 
 
-    if (debug&meshRefinement::MESH)
-    {
-        // Dump current nearest opposite surfaces
-        volScalarField distance
-        (
-            IOobject
-            (
-                "gapSize",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                IOobject::NO_REGISTER
-            ),
-            mesh_,
-            dimensionedScalar
-            (
-                "zero",
-                dimLength,  //dimensionSet(0, 1, 0, 0, 0),
-                -1
-            )
-        );
+    // Extract faces of visited cells
 
-        forAll(allCellInfo, celli)
-        {
-            if (allCellInfo[celli].valid(wallDistCalc.data()))
-            {
-                const point& cc = mesh_.cellCentres()[celli];
-                // Nearest surface points
-                const List<point>& origin = allCellInfo[celli].origin();
-
-                // Find 'opposite' pair with minimum distance
-                for (label i = 0; i < origin.size(); i++)
-                {
-                    for (label j = i + 1; j < origin.size(); j++)
-                    {
-                        if (((cc-origin[i]) & (cc-origin[j])) < 0)
-                        {
-                            const scalar d(mag(origin[i]-origin[j]));
-                            if (distance[celli] < 0)
-                            {
-                                distance[celli] = d;
-                            }
-                            else
-                            {
-                                distance[celli] = min(distance[celli], d);
-                            }
-                        }
-                    }
-                }
-            }
-        }
-        distance.correctBoundaryConditions();
-
-        Info<< "Writing measured gap distance to "
-            << distance.name() << endl;
-        distance.write();
-    }
-
-
-
-    // Detect tight gaps:
-    // - cell is inbetween the two surfaces
-    // - two surfaces are planarish
-    // - two surfaces are not too far apart
-    //   (number of walking iterations is a too-coarse measure)
-
-    scalarField smallGapDistance(mesh_.nCells(), 0.0);
-    label nMulti = 0;
-    label nSmallGap = 0;
-
-    //OBJstream str(mesh_.time().timePath()/"multiRegion.obj");
+    bitSet isProximityFace(mesh_.nFaces(), false);
 
+    // Make sure to include initial intersected ones
+    isProximityFace.set(testFaces);
 
     forAll(allCellInfo, celli)
     {
         if (allCellInfo[celli].valid(wallDistCalc.data()))
         {
-            const point& cc = mesh_.cellCentres()[celli];
+            isProximityFace.set(mesh_.cells()[celli]);
+        }
+    }
 
-            const List<point>& origin = allCellInfo[celli].origin();
-            const List<FixedList<label, 3>>& surface =
-                allCellInfo[celli].surface();
+    syncTools::syncFaceList
+    (
+        mesh_,
+        isProximityFace,
+        orEqOp<unsigned int>(),
+        0u
+    );
 
-            // Find pair with minimum distance
-            for (label i = 0; i < origin.size(); i++)
+    testFaces = isProximityFace.toc();
+}
+
+
+Foam::label Foam::meshRefinement::markFakeGapRefinement
+(
+    const scalar planarCos,
+    const labelList& blockedSurfaces,
+    const label nAllowRefine,
+    const labelList& neiLevel,
+    const pointField& neiCc,
+
+    labelList& refineCell,
+    label& nRefine
+) const
+{
+    // Re-work the blockLevel of the blockedSurfaces into a length scale
+    // and a number of cells to cross
+    List<scalarList> regionToBlockSize(surfaces_.surfaces().size());
+    scalar maxBlockSize(-1);
+    for (const label surfi : blockedSurfaces)
+    {
+        const label geomi = surfaces_.surfaces()[surfi];
+        const searchableSurface& s = surfaces_.geometry()[geomi];
+        const label nRegions = s.regions().size();
+
+        regionToBlockSize[surfi].setSize(nRegions, -1);
+        for (label regioni = 0; regioni < nRegions; regioni++)
+        {
+            const label globalRegioni = surfaces_.globalRegion(surfi, regioni);
+            const label bLevel = surfaces_.blockLevel()[globalRegioni];
+
+            // If valid blockLevel specified
+            if (bLevel != labelMax)
             {
-                for (label j = i + 1; j < origin.size(); j++)
+                regionToBlockSize[surfi][regioni] =
+                    meshCutter_.level0EdgeLength()/pow(2.0, bLevel);
+            }
+
+            maxBlockSize = max(maxBlockSize, regionToBlockSize[surfi][regioni]);
+
+            //const label mLevel = surfaces_.maxLevel()[globalRegioni];
+            //// TBD: check for higher cached level of surface due to vol
+            ////      refinement. Problem: might still miss refinement bubble
+            ////      fully inside thin channel
+            //if (isA<triSurfaceMesh>(s) && !isA<distributedTriSurfaceMesh>(s))
+            //{
+            //    const triSurfaceMesh& surf = refCast<const triSurfaceMesh>(s);
+            //}
+        }
+    }
+
+    // Collect candidate faces (i.e. intersecting any surface and
+    // owner/neighbour not yet refined) and their cells
+    labelList cellMap;
+    {
+        labelList testFaces(getRefineCandidateFaces(refineCell));
+
+        // Extend the set by faceCellFace wave upto given 3*blockLevel size
+        markProximityCandidateFaces
+        (
+            blockedSurfaces,
+            regionToBlockSize,
+            neiLevel,
+            neiCc,
+            testFaces
+        );
+
+
+        bitSet isTestCell(mesh_.nCells());
+        for (const label facei : testFaces)
+        {
+            const label own = mesh_.faceOwner()[facei];
+            if (refineCell[own] == -1)
+            {
+                isTestCell.set(own);
+            }
+            if (mesh_.isInternalFace(facei))
+            {
+                const label nei = mesh_.faceNeighbour()[facei];
+                if (refineCell[nei] == -1)
                 {
-                    //if (isMultiRegion[celli])
-                    //{
-                    //    // The intersection locations are too inaccurate
-                    //    // (since not proper nearest, just a cell-cell ray
-                    //    //  intersection) so include always
-                    //
-                    //    smallGapDistance[celli] =
-                    //        max(smallGapDistance[celli], maxDist);
-                    //
-                    //
-                    //    str.writeLine(cc, origin[i]);
-                    //    str.writeLine(cc, origin[j]);
-                    //
-                    //    nMulti++;
-                    //}
-                    //else
-                    if (((cc-origin[i]) & (cc-origin[j])) < 0)
-                    {
-                        const label surfi = surface[i][0];
-                        const label regioni = surface[i][1];
+                    isTestCell.set(nei);
+                }
+            }
+        }
+        cellMap = isTestCell.sortedToc();
+    }
 
-                        const label surfj = surface[j][0];
-                        const label regionj = surface[j][1];
 
-                        const scalar maxSize = max
-                        (
-                            regionToBlockSize[surfi][regioni],
-                            regionToBlockSize[surfj][regionj]
-                        );
+    // Shoot rays in all 6 directions
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    //
+    // Almost exact copy of meshRefinement::markInternalGapRefinement.
+    // Difference is only in the length scale (now based on blockLevel)
 
-                        if
-                        (
-                            magSqr(origin[i]-origin[j])
-                          < Foam::sqr(2*maxSize)
-                        )
-                        {
-                            const scalar maxDist
-                            (
-                                max
-                                (
-                                    mag(cc-origin[i]),
-                                    mag(cc-origin[j])
-                                )
-                            );
 
-                            smallGapDistance[celli] =
-                                max(smallGapDistance[celli], maxDist);
-                            nSmallGap++;
-                        }
-                    }
+    // Find per cell centre the nearest point and normal on the surfaces
+    List<pointIndexHit> nearInfo;
+    vectorField nearNormal;
+    labelList nearSurface;
+    labelList nearRegion;
+    {
+        surfaces_.findNearestRegion
+        (
+            blockedSurfaces,
+            pointField(mesh_.cellCentres(), cellMap),
+            scalarField(cellMap.size(), maxBlockSize),  // use uniform for now
+            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())
+        {
+            const label globalRegioni = surfaces_.globalRegion
+            (
+                nearSurface[i],
+                nearRegion[i]
+            );
+
+            // Combine info from shell and surface
+            const label bLevel = surfaces_.blockLevel()[globalRegioni];
+            const FixedList<label, 3> gapInfo
+            ({
+                100,
+                bLevel,
+                100
+            });
+            const volumeType gapMode(volumeType::MIXED);
+
+            // Store wanted number of cells in gap
+            const label celli = cellMap[i];
+            const label cLevel = meshCutter_.cellLevel()[celli];
+
+            // Construct one or more rays to test for oppositeness
+            const label nRays = generateRays
+            (
+                false,                      // ignore actual surface normal
+                nearInfo[i].hitPoint(),
+                nearNormal[i],
+                gapInfo,
+                gapMode,
+
+                mesh_.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();
 
-    if (debug)
-    {
-        Info<< "Marked for blocking due to intersecting multiple surfaces  : "
-            << returnReduce(nMulti, sumOp<label>()) << " cells." << endl;
-        Info<< "Marked for blocking due to close opposite surfaces         : "
-            << returnReduce(nSmallGap, sumOp<label>()) << " cells." << endl;
-    }
+    rayStart2.shrink();
+    rayEnd2.shrink();
+    gapSize2.shrink();
 
-    if (debug&meshRefinement::MESH)
-    {
-        volScalarField distance
-        (
-            IOobject
-            (
-                "smallGapDistance",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                IOobject::NO_REGISTER
-            ),
-            mesh_,
-            dimensionedScalar(dimLength, Zero)
-        );
-        distance.field() = smallGapDistance;
-        distance.correctBoundaryConditions();
-
-        Info<< "Writing all small-gap cells to "
-            << distance.name() << endl;
-        distance.write();
-    }
+    cellMap = labelUIndList(cellMap, map)();
+    nearNormal = UIndirectList<vector>(nearNormal, map)();
+    nearInfo.clear();
+    nearSurface.clear();
+    nearRegion.clear();
 
 
-    // Mark refinement
+    // Do intersections in pairs
+    labelList surf1;
+    labelList region1;
+    List<pointIndexHit> hit1;
+    vectorField normal1;
+    surfaces_.findNearestIntersection
+    (
+        rayStart,
+        rayEnd,
+        surf1,
+        region1,
+        hit1,
+        normal1
+    );
+
+    labelList surf2;
+    labelList region2;
+    List<pointIndexHit> hit2;
+    vectorField normal2;
+    surfaces_.findNearestIntersection
+    (
+        rayStart2,
+        rayEnd2,
+        surf2,
+        region2,
+        hit2,
+        normal2
+    );
+
+    // Extract cell based gap size
     const label oldNRefine = nRefine;
-    forAll(smallGapDistance, celli)
+    forAll(surf1, i)
     {
-        if (smallGapDistance[celli] > SMALL)
+        if (surf1[i] != -1 && surf2[i] != -1)
         {
-            if
-            (
-                !markForRefine
-                (
-                    0,                      // mark level
-                    nAllowRefine,
-                    refineCell[celli],
-                    nRefine
-                )
-            )
+            // Found intersections with surface. Check for
+            // - small gap
+            // - coplanar normals
+
+            const label celli = cellMap[i];
+            if (celli != -1 && refineCell[celli] == -1)
             {
-                if (debug)
+                const scalar d2 = magSqr(hit1[i].hitPoint()-hit2[i].hitPoint());
+
+                const scalar maxGapSize
+                (
+                    max
+                    (
+                        regionToBlockSize[surf1[i]][region1[i]],
+                        regionToBlockSize[surf2[i]][region2[i]]
+                    )
+                );
+
+                if
+                (
+                    (mag(normal1[i]&normal2[i]) > planarCos)
+                 && (d2 < Foam::sqr(maxGapSize))
+                )
                 {
-                    Pout<< "Stopped refining since reaching my cell"
-                        << " limit of " << mesh_.nCells()+7*nRefine
-                        << endl;
+                    if
+                    (
+                        !markForRefine
+                        (
+                            0,                      // mark level
+                            nAllowRefine,
+                            refineCell[celli],
+                            nRefine
+                        )
+                    )
+                    {
+                        if (debug)
+                        {
+                            Pout<< "Stopped refining since reaching my cell"
+                                << " limit of " << mesh_.nCells()+7*nRefine
+                                << endl;
+                        }
+                        break;
+                    }
                 }
-                break;
             }
         }
     }
@@ -1006,7 +1076,20 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::removeGapCells
     }
 
 
-    markProximityRefinementWave
+    //markProximityRefinementWave
+    //(
+    //    Foam::cos(degToRad(planarAngle)),
+    //    surfToBlockLevel.sortedToc(),
+    //
+    //    labelMax/Pstream::nProcs(), //nAllowRefine,
+    //    neiLevel,
+    //    neiCc,
+    //
+    //    refineCell,
+    //    nRefine
+    //);
+
+    markFakeGapRefinement
     (
         Foam::cos(degToRad(planarAngle)),
         surfToBlockLevel.sortedToc(),
@@ -1020,20 +1103,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::removeGapCells
     );
 
 
-    //// Mark big-gap refinement
-    //markFakeGapRefinement
-    //(
-    //    Foam::cos(degToRad(planarAngle)),
-    //
-    //    labelMax/Pstream::nProcs(), //nAllowRefine,
-    //    neiLevel,
-    //    neiCc,
-    //
-    //    refineCell,
-    //    nRefine
-    //);
-
-
     Info<< "Marked for blocking due to close opposite surfaces         : "
         << returnReduce(nRefine, sumOp<label>()) << " cells." << endl;
 
@@ -2147,4 +2216,499 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::blockLeakFaces
 }
 
 
+
+/*
+Foam::label Foam::meshRefinement::markProximityRefinementWave
+(
+    const scalar planarCos,
+    const labelList& blockedSurfaces,
+    const label nAllowRefine,
+    const labelList& neiLevel,
+    const pointField& neiCc,
+
+    labelList& refineCell,
+    label& nRefine
+) const
+{
+    labelListList faceZones(surfaces_.surfaces().size());
+    {
+        // If triSurface do additional zoning based on connectivity
+        for (const label surfi : blockedSurfaces)
+        {
+            const label geomi = surfaces_.surfaces()[surfi];
+            const searchableSurface& s = surfaces_.geometry()[geomi];
+            if (isA<triSurfaceMesh>(s) && !isA<distributedTriSurfaceMesh>(s))
+            {
+                const triSurfaceMesh& surf = refCast<const triSurfaceMesh>(s);
+                const labelListList& edFaces = surf.edgeFaces();
+                boolList isOpenEdge(edFaces.size(), false);
+                forAll(edFaces, i)
+                {
+                    if (edFaces[i].size() == 1)
+                    {
+                        isOpenEdge[i] = true;
+                    }
+                }
+
+                labelList faceZone;
+                const label nZones = surf.markZones(isOpenEdge, faceZone);
+                if (nZones > 1)
+                {
+                    faceZones[surfi].transfer(faceZone);
+                }
+            }
+        }
+    }
+
+
+    // Re-work the blockLevel of the blockedSurfaces into a length scale
+    // and a number of cells to cross
+    List<scalarList> regionToBlockSize(surfaces_.surfaces().size());
+    for (const label surfi : blockedSurfaces)
+    {
+        const label geomi = surfaces_.surfaces()[surfi];
+        const searchableSurface& s = surfaces_.geometry()[geomi];
+        const label nRegions = s.regions().size();
+        regionToBlockSize[surfi].setSize(nRegions);
+        for (label regioni = 0; regioni < nRegions; regioni++)
+        {
+            const label globalRegioni = surfaces_.globalRegion(surfi, regioni);
+            const label bLevel = surfaces_.blockLevel()[globalRegioni];
+            regionToBlockSize[surfi][regioni] =
+                meshCutter_.level0EdgeLength()/pow(2.0, bLevel);
+
+            //const label mLevel = surfaces_.maxLevel()[globalRegioni];
+            //// TBD: check for higher cached level of surface due to vol
+            ////      refinement. Problem: might still miss refinement bubble
+            ////      fully inside thin channel
+            //if (isA<triSurfaceMesh>(s) && !isA<distributedTriSurfaceMesh>(s))
+            //{
+            //    const triSurfaceMesh& surf = refCast<const triSurfaceMesh>(s);
+            //}
+
+            //nIters = max(nIters, (2<<(mLevel-bLevel)));
+        }
+    }
+
+    // Clever limiting of the number of iterations (= max cells in the channel)
+    // since it has too many problematic issues, e.g. with volume refinement
+    // and the real check uses the proper distance anyway just disable.
+    const label nIters = mesh_.globalData().nTotalCells();
+
+
+    // Collect candidate faces (i.e. intersecting any surface and
+    // owner/neighbour not yet refined)
+    const 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
+    );
+    // TBD. correct nIters for actual cellLevel (since e.g. volume refinement
+    //      might add to cell level). Unfortunately we only have minLevel,
+    //      not maxLevel. Problem: what if volume refinement only in middle of
+    //      channel? Say channel 1m wide with a 0.1m sphere of refinement
+    //      Workaround: have dummy surface with e.g. maxLevel 100 to
+    //      force nIters to be high enough.
+
+
+    // Test for all intersections (with surfaces of higher gap level than
+    // minLevel) and cache per cell the max surface level and the local normal
+    // on that surface.
+
+    labelList surface1;
+    List<pointIndexHit> hit1;
+    labelList region1;
+    vectorField normal1;
+
+    labelList surface2;
+    List<pointIndexHit> hit2;
+    labelList region2;
+    vectorField normal2;
+
+    surfaces_.findNearestIntersection
+    (
+        blockedSurfaces,
+        start,
+        end,
+
+        surface1,
+        hit1,
+        region1,    // local region
+        normal1,
+
+        surface2,
+        hit2,
+        region2,    // local region
+        normal2
+    );
+
+
+    // Detect cells that are using multiple surface regions
+    //bitSet isMultiRegion;
+    //detectMultiRegionCells
+    //(
+    //    faceZones,
+    //    testFaces,
+    //
+    //    surface1,
+    //    hit1,
+    //    region1,
+    //
+    //    surface2,
+    //    hit2,
+    //    region2,
+    //
+    //    isMultiRegion
+    //);
+
+
+    label n = 0;
+    forAll(testFaces, i)
+    {
+        if (hit1[i].hit())
+        {
+            n++;
+        }
+    }
+
+    List<wallPoints> faceDist(n);
+    labelList changedFaces(n);
+    n = 0;
+
+    DynamicList<point> originLocation(2);
+    DynamicList<scalar> originDistSqr(2);
+    DynamicList<FixedList<label, 3>> originSurface(2);
+    //DynamicList<point> originNormal(2);
+
+
+    //- To avoid walking through surfaces we mark all faces that have been
+    //  intersected. We can either mark only those faces intersecting
+    //  blockedSurfaces (i.e. with a 'blockLevel') or mark faces intersecting
+    //  any (refinement) surface (this includes e.g. faceZones). This is
+    //  much easier since that information is already cached
+    //  (meshRefinement::intersectedFaces())
+
+    //bitSet isBlockedFace(mesh_.nFaces());
+    forAll(testFaces, i)
+    {
+        if (hit1[i].hit())
+        {
+            const label facei = testFaces[i];
+            //isBlockedFace.set(facei);
+            const point& fc = mesh_.faceCentres()[facei];
+            const labelList& fz1 = faceZones[surface1[i]];
+
+            originLocation.clear();
+            originDistSqr.clear();
+            //originNormal.clear();
+            originSurface.clear();
+
+            originLocation.append(hit1[i].point());
+            originDistSqr.append(hit1[i].point().distSqr(fc));
+            //originNormal.append(normal1[i]);
+            originSurface.append
+            (
+                FixedList<label, 3>
+                ({
+                    surface1[i],
+                    region1[i],
+                    (fz1.size() ? fz1[hit1[i].index()] : region1[i])
+                })
+            );
+
+            if (hit2[i].hit() && hit1[i] != hit2[i])
+            {
+                const labelList& fz2 = faceZones[surface2[i]];
+                originLocation.append(hit2[i].point());
+                originDistSqr.append(hit2[i].point().distSqr(fc));
+                //originNormal.append(normal2[i]);
+                originSurface.append
+                (
+                    FixedList<label, 3>
+                    ({
+                        surface2[i],
+                        region2[i],
+                        (fz2.size() ? fz2[hit2[i].index()] : region2[i])
+                    })
+                );
+            }
+
+            // Collect all seed data. Currently walking does not look at
+            // surface direction - if so pass in surface normal as well
+            faceDist[n] = wallPoints
+            (
+                originLocation,     // origin
+                originDistSqr,      // distance to origin
+                originSurface       // surface+region+zone
+                //originNormal        // normal at origin
+            );
+            changedFaces[n] = facei;
+            n++;
+        }
+    }
+
+
+    // Clear intersection info
+    surface1.clear();
+    hit1.clear();
+    region1.clear();
+    normal1.clear();
+    surface2.clear();
+    hit2.clear();
+    region2.clear();
+    normal2.clear();
+
+
+    List<wallPoints> allFaceInfo(mesh_.nFaces());
+    List<wallPoints> allCellInfo(mesh_.nCells());
+
+    // Any refinement surface (even a faceZone) should stop the gap walking.
+    // This is exactly the information which is cached in the surfaceIndex_
+    // field.
+    const bitSet isBlockedFace(intersectedFaces());
+
+    wallPoints::trackData td(isBlockedFace, regionToBlockSize);
+    FaceCellWave<wallPoints, wallPoints::trackData> wallDistCalc
+    (
+        mesh_,
+        changedFaces,
+        faceDist,
+        allFaceInfo,
+        allCellInfo,
+        0,            // max iterations
+        td
+    );
+    wallDistCalc.iterate(nIters);
+
+
+    if (debug&meshRefinement::MESH)
+    {
+        // Dump current nearest opposite surfaces
+        volScalarField distance
+        (
+            IOobject
+            (
+                "gapSize",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh_,
+            dimensionedScalar
+            (
+                "zero",
+                dimLength,  //dimensionSet(0, 1, 0, 0, 0),
+                -1
+            )
+        );
+
+        forAll(allCellInfo, celli)
+        {
+            if (allCellInfo[celli].valid(wallDistCalc.data()))
+            {
+                const point& cc = mesh_.cellCentres()[celli];
+                // Nearest surface points
+                const List<point>& origin = allCellInfo[celli].origin();
+
+                // Find 'opposite' pair with minimum distance
+                for (label i = 0; i < origin.size(); i++)
+                {
+                    for (label j = i + 1; j < origin.size(); j++)
+                    {
+                        if (((cc-origin[i]) & (cc-origin[j])) < 0)
+                        {
+                            const scalar d(mag(origin[i]-origin[j]));
+                            if (distance[celli] < 0)
+                            {
+                                distance[celli] = d;
+                            }
+                            else
+                            {
+                                distance[celli] = min(distance[celli], d);
+                            }
+                        }
+                    }
+                }
+            }
+        }
+        distance.correctBoundaryConditions();
+
+        Info<< "Writing measured gap distance to "
+            << distance.name() << endl;
+        distance.write();
+    }
+
+
+
+    // Detect tight gaps:
+    // - cell is inbetween the two surfaces
+    // - two surfaces are planarish
+    // - two surfaces are not too far apart
+    //   (number of walking iterations is a too-coarse measure)
+
+    scalarField smallGapDistance(mesh_.nCells(), 0.0);
+    label nMulti = 0;
+    label nSmallGap = 0;
+
+    //OBJstream str(mesh_.time().timePath()/"multiRegion.obj");
+
+
+    forAll(allCellInfo, celli)
+    {
+        if (allCellInfo[celli].valid(wallDistCalc.data()))
+        {
+            const point& cc = mesh_.cellCentres()[celli];
+
+            const List<point>& origin = allCellInfo[celli].origin();
+            const List<FixedList<label, 3>>& surface =
+                allCellInfo[celli].surface();
+
+            // Find pair with minimum distance
+            for (label i = 0; i < origin.size(); i++)
+            {
+                for (label j = i + 1; j < origin.size(); j++)
+                {
+                    //if (isMultiRegion[celli])
+                    //{
+                    //    // The intersection locations are too inaccurate
+                    //    // (since not proper nearest, just a cell-cell ray
+                    //    //  intersection) so include always
+                    //
+                    //    smallGapDistance[celli] =
+                    //        max(smallGapDistance[celli], maxDist);
+                    //
+                    //
+                    //    str.writeLine(cc, origin[i]);
+                    //    str.writeLine(cc, origin[j]);
+                    //
+                    //    nMulti++;
+                    //}
+                    //else
+                    if (((cc-origin[i]) & (cc-origin[j])) < 0)
+                    {
+                        const label surfi = surface[i][0];
+                        const label regioni = surface[i][1];
+
+                        const label surfj = surface[j][0];
+                        const label regionj = surface[j][1];
+
+                        const scalar maxSize = max
+                        (
+                            regionToBlockSize[surfi][regioni],
+                            regionToBlockSize[surfj][regionj]
+                        );
+
+                        if
+                        (
+                            magSqr(origin[i]-origin[j])
+                          < Foam::sqr(2*maxSize)
+                        )
+                        {
+                            const scalar maxDist
+                            (
+                                max
+                                (
+                                    mag(cc-origin[i]),
+                                    mag(cc-origin[j])
+                                )
+                            );
+
+                            smallGapDistance[celli] =
+                                max(smallGapDistance[celli], maxDist);
+                            nSmallGap++;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+
+    if (debug)
+    {
+        Info<< "Marked for blocking due to intersecting multiple surfaces  : "
+            << returnReduce(nMulti, sumOp<label>()) << " cells." << endl;
+        Info<< "Marked for blocking due to close opposite surfaces         : "
+            << returnReduce(nSmallGap, sumOp<label>()) << " cells." << endl;
+    }
+
+    if (debug&meshRefinement::MESH)
+    {
+        volScalarField distance
+        (
+            IOobject
+            (
+                "smallGapDistance",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh_,
+            dimensionedScalar(dimLength, Zero)
+        );
+        distance.field() = smallGapDistance;
+        distance.correctBoundaryConditions();
+
+        Info<< "Writing all small-gap cells to "
+            << distance.name() << endl;
+        distance.write();
+    }
+
+
+    // Mark refinement
+    const label oldNRefine = nRefine;
+    forAll(smallGapDistance, celli)
+    {
+        if (smallGapDistance[celli] > SMALL)
+        {
+            if
+            (
+                !markForRefine
+                (
+                    0,                      // mark level
+                    nAllowRefine,
+                    refineCell[celli],
+                    nRefine
+                )
+            )
+            {
+                if (debug)
+                {
+                    Pout<< "Stopped refining since reaching my cell"
+                        << " limit of " << mesh_.nCells()+7*nRefine
+                        << endl;
+                }
+                break;
+            }
+        }
+    }
+
+    if
+    (
+        returnReduce(nRefine, sumOp<label>())
+      > returnReduce(nAllowRefine, sumOp<label>())
+    )
+    {
+        Info<< "Reached refinement limit." << endl;
+    }
+
+    return returnReduce(nRefine-oldNRefine, sumOp<label>());
+}
+*/
+
 // ************************************************************************* //
diff --git a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C
index 81b70f72e9..5b32ac72cf 100644
--- a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -1954,6 +1954,59 @@ void Foam::refinementSurfaces::findNearestIntersection
 }
 
 
+void Foam::refinementSurfaces::findNearestIntersection
+(
+//    const labelList& surfacesToTest,
+    const pointField& start,
+    const pointField& end,
+
+    labelList& surface1,
+    labelList& region1,
+    List<pointIndexHit>& hitInfo1,
+    vectorField& normal1
+) const
+{
+    // Initialize arguments
+    surface1.setSize(start.size());
+    surface1 = -1;
+    region1.setSize(start.size());
+    region1 = -1;
+    hitInfo1.setSize(start.size());
+    hitInfo1 = pointIndexHit();
+    normal1.setSize(start.size());
+    normal1 = 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.getRegion(nearestInfo, region);
+        geom.getNormal(nearestInfo, normal);
+
+        forAll(nearestInfo, pointI)
+        {
+            if (nearestInfo[pointI].hit())
+            {
+                surface1[pointI] = surfI;
+                region1[pointI] = region[pointI];
+                hitInfo1[pointI] = nearestInfo[pointI];
+                normal1[pointI] = normal[pointI];
+                nearest[pointI] = nearestInfo[pointI].point();
+            }
+        }
+    }
+}
+
+
 void Foam::refinementSurfaces::findNearestIntersection
 (
     const pointField& start,
diff --git a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H
index 2c4ff38380..8f02203358 100644
--- a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -437,6 +437,17 @@ public:
                 vectorField& normal
             ) const;
 
+            //- Find nearest (to start only) intersection of edge
+            void findNearestIntersection
+            (
+                const pointField& start,
+                const pointField& end,
+                labelList& surfaces,
+                labelList& regions,
+                List<pointIndexHit>&,
+                vectorField& normal
+            ) const;
+
             //- Find nearest (to start only) intersection of edge
             void findNearestIntersection
             (