BUG: snappyHexMesh: use nearest meshed patch when exposing internal faces

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -377,14 +377,6 @@ private:
                 labelList& neiPatch
             ) const;
 
-            //- Determine patch for baffle using some heuristic (and not
-            //  surface)
-            label getBafflePatch
-            (
-                const labelList& facePatch,
-                const label faceI
-            ) const;
-
             //- Repatches external face or creates baffle for internal face
             //  with user specified patches (might be different for both sides).
             //  Returns label of added face.
@@ -450,6 +442,11 @@ private:
                 const labelList& globalToMasterPatch
             ) const;
 
+            //- Returns list with for every face the label of the nearest
+            //  patch. Any unreached face (disconnected mesh?) becomes
+            //  adaptPatchIDs[0]
+            labelList nearestPatch(const labelList& adaptPatchIDs) const;
+
             //- Returns list with for every internal face -1 or the patch
             //  they should be baffled into.
             labelList markFacesOnProblemCellsGeometric

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -131,96 +131,6 @@ Foam::label Foam::meshRefinement::createBaffle
 }
 
 
-// Get an estimate for the patch the internal face should get. Bit heuristic.
-Foam::label Foam::meshRefinement::getBafflePatch
-(
-    const labelList& facePatch,
-    const label faceI
-) const
-{
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-    // Loop over face points
-    // for each point check all faces patch IDs
-    // as soon as an ID >= 0 is found, break and assign that ID
-    // to the current face.
-    // Check first for real patch (so proper surface intersection and then
-    // in facePatch array for patches to block off faces
-
-    forAll(mesh_.faces()[faceI], fp)
-    {
-        label pointI = mesh_.faces()[faceI][fp];
-
-        forAll(mesh_.pointFaces()[pointI], pf)
-        {
-            label pFaceI = mesh_.pointFaces()[pointI][pf];
-
-            label patchI = patches.whichPatch(pFaceI);
-
-            if (patchI != -1 && !patches[patchI].coupled())
-            {
-                return patchI;
-            }
-            else if (facePatch[pFaceI] != -1)
-            {
-                return facePatch[pFaceI];
-            }
-        }
-    }
-
-    // Loop over owner and neighbour cells, looking for the first face with a
-    // valid patch number
-    const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
-
-    forAll(ownFaces, i)
-    {
-        label cFaceI = ownFaces[i];
-
-        label patchI = patches.whichPatch(cFaceI);
-
-        if (patchI != -1 && !patches[patchI].coupled())
-        {
-            return patchI;
-        }
-        else if (facePatch[cFaceI] != -1)
-        {
-            return facePatch[cFaceI];
-        }
-    }
-
-    if (mesh_.isInternalFace(faceI))
-    {
-        const cell& neiFaces = mesh_.cells()[mesh_.faceNeighbour()[faceI]];
-
-        forAll(neiFaces, i)
-        {
-            label cFaceI = neiFaces[i];
-
-            label patchI = patches.whichPatch(cFaceI);
-
-            if (patchI != -1 && !patches[patchI].coupled())
-            {
-                return patchI;
-            }
-            else if (facePatch[cFaceI] != -1)
-            {
-                return facePatch[cFaceI];
-            }
-        }
-    }
-
-    WarningIn
-    (
-        "meshRefinement::getBafflePatch(const labelList&, const label)"
-    )   << "Could not find boundary face neighbouring internal face "
-        << faceI << " with face centre " << mesh_.faceCentres()[faceI]
-        << nl
-        << "Using arbitrary patch " << 0 << " instead." << endl;
-
-    return 0;
-}
-
-
 //// Check if we are a boundary face and normal of surface does
 //// not align with test vector. In this case there'd probably be
 //// a freestanding 'baffle' so we might as well not create it.

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,6 +40,8 @@ License
 
 #include "snapParameters.H"
 #include "motionSmoother.H"
+#include "topoDistanceData.H"
+#include "FaceCellWave.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -372,6 +374,82 @@ bool Foam::meshRefinement::isCollapsedCell
 }
 
 
+Foam::labelList Foam::meshRefinement::nearestPatch
+(
+    const labelList& adaptPatchIDs
+) const
+{
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    // 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
+    labelList nearestAdaptPatch(mesh_.nFaces(), adaptPatchIDs[0]);
+
+    bool haveWarned = false;
+    forAll(faceData, faceI)
+    {
+        if (!faceData[faceI].valid(deltaCalc.data()))
+        {
+            if (!haveWarned)
+            {
+                WarningIn("meshRefinement::nearestPatch(..)")
+                    << "Did not visit some faces, e.g. face " << faceI
+                    << " at " << mesh_.faceCentres()[faceI] << endl
+                    << "Assigning  these cells to patch "
+                    << adaptPatchIDs[0]
+                    << endl;
+                haveWarned = true;
+            }
+        }
+        else
+        {
+            nearestAdaptPatch[faceI] = faceData[faceI].data();
+        }
+    }
+
+    return nearestAdaptPatch;
+}
+
+
 // Returns list with for every internal face -1 or the patch they should
 // be baffled into. Gets run after createBaffles so all the unzoned surface
 // intersections have already been turned into baffles. (Note: zoned surfaces
@@ -390,9 +468,6 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     const labelList& pointLevel = meshCutter_.pointLevel();
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    // Per internal face (boundary faces not used) the patch that the
-    // baffle should get (or -1)
-    labelList facePatch(mesh_.nFaces(), -1);
 
     // Mark all points and edges on baffle patches (so not on any inlets,
     // outlets etc.)
@@ -406,9 +481,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 
     forAll(adaptPatchIDs, i)
     {
-        label patchI = adaptPatchIDs[i];
+        const polyPatch& pp = patches[adaptPatchIDs[i]];
-
-        const polyPatch& pp = patches[patchI];
 
         label faceI = pp.start();
 
@@ -426,6 +499,15 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
         }
     }
 
+
+    // Per face the nearest adaptPatch
+    const labelList nearestAdaptPatch(nearestPatch(adaptPatchIDs));
+
+
+    // Per internal face (boundary faces not used) the patch that the
+    // baffle should get (or -1)
+    labelList facePatch(mesh_.nFaces(), -1);
+
     // Swap neighbouring cell centres and cell level
     labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
     pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
@@ -467,7 +549,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 
                 if (facePatch[faceI] == -1 && mesh_.isInternalFace(faceI))
                 {
-                    facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                    facePatch[faceI] = nearestAdaptPatch[faceI];
                     nBaffleFaces++;
 
                     // Mark face as a 'boundary'
@@ -708,7 +790,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
                          && mesh_.isInternalFace(faceI)
                         )
                         {
-                            facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                            facePatch[faceI] = nearestAdaptPatch[faceI];
                             nBaffleFaces++;
 
                             // Mark face as a 'boundary'
@@ -795,11 +877,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
                              && mesh_.isInternalFace(faceI)
                             )
                             {
-                                facePatch[faceI] = getBafflePatch
+                                facePatch[faceI] = nearestAdaptPatch[faceI];
-                                (
-                                    facePatch,
-                                    faceI
-                                );
                                 nBaffleFaces++;
 
                                 // Mark face as a 'boundary'
@@ -886,7 +964,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
                 }
                 else
                 {
-                    facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                    facePatch[faceI] = nearestAdaptPatch[faceI];
                     nBaffleFaces++;
 
                     // Do NOT update boundary data since this would grow blocked
@@ -943,7 +1021,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
                         }
                         else
                         {
-                            facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                            facePatch[faceI] = nearestAdaptPatch[faceI];
                             if (isMasterFace[faceI])
                             {
                                 nBaffleFaces++;
@@ -1091,6 +1169,9 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
     }
 
 
+    // Per face the nearest adaptPatch
+    const labelList nearestAdaptPatch(nearestPatch(meshedPatches()));
+
     // Per face (internal or coupled!) the patch that the
     // baffle should get (or -1).
     labelList facePatch(mesh_.nFaces(), -1);
@@ -1199,7 +1280,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
 
             if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
             {
-                facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
+                facePatch[iter.key()] = nearestAdaptPatch[iter.key()];
                 nBaffleFaces++;
 
                 //Pout<< "    " << iter.key()