ENH: snappyHexMesh: proximity check

This adds automatic deletion of cells inside small gaps. This is
generally used to avoid having excessive numbers of cells in irrelevant
areas of a geometry. It is nearly the opposite of automatic gap refinement
 - that refines cells to resolve the gap; this functionality removes cells
to not mesh the gap.

The proximity handling will remove those cells which are inside 'thin' gaps
where 'thin' is defined as a distance of 2*'blockLevel'
It will
- detect surfaces which have the new 'blockLevel' specification
- convert this to a minimum gap distance
- detect cells which are inside this gap
- remove these cells and add exposed faces to the nearest 'real' patch

src/mesh/snappyHexMesh/Make/files

@@ -15,6 +15,7 @@ meshRefinement/meshRefinementProblemCells.C
 meshRefinement/meshRefinementRefine.C
 meshRefinement/meshRefinementGapRefine.C
 meshRefinement/meshRefinementBlock.C
+meshRefinement/wallPoints.C
 meshRefinement/patchFaceOrientation.C
 
 refinementFeatures/refinementFeatures.C

src/mesh/snappyHexMesh/Make/options

@@ -7,7 +7,8 @@ EXE_INC = \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/overset/lnInclude \
-    -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude
+    -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
+    -I$(LIB_SRC)/parallel/distributed/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \

src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H

@@ -533,6 +533,47 @@ private:
                 label& nRefine
             ) const;
 
+            void markMultiRegionCell
+            (
+                const label celli,
+                //const label globalRegion,
+                //const label zonei,
+                const FixedList<label, 3>& surface,
+
+                Map<FixedList<label, 3>>& cellToRegions,
+                bitSet& isMultiRegion
+            ) const;
+
+            void detectMultiRegionCells
+            (
+                const labelListList& faceZones,
+                const labelList& testFaces,
+
+                const labelList& surface1,
+                const List<pointIndexHit>& hit1,
+                const labelList& region1,
+
+                const labelList& surface2,
+                const List<pointIndexHit>& hit2,
+                const labelList& region2,
+
+                bitSet& isMultiRegion
+            ) const;
+
+            //- Mark cells for surface proximity based refinement.
+            label markProximityRefinementWave
+            (
+                const scalar planarCos,
+                const labelList& blockedSurfaces,
+                const label nAllowRefine,
+                const labelList& neiLevel,
+                const pointField& neiCc,
+
+                labelList& refineCell,
+                label& nRefine
+            ) const;
+
+
         // Baffle handling
 
             //- Get faces to repatch. Returns map from face to patch.

src/mesh/snappyHexMesh/meshRefinement/meshRefinementBlock.C

@@ -33,6 +33,12 @@ License
 #include "unitConversion.H"
 #include "bitSet.H"
 
+#include "FaceCellWave.H"
+#include "volFields.H"
+#include "wallPoints.H"
+#include "searchableSurfaces.H"
+#include "distributedTriSurfaceMesh.H"
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 //Foam::label Foam::meshRefinement::markFakeGapRefinement
@@ -288,6 +294,636 @@ void Foam::meshRefinement::growSet
 }
 
 
+void Foam::meshRefinement::markMultiRegionCell
+(
+    const label celli,
+    const FixedList<label, 3>& surface,
+
+    Map<FixedList<label, 3>>& cellToRegions,
+    bitSet& isMultiRegion
+) const
+{
+    if (!isMultiRegion[celli])
+    {
+        Map<FixedList<label, 3>>::iterator fnd = cellToRegions.find(celli);
+        if (!fnd.found())
+        {
+            cellToRegions.insert(celli, surface);
+        }
+        else if (fnd() != surface)
+        {
+            // Found multiple intersections on cell
+            isMultiRegion.set(celli);
+        }
+    }
+}
+
+
+void Foam::meshRefinement::detectMultiRegionCells
+(
+    const labelListList& faceZones,
+    const labelList& testFaces,
+
+    const labelList& surface1,
+    const List<pointIndexHit>& hit1,
+    const labelList& region1,
+
+    const labelList& surface2,
+    const List<pointIndexHit>& hit2,
+    const labelList& region2,
+
+    bitSet& isMultiRegion
+) const
+{
+    isMultiRegion.clear();
+    isMultiRegion.setSize(mesh_.nCells());
+
+    Map<FixedList<label, 3>> cellToRegions(testFaces.size());
+
+    forAll(testFaces, i)
+    {
+        const pointIndexHit& info1 = hit1[i];
+        if (info1.hit())
+        {
+            const label facei = testFaces[i];
+            const labelList& fz1 = faceZones[surface1[i]];
+            const FixedList<label, 3> surfaceInfo1
+            ({
+                surface1[i],
+                region1[i],
+                (fz1.size() ? fz1[info1.index()] : region1[i])
+            });
+
+            markMultiRegionCell
+            (
+                mesh_.faceOwner()[facei],
+                surfaceInfo1,
+                cellToRegions,
+                isMultiRegion
+            );
+            if (mesh_.isInternalFace(facei))
+            {
+                markMultiRegionCell
+                (
+                    mesh_.faceNeighbour()[facei],
+                    surfaceInfo1,
+                    cellToRegions,
+                    isMultiRegion
+                );
+            }
+
+            const pointIndexHit& info2 = hit2[i];
+
+            if (info2.hit() && info1 != info2)
+            {
+                const labelList& fz2 = faceZones[surface2[i]];
+                const FixedList<label, 3> surfaceInfo2
+                ({
+                    surface2[i],
+                    region2[i],
+                    (fz2.size() ? fz2[info2.index()] : region2[i])
+                });
+
+                markMultiRegionCell
+                (
+                    mesh_.faceOwner()[facei],
+                    surfaceInfo2,
+                    cellToRegions,
+                    isMultiRegion
+                );
+                if (mesh_.isInternalFace(facei))
+                {
+                    markMultiRegionCell
+                    (
+                        mesh_.faceNeighbour()[facei],
+                        surfaceInfo2,
+                        cellToRegions,
+                        isMultiRegion
+                    );
+                }
+            }
+        }
+    }
+
+
+    if (debug&meshRefinement::MESH)
+    {
+        volScalarField multiCell
+        (
+            IOobject
+            (
+                "multiCell",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar
+            (
+                "zero",
+                dimensionSet(0, 1, 0, 0, 0),
+                0.0
+            )
+        );
+        forAll(isMultiRegion, celli)
+        {
+            if (isMultiRegion[celli])
+            {
+                multiCell[celli] = 1.0;
+            }
+        }
+
+        Info<< "Writing all multi cells to " << multiCell.name() << endl;
+        multiCell.write();
+    }
+}
+
+
+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
+    scalarField regionToBlockSize(surfaces_.blockLevel().size(), 0);
+    //label nIters = 2;
+
+    for (const label surfi : blockedSurfaces)
+    {
+        const label geomi = surfaces_.surfaces()[surfi];
+        const searchableSurface& s = surfaces_.geometry()[geomi];
+        const label nRegions = s.regions().size();
+        for (label regioni = 0; regioni < nRegions; regioni++)
+        {
+            const label globalRegioni = surfaces_.globalRegion(surfi, regioni);
+            const label bLevel = surfaces_.blockLevel()[globalRegioni];
+            regionToBlockSize[globalRegioni] =
+                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.
+    // Since 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? 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);
+
+    forAll(testFaces, i)
+    {
+        if (hit1[i].hit())
+        {
+            const label facei = testFaces[i];
+            const point& fc = mesh_.faceCentres()[facei];
+            const labelList& fz1 = faceZones[surface1[i]];
+
+            originLocation.clear();
+            originDistSqr.clear();
+            //originNormal.clear();
+            originSurface.clear();
+
+            originLocation.append(hit1[i].hitPoint());
+            originDistSqr.append(magSqr(fc-hit1[i].hitPoint()));
+            //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].hitPoint());
+                originDistSqr.append(magSqr(fc-hit2[i].hitPoint()));
+                //originNormal.append(normal2[i]);
+                originSurface.append
+                (
+                    FixedList<label, 3>
+                    ({
+                        surface2[i],
+                        region2[i],
+                        (fz2.size() ? fz2[hit2[i].index()] : region2[i])
+                    })
+                );
+            }
+
+            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());
+
+    FaceCellWave<wallPoints> wallDistCalc
+    (
+        mesh_,
+        changedFaces,
+        faceDist,
+        allFaceInfo,
+        allCellInfo,
+        0             // max iterations
+    );
+    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,
+                false
+            ),
+            mesh_,
+            dimensionedScalar
+            (
+                "zero",
+                dimLength,  //dimensionSet(0, 1, 0, 0, 0),
+                -1
+            )
+        );
+
+        forAll(allCellInfo, celli)
+        {
+            if (allCellInfo[celli].valid(wallDistCalc.data()))
+            {
+                const point& cc = mesh_.C()[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_.C()[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++)
+                {
+                    scalar maxDist
+                    (
+                        max
+                        (
+                            mag(cc-origin[i]),
+                            mag(cc-origin[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.write(linePointRef(cc, origin[i]));
+                    //    str.write(linePointRef(cc, origin[j]));
+                    //
+                    //    nMulti++;
+                    //}
+                    //else
+                    if (((cc-origin[i]) & (cc-origin[j])) < 0)
+                    {
+                        const label globalRegioni = surfaces_.globalRegion
+                        (
+                            surface[i][0],
+                            surface[i][1]
+                        );
+                        const label globalRegionj = surfaces_.globalRegion
+                        (
+                            surface[j][0],
+                            surface[j][1]
+                        );
+
+                        const scalar maxSize = max
+                        (
+                            regionToBlockSize[globalRegioni],
+                            regionToBlockSize[globalRegionj]
+                        );
+
+                        if
+                        (
+                            magSqr(origin[i]-origin[j])
+                          < Foam::sqr(2*maxSize)
+                        )
+                        {
+                            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,
+                false
+            ),
+            mesh_,
+            dimensionedScalar
+            (
+                "zero",
+                dimensionSet(0, 1, 0, 0, 0),
+                0.0
+            )
+        );
+        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>());
+}
+
+
 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::removeGapCells
 (
     const scalar planarAngle,
@@ -303,12 +939,53 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::removeGapCells
 
     labelList refineCell(mesh_.nCells(), -1);
     label nRefine = 0;
-    markProximityRefinement
+    //markProximityRefinement
+    //(
+    //    Foam::cos(degToRad(planarAngle)),
+    //
+    //    minSurfaceLevel,                                // surface min level
+    //    labelList(minSurfaceLevel.size(), labelMax),    // surfaceGapLevel
+    //
+    //    labelMax/Pstream::nProcs(), //nAllowRefine,
+    //    neiLevel,
+    //    neiCc,
+    //
+    //    refineCell,
+    //    nRefine
+    //);
+
+
+    // Determine minimum blockLevel per surface
+    Map<label> surfToBlockLevel;
+
+    forAll(surfaces_.surfaces(), surfi)
+    {
+        const label geomi = surfaces_.surfaces()[surfi];
+        const searchableSurface& s = surfaces_.geometry()[geomi];
+        const label nRegions = s.regions().size();
+
+        label minBlockLevel = labelMax;
+        for (label regioni = 0; regioni < nRegions; regioni++)
+        {
+            const label globalRegioni = surfaces_.globalRegion(surfi, regioni);
+            minBlockLevel = min
+            (
+                minBlockLevel,
+                surfaces_.blockLevel()[globalRegioni]
+            );
+        }
+
+        if (minBlockLevel < labelMax)
+        {
+            surfToBlockLevel.insert(surfi, minBlockLevel);
+        }
+    }
+
+
+    markProximityRefinementWave
     (
         Foam::cos(degToRad(planarAngle)),
-
+        surfToBlockLevel.sortedToc(),
-        minSurfaceLevel,                                // surface min level
-        labelList(minSurfaceLevel.size(), labelMax),    // surfaceGapLevel
 
         labelMax/Pstream::nProcs(), //nAllowRefine,
         neiLevel,

src/mesh/snappyHexMesh/meshRefinement/wallPoints.C

@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2018-2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "wallPoints.H"
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Foam::Ostream& os,
+    const Foam::wallPoints& wDist
+)
+{
+    return os
+        << wDist.origin_ << token::SPACE << wDist.distSqr_
+        << token::SPACE << wDist.surface_;  // << token::SPACE << wDist.normal_;
+}
+
+
+Foam::Istream& Foam::operator>>
+(
+    Foam::Istream& is,
+    Foam::wallPoints& wDist
+)
+{
+    return is
+        >> wDist.origin_ >> wDist.distSqr_>> wDist.surface_;
+        //>> wDist.normal_;
+}
+
+
+// ************************************************************************* //

src/mesh/snappyHexMesh/meshRefinement/wallPoints.H

@@ -0,0 +1,260 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2018-2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::wallPoints
+
+Description
+    For use with FaceCellWave. Determines topological distance to starting faces
+
+SourceFiles
+    wallPointsI.H
+    wallPoints.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wallPoints_H
+#define wallPoints_H
+
+#include "point.H"
+#include "tensor.H"
+#include "DynamicList.H"
+#include "labelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyPatch;
+class polyMesh;
+
+// Forward declaration of friend functions and operators
+
+class wallPoints;
+
+Istream& operator>>(Istream&, wallPoints&);
+Ostream& operator<<(Ostream&, const wallPoints&);
+
+
+/*---------------------------------------------------------------------------*\
+                           Class wallPoints Declaration
+\*---------------------------------------------------------------------------*/
+
+class wallPoints
+{
+protected:
+
+    // Protected data
+
+        //- Starting points
+        DynamicList<point> origin_;
+
+        //- Distance (squared) from cellcenter to origin
+        DynamicList<scalar> distSqr_;
+
+        //- Originating surface,region and topological region
+        DynamicList<FixedList<label, 3>> surface_;
+
+        //- Originating normal
+        //DynamicList<vector> normal_;
+
+
+    // Protected 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& pt,
+            const label index1,
+            const wallPoints& w2,
+            const label index2,
+
+            const scalar tol,
+            TrackingData& td
+        );
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline wallPoints();
+
+        //- Construct from count
+        inline wallPoints
+        (
+            const UList<point>& origin,
+            const UList<scalar>& distSqr,
+            const UList<FixedList<label, 3>>& surface
+            //const UList<vector>& normal
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            inline const List<point>& origin() const
+            {
+                return origin_;
+            }
+            inline const List<scalar>& distSqr() const
+            {
+                return distSqr_;
+            }
+
+            inline const List<FixedList<label, 3>>& surface() const
+            {
+                return surface_;
+            }
+
+            //inline const List<vector>& normal() const
+            //{
+            //    return normal_;
+            //}
+
+
+        // Needed by FaceCellWave
+
+            //- Check whether origin has been changed at all or
+            //  still contains original (invalid) value.
+            template<class TrackingData>
+            inline bool valid(TrackingData& td) const;
+
+            //- Check for identical geometrical data. Used for cyclics checking.
+            template<class TrackingData>
+            inline bool sameGeometry
+            (
+                const polyMesh&,
+                const wallPoints&,
+                const scalar,
+                TrackingData& td
+            ) const;
+
+            //- Convert any absolute coordinates into relative to (patch)face
+            //  centre
+            template<class TrackingData>
+            inline void leaveDomain
+            (
+                const polyMesh&,
+                const polyPatch&,
+                const label patchFacei,
+                const point& faceCentre,
+                TrackingData& td
+            );
+
+            //- Reverse of leaveDomain
+            template<class TrackingData>
+            inline void enterDomain
+            (
+                const polyMesh&,
+                const polyPatch&,
+                const label patchFacei,
+                const point& faceCentre,
+                TrackingData& td
+            );
+
+            //- Apply rotation matrix to any coordinates
+            template<class TrackingData>
+            inline void transform
+            (
+                const polyMesh&,
+                const tensor&,
+                TrackingData& td
+            );
+
+            //- Influence of neighbouring face.
+            template<class TrackingData>
+            inline bool updateCell
+            (
+                const polyMesh&,
+                const label thisCelli,
+                const label neighbourFacei,
+                const wallPoints& neighbourInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of neighbouring cell.
+            template<class TrackingData>
+            inline bool updateFace
+            (
+                const polyMesh&,
+                const label thisFacei,
+                const label neighbourCelli,
+                const wallPoints& neighbourInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of different value on same face.
+            template<class TrackingData>
+            inline bool updateFace
+            (
+                const polyMesh&,
+                const label thisFacei,
+                const wallPoints& neighbourInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Same (like operator==)
+            template<class TrackingData>
+            inline bool equal(const wallPoints&, TrackingData&) const;
+
+
+    // Member Operators
+
+        // Needed for List IO
+        inline bool operator==(const wallPoints&) const;
+
+        inline bool operator!=(const wallPoints&) const;
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const wallPoints&);
+        friend Istream& operator>>(Istream&, wallPoints&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "wallPointsI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/mesh/snappyHexMesh/meshRefinement/wallPointsI.H

@@ -0,0 +1,352 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2018-2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "polyMesh.H"
+#include "transform.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class TrackingData>
+inline bool Foam::wallPoints::update
+(
+    const point& pt,
+    const label index1,
+    const wallPoints& w2,
+    const label index2,
+
+    const scalar tol,
+    TrackingData& td
+)
+{
+    scalar dist2 = magSqr(pt - w2.origin_[index2]);
+
+    if (!valid(td))
+    {
+        // currently not yet set so use any value
+        distSqr_[index1] = dist2;
+        origin_[index1] = w2.origin_[index2];
+        surface_[index1] = w2.surface_[index2];
+        //normal_[index1] = w2.normal_[index2];
+
+        return true;
+    }
+
+    scalar diff = distSqr_[index1] - dist2;
+
+    if (diff < 0)
+    {
+        // already nearer to pt
+        return false;
+    }
+
+    if
+    (
+        (diff < SMALL)
+     || ((distSqr_[index1] > SMALL) && (diff/distSqr_[index1] < tol))
+    )
+    {
+        // don't propagate small changes
+        return false;
+    }
+    else
+    {
+        // update with new values
+        distSqr_[index1] = dist2;
+        origin_[index1] = w2.origin_[index2];
+        surface_[index1] = w2.surface_[index2];
+        //normal_[index1] = w2.normal_[index2];
+
+        return true;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::wallPoints::wallPoints()
+:
+    origin_(0),
+    distSqr_(0),
+    surface_(0)
+    //normal_(0)
+{}
+
+
+inline Foam::wallPoints::wallPoints
+(
+    const UList<point>& origin,
+    const UList<scalar>& distSqr,
+    const UList<FixedList<label, 3>>& surface
+    //const UList<vector>& normal
+)
+:
+    origin_(origin),
+    distSqr_(distSqr),
+    surface_(surface)
+    //normal_(normal)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class TrackingData>
+inline bool Foam::wallPoints::valid(TrackingData& td) const
+{
+    return origin_.size();
+}
+
+
+// No geometric data so never any problem on cyclics
+template<class TrackingData>
+inline bool Foam::wallPoints::sameGeometry
+(
+    const polyMesh&,
+    const wallPoints&,
+    const scalar,
+    TrackingData&
+) const
+{
+    return true;
+}
+
+
+// No geometric data.
+template<class TrackingData>
+inline void Foam::wallPoints::leaveDomain
+(
+    const polyMesh&,
+    const polyPatch& patch,
+    const label patchFacei,
+    const point& faceCentre,
+    TrackingData&
+)
+{
+    for (auto& o : origin_)
+    {
+        o -= faceCentre;
+    }
+}
+
+
+// No geometric data.
+template<class TrackingData>
+inline void Foam::wallPoints::transform
+(
+    const polyMesh&,
+    const tensor& rotTensor,
+    TrackingData&
+)
+{
+    for (auto& o : origin_)
+    {
+        o = Foam::transform(rotTensor, o);
+    }
+}
+
+
+// No geometric data.
+template<class TrackingData>
+inline void Foam::wallPoints::enterDomain
+(
+    const polyMesh&,
+    const polyPatch& patch,
+    const label patchFacei,
+    const point& faceCentre,
+    TrackingData&
+)
+{
+    // back to absolute form
+    for (auto& o : origin_)
+    {
+        o += faceCentre;
+    }
+}
+
+
+// Update cell with neighbouring face information
+template<class TrackingData>
+inline bool Foam::wallPoints::updateCell
+(
+    const polyMesh& mesh,
+    const label thisCelli,
+    const label neighbourFacei,
+    const wallPoints& neighbourInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    const point& cc = mesh.cellCentres()[thisCelli];
+
+    bool hasChanged = false;
+    forAll(neighbourInfo.surface_, i)
+    {
+        const FixedList<label, 3>& nbrSurface = neighbourInfo.surface_[i];
+
+        // Find in my surfaces
+        label index = surface_.find(nbrSurface);
+        if (index == -1)
+        {
+            // Append
+            origin_.append(neighbourInfo.origin_[i]);
+            distSqr_.append(magSqr(cc-neighbourInfo.origin_[i]));
+            surface_.append(nbrSurface);
+            //normal_.append(neighbourInfo.normal_[i]);
+            hasChanged = true;
+        }
+        else
+        {
+            hasChanged =
+                update(cc, index, neighbourInfo, i, tol, td)
+             || hasChanged;
+        }
+    }
+
+    return hasChanged;
+}
+
+
+// Update face with neighbouring cell information
+template<class TrackingData>
+inline bool Foam::wallPoints::updateFace
+(
+    const polyMesh& mesh,
+    const label thisFacei,
+    const label neighbourCelli,
+    const wallPoints& neighbourInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    const point& fc = mesh.faceCentres()[thisFacei];
+
+    // From cell to its faces.
+    bool hasChanged = false;
+    forAll(neighbourInfo.surface_, i)
+    {
+        const FixedList<label, 3>& nbrSurface = neighbourInfo.surface_[i];
+
+        // Find in my surfaces
+        label index = surface_.find(nbrSurface);
+        if (index == -1)
+        {
+            // Append
+            origin_.append(neighbourInfo.origin_[i]);
+            distSqr_.append(magSqr(fc-neighbourInfo.origin_[i]));
+            surface_.append(nbrSurface);
+            //normal_.append(neighbourInfo.normal_[i]);
+            hasChanged = true;
+        }
+        else
+        {
+            hasChanged =
+                update(fc, index, neighbourInfo, i, tol, td)
+             || hasChanged;
+        }
+    }
+
+    return hasChanged;
+}
+
+
+// Update face with coupled face information
+template<class TrackingData>
+inline bool Foam::wallPoints::updateFace
+(
+    const polyMesh& mesh,
+    const label thisFacei,
+    const wallPoints& neighbourInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    const point& fc = mesh.faceCentres()[thisFacei];
+
+    // From face to face (e.g. coupled faces)
+    bool hasChanged = false;
+    forAll(neighbourInfo.surface_, i)
+    {
+        const FixedList<label, 3>& nbrSurface = neighbourInfo.surface_[i];
+
+        // Find in my surfaces
+        label index = surface_.find(nbrSurface);
+        if (index == -1)
+        {
+            // Append
+            origin_.append(neighbourInfo.origin_[i]);
+            distSqr_.append(magSqr(fc-neighbourInfo.origin_[i]));
+            surface_.append(nbrSurface);
+            //normal_.append(neighbourInfo.normal_[i]);
+            hasChanged = true;
+        }
+        else
+        {
+            hasChanged =
+                update(fc, index, neighbourInfo, i, tol, td)
+             || hasChanged;
+        }
+    }
+
+    return hasChanged;
+}
+
+
+template<class TrackingData>
+inline bool Foam::wallPoints::equal
+(
+    const wallPoints& rhs,
+    TrackingData& td
+) const
+{
+    return operator==(rhs);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline bool Foam::wallPoints::operator==
+(
+    const Foam::wallPoints& rhs
+) const
+{
+    return
+        surface_ == rhs.surface_
+     && distSqr_ == rhs.distSqr_
+     && origin_ == rhs.origin_;
+     //&& normal_ == rhs.normal_;
+}
+
+
+inline bool Foam::wallPoints::operator!=
+(
+    const Foam::wallPoints& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
+// ************************************************************************* //

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C

@@ -777,6 +777,19 @@ Foam::label Foam::snappyRefineDriver::surfaceProximityBlock
         if (debug)
         {
             const_cast<Time&>(mesh.time())++;
+            Pout<< "Writing gap blocking iteration "
+                << iter << " mesh to time " << meshRefiner_.timeName()
+                << endl;
+            meshRefiner_.write
+            (
+                meshRefinement::debugType(debug),
+                meshRefinement::writeType
+                (
+                    meshRefinement::writeLevel()
+                  | meshRefinement::WRITEMESH
+                ),
+                mesh.time().path()/meshRefiner_.timeName()
+            );
         }
     }
     return iter;
@@ -3011,6 +3024,18 @@ void Foam::snappyRefineDriver::doRefine
         100    // maxIter
     );
 
+    //// Re-remove cells inbetween two surfaces. The shell refinement/
+    //// directional shell refinement might have caused new small
+    //// gaps to be resolved. This is currently disabled since it finds
+    //// gaps just because it very occasionally walks around already removed
+    //// gaps and still finds 'opposite' surfaces. This probably need additional
+    //// path-length counting to avoid walking huge distances.
+    //surfaceProximityBlock
+    //(
+    //    refineParams,
+    //    1  //100     // maxIter
+    //);
+
     if
     (
         max(meshRefiner_.shells().nSmoothExpansion()) > 0