Snappy hex mesh proximity check

src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C

@@ -58,6 +58,8 @@ License
 #include "fvMeshTools.H"
 #include "motionSmoother.H"
 #include "faceSet.H"
+#include "topoDistanceData.H"
+#include "FaceCellWave.H"
 
 // Leak path
 #include "shortestPathSet.H"
@@ -358,6 +360,218 @@ void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
 }
 
 
+Foam::labelList Foam::meshRefinement::nearestPatch
+(
+    const labelList& adaptPatchIDs
+) const
+{
+    // Determine nearest patch for all mesh faces. Used when removing cells
+    // to give some reasonable patch to exposed faces.
+
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    labelList nearestAdaptPatch;
+
+    if (adaptPatchIDs.size())
+    {
+        nearestAdaptPatch.setSize(mesh_.nFaces(), adaptPatchIDs[0]);
+
+
+        // Count number of faces in adaptPatchIDs
+        label nFaces = 0;
+        forAll(adaptPatchIDs, i)
+        {
+            const polyPatch& pp = patches[adaptPatchIDs[i]];
+            nFaces += pp.size();
+        }
+
+        // Field on cells and faces.
+        List<topoDistanceData> cellData(mesh_.nCells());
+        List<topoDistanceData> faceData(mesh_.nFaces());
+
+        // Start of changes
+        labelList patchFaces(nFaces);
+        List<topoDistanceData> patchData(nFaces);
+        nFaces = 0;
+        forAll(adaptPatchIDs, i)
+        {
+            label patchi = adaptPatchIDs[i];
+            const polyPatch& pp = patches[patchi];
+
+            forAll(pp, i)
+            {
+                patchFaces[nFaces] = pp.start()+i;
+                patchData[nFaces] = topoDistanceData(patchi, 0);
+                nFaces++;
+            }
+        }
+
+        // Propagate information inwards
+        FaceCellWave<topoDistanceData> deltaCalc
+        (
+            mesh_,
+            patchFaces,
+            patchData,
+            faceData,
+            cellData,
+            mesh_.globalData().nTotalCells()+1
+        );
+
+        // And extract
+
+        bool haveWarned = false;
+        forAll(faceData, facei)
+        {
+            if (!faceData[facei].valid(deltaCalc.data()))
+            {
+                if (!haveWarned)
+                {
+                    WarningInFunction
+                        << "Did not visit some faces, e.g. face " << facei
+                        << " at " << mesh_.faceCentres()[facei] << endl
+                        << "Assigning these faces to patch "
+                        << adaptPatchIDs[0]
+                        << endl;
+                    haveWarned = true;
+                }
+            }
+            else
+            {
+                nearestAdaptPatch[facei] = faceData[facei].data();
+            }
+        }
+    }
+    else
+    {
+        // Use patch 0
+        nearestAdaptPatch.setSize(mesh_.nFaces(), 0);
+    }
+
+    return nearestAdaptPatch;
+}
+
+
+Foam::labelList Foam::meshRefinement::nearestIntersection
+(
+    const labelList& surfacesToTest,
+    const label defaultRegion
+) const
+{
+    // Determine nearest intersection for all mesh faces. Used when removing
+    // cells to give some reasonable patch to exposed faces. Use this
+    // function instead of nearestPatch if you don't have patches yet.
+
+
+    // Swap neighbouring cell centres and cell level
+    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
+    pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
+    calcNeighbourData(neiLevel, neiCc);
+
+
+    // Collect segments
+    // ~~~~~~~~~~~~~~~~
+
+    const labelList testFaces(intersectedFaces());
+
+    pointField start(testFaces.size());
+    pointField end(testFaces.size());
+    labelList minLevel(testFaces.size());
+
+    calcCellCellRays
+    (
+        neiCc,
+        neiLevel,
+        testFaces,
+        start,
+        end,
+        minLevel
+    );
+
+    // Do tests in one go
+    labelList surface1;
+    List<pointIndexHit> hit1;
+    labelList region1;
+    labelList surface2;
+    List<pointIndexHit> hit2;
+    labelList region2;
+    surfaces_.findNearestIntersection
+    (
+        surfacesToTest,
+        start,
+        end,
+
+        surface1,
+        hit1,
+        region1,
+        surface2,
+        hit2,
+        region2
+    );
+
+    labelList nearestRegion(mesh_.nFaces(), defaultRegion);
+
+    // Field on cells and faces.
+    List<topoDistanceData> cellData(mesh_.nCells());
+    List<topoDistanceData> faceData(mesh_.nFaces());
+
+    // Start walking from all intersected faces
+    DynamicList<label> patchFaces(start.size());
+    DynamicList<topoDistanceData> patchData(start.size());
+    forAll(start, i)
+    {
+        label facei = testFaces[i];
+        if (surface1[i] != -1)
+        {
+            patchFaces.append(facei);
+            label regioni = surfaces_.globalRegion(surface1[i], region1[i]);
+            patchData.append(topoDistanceData(regioni, 0));
+        }
+        else if (surface2[i] != -1)
+        {
+            patchFaces.append(facei);
+            label regioni = surfaces_.globalRegion(surface2[i], region2[i]);
+            patchData.append(topoDistanceData(regioni, 0));
+        }
+    }
+
+    // Propagate information inwards
+    FaceCellWave<topoDistanceData> deltaCalc
+    (
+        mesh_,
+        patchFaces,
+        patchData,
+        faceData,
+        cellData,
+        mesh_.globalData().nTotalCells()+1
+    );
+
+    // And extract
+
+    bool haveWarned = false;
+    forAll(faceData, facei)
+    {
+        if (!faceData[facei].valid(deltaCalc.data()))
+        {
+            if (!haveWarned)
+            {
+                WarningInFunction
+                    << "Did not visit some faces, e.g. face " << facei
+                    << " at " << mesh_.faceCentres()[facei] << endl
+                    << "Assigning these faces to global region "
+                    << defaultRegion << endl;
+                haveWarned = true;
+            }
+        }
+        else
+        {
+            nearestRegion[facei] = faceData[facei].data();
+        }
+    }
+
+    return nearestRegion;
+}
+
+
 void Foam::meshRefinement::testSyncPointList
 (
     const string& msg,
@@ -3147,7 +3361,7 @@ void Foam::meshRefinement::write
         dumpRefinementLevel();
     }
 
-    if (debugFlags & OBJINTERSECTIONS && prefix.size())
+    if ((debugFlags & OBJINTERSECTIONS) && prefix.size())
     {
         dumpIntersections(prefix);
     }