ENH: shanppyHexMesh: disable gap detection on same surface. See #1463.

etc/caseDicts/annotated/snappyHexMeshDict

@@ -282,6 +282,7 @@ castellatedMeshControls
         //    //           cell selection. Default is 'mixed' i.e. keep cells
         //    //           whilst doing the gap-level refinement.
         //    //gapMode inside;  // inside/outside/mixed
+        //    //gapSelf false;   // ignore gaps in same surface
         //}
 
         //wakeBox

src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H

@@ -434,6 +434,7 @@ private:
                 const label nRefine,
 
                 labelList& cellMap,
+                labelList& gapShell,
                 List<FixedList<label, 3>>& shellGapInfo,
                 List<volumeType>& shellGapMode
             ) const;

src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2015-2016 OpenCFD Ltd.
+    Copyright (C) 2015-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -94,6 +94,7 @@ Foam::label Foam::meshRefinement::markSurfaceGapRefinement
         // Collect cells to test for inside/outside in shell
         labelList cellToCompact(mesh_.nCells(), -1);
         labelList bFaceToCompact(mesh_.nBoundaryFaces(), -1);
+        labelList gapShell;
         List<FixedList<label, 3>> shellGapInfo;
         List<volumeType> shellGapMode;
         {
@@ -135,22 +136,69 @@ Foam::label Foam::meshRefinement::markSurfaceGapRefinement
             (
                 compactToCc,
                 compactToLevel,
+
+                gapShell,
                 shellGapInfo,
                 shellGapMode
             );
         }
 
 
+        //const fileName dir(mesh_.time().path()/timeName());
+        //if (debug)
+        //{
+        //    mkDir(dir);
+        //    OBJstream insideStr(dir/"insideShell.obj");
+        //    OBJstream outsideStr(dir/"outsideShell.obj");
+        //    Pout<< "Writing points to:" << nl
+        //        << "    inside : " << insideStr.name() << nl
+        //        << "    outside: " << outsideStr.name() << nl
+        //        << endl;
+        //
+        //    forAll(cellToCompact, celli)
+        //    {
+        //        const label compacti = cellToCompact[celli];
+        //
+        //        if (compacti != -1)
+        //        {
+        //            if (gapShell[compacti] != -1)
+        //            {
+        //                insideStr.write(mesh_.cellCentres()[celli]);
+        //            }
+        //            else
+        //            {
+        //                outsideStr.write(mesh_.cellCentres()[celli]);
+        //            }
+        //        }
+        //    }
+        //    forAll(bFaceToCompact, bFacei)
+        //    {
+        //        const label compacti = bFaceToCompact[bFacei];
+        //        if (compacti != -1)
+        //        {
+        //            if (gapShell[compacti] != -1)
+        //            {
+        //                insideStr.write(neiCc[bFacei]);
+        //            }
+        //            else
+        //            {
+        //                outsideStr.write(neiCc[bFacei]);
+        //            }
+        //        }
+        //    }
+        //}
+
+
         const List<FixedList<label, 3>>& extendedGapLevel =
             surfaces_.extendedGapLevel();
         const List<volumeType>& extendedGapMode =
             surfaces_.extendedGapMode();
+        const boolList& extendedGapSelf = surfaces_.gapSelf();
 
         labelList ccSurface1;
         List<pointIndexHit> ccHit1;
         labelList ccRegion1;
         vectorField ccNormal1;
-
         {
             labelList ccSurface2;
             List<pointIndexHit> ccHit2;
@@ -218,6 +266,7 @@ Foam::label Foam::meshRefinement::markSurfaceGapRefinement
                         shellGapMode[compactI],
                         extendedGapLevel[globalRegionI],
                         extendedGapMode[globalRegionI],
+
                         gapInfo,
                         gapMode
                     );
@@ -266,6 +315,7 @@ Foam::label Foam::meshRefinement::markSurfaceGapRefinement
                             shellGapMode[compactI],
                             extendedGapLevel[globalRegionI],
                             extendedGapMode[globalRegionI],
+
                             gapInfo,
                             gapMode
                         );
@@ -320,6 +370,7 @@ Foam::label Foam::meshRefinement::markSurfaceGapRefinement
                             shellGapMode[compactI],
                             extendedGapLevel[globalRegionI],
                             extendedGapMode[globalRegionI],
+
                             gapInfo,
                             gapMode
                         );
@@ -402,7 +453,33 @@ Foam::label Foam::meshRefinement::markSurfaceGapRefinement
 
         forAll(surf1, i)
         {
-            if (surf1[i] != -1 && surf2[i] != -1)
+            // Combine selfProx of shell and surfaces. Ignore regions for
+            // now
+            const label cellI = cellMap[i];
+            const label shelli =
+            (
+                cellToCompact[cellI] != -1
+              ? gapShell[cellToCompact[cellI]]
+              : -1
+            );
+
+            bool selfProx = true;
+            if (shelli != -1)
+            {
+                selfProx = shells_.gapSelf()[shelli][0];
+            }
+            if (surf1[i] != -1 && selfProx)
+            {
+                const label globalRegioni = surfaces_.globalRegion(surf1[i], 0);
+                selfProx = extendedGapSelf[globalRegioni];
+            }
+
+            if
+            (
+                surf1[i] != -1
+             && surf2[i] != -1
+             && (surf2[i] != surf1[i] || selfProx)
+            )
             {
                 // Found intersection with surface. Check opposite normal.
                 label cellI = cellMap[i];
@@ -1000,6 +1077,7 @@ void Foam::meshRefinement::selectGapCandidates
     const label nRefine,
 
     labelList& cellMap,
+    labelList& gapShell,
     List<FixedList<label, 3>>& shellGapInfo,
     List<volumeType>& shellGapMode
 ) const
@@ -1029,6 +1107,8 @@ void Foam::meshRefinement::selectGapCandidates
     (
         pointField(cellCentres, cellMap),
         labelUIndList(cellLevel, cellMap)(),
+
+        gapShell,
         shellGapInfo,
         shellGapMode
     );
@@ -1051,6 +1131,7 @@ void Foam::meshRefinement::selectGapCandidates
 
     map.setSize(compactI);
     cellMap = labelUIndList(cellMap, map)();
+    gapShell = labelUIndList(gapShell, map)();
     shellGapInfo = UIndirectList<FixedList<label, 3>>(shellGapInfo, map)();
     shellGapMode = UIndirectList<volumeType>(shellGapMode, map)();
 }
@@ -1062,6 +1143,7 @@ void Foam::meshRefinement::mergeGapInfo
     const volumeType shellGapMode,
     const FixedList<label, 3>& surfGapInfo,
     const volumeType surfGapMode,
+
     FixedList<label, 3>& gapInfo,
     volumeType& gapMode
 ) const
@@ -1114,6 +1196,7 @@ Foam::label Foam::meshRefinement::markInternalGapRefinement
     const List<FixedList<label, 3>>& extendedGapLevel =
         surfaces_.extendedGapLevel();
     const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
+    const boolList& extendedGapSelf = surfaces_.gapSelf();
 
     // Get the gap level for the shells
     const labelList maxLevel(shells_.maxGapLevel());
@@ -1124,6 +1207,7 @@ Foam::label Foam::meshRefinement::markInternalGapRefinement
     {
         // Collect cells to test
         labelList cellMap;
+        labelList gapShell;
         List<FixedList<label, 3>> shellGapInfo;
         List<volumeType> shellGapMode;
         selectGapCandidates
@@ -1132,6 +1216,7 @@ Foam::label Foam::meshRefinement::markInternalGapRefinement
             nRefine,
 
             cellMap,
+            gapShell,
             shellGapInfo,
             shellGapMode
         );
@@ -1194,8 +1279,10 @@ Foam::label Foam::meshRefinement::markInternalGapRefinement
                 (
                     shellGapInfo[i],
                     shellGapMode[i],
+
                     extendedGapLevel[globalRegionI],
                     extendedGapMode[globalRegionI],
+
                     gapInfo,
                     gapMode
                 );
@@ -1288,15 +1375,35 @@ Foam::label Foam::meshRefinement::markInternalGapRefinement
         // Extract cell based gap size
         forAll(surf1, i)
         {
-            if (surf1[i] != -1 && surf2[i] != -1)
+            // Combine selfProx of shell and surfaces. Ignore regions for
+            // now
+            const label shelli = gapShell[map[i]];
+
+            bool selfProx = true;
+            if (shelli != -1)
+            {
+                selfProx = shells_.gapSelf()[shelli][0];
+            }
+            if (surf1[i] != -1 && selfProx)
+            {
+                const label globalRegioni = surfaces_.globalRegion(surf1[i], 0);
+                selfProx = extendedGapSelf[globalRegioni];
+            }
+
+            if
+            (
+                surf1[i] != -1
+             && surf2[i] != -1
+             && (surf2[i] != surf1[i] || selfProx)
+            )
             {
                 // Found intersections with surface. Check for
                 // - small gap
                 // - coplanar normals
 
-                label cellI = cellMap[i];
+                const label cellI = cellMap[i];
 
-                scalar d2 = magSqr(hit1[i].hitPoint()-hit2[i].hitPoint());
+                const scalar d2 = magSqr(hit1[i].hitPoint()-hit2[i].hitPoint());
 
                 if
                 (
@@ -1479,6 +1586,7 @@ Foam::label Foam::meshRefinement::markSmallFeatureRefinement
     const List<FixedList<label, 3>>& extendedGapLevel =
         surfaces_.extendedGapLevel();
     const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
+    const boolList& extendedGapSelf = surfaces_.gapSelf();
 
     label oldNRefine = nRefine;
 
@@ -1535,10 +1643,13 @@ Foam::label Foam::meshRefinement::markSmallFeatureRefinement
         // applicable specification (minLevel <= celllevel < maxLevel)
         List<FixedList<label, 3>> shellGapInfo;
         List<volumeType> shellGapMode;
+        labelList gapShell;
         shells_.findHigherGapLevel
         (
             ctrs,
             labelList(ctrs.size(), Zero),
+
+            gapShell,
             shellGapInfo,
             shellGapMode
         );
@@ -1566,8 +1677,10 @@ Foam::label Foam::meshRefinement::markSmallFeatureRefinement
                 (
                     shellGapInfo[i],
                     shellGapMode[i],
+
                     extendedGapLevel[globalRegionI],
                     extendedGapMode[globalRegionI],
+
                     gapInfo,
                     gapMode
                 );
@@ -1647,7 +1760,25 @@ Foam::label Foam::meshRefinement::markSmallFeatureRefinement
 
         forAll(surfaceHit, i)
         {
-            if (surfaceHit[i] != -1) // && surf2[i] != -1)
+            // Combine selfProx of shell and surfaces. Ignore regions for
+            // now
+            const label shelli = gapShell[map[i]];
+            bool selfProx = true;
+            if (shelli != -1)
+            {
+                selfProx = shells_.gapSelf()[shelli][0];
+            }
+            if (surfI != -1 && selfProx)
+            {
+                const label globalRegioni = surfaces_.globalRegion(surfI, 0);
+                selfProx = extendedGapSelf[globalRegioni];
+            }
+
+            if
+            (
+                surfaceHit[i] != -1
+             && (surfaceHit[i] != surfI || selfProx)
+            )
             {
                 // Found intersection with surface. Check coplanar normals.
                 label cellI = cellMap[i];

src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -201,6 +201,7 @@ Foam::refinementSurfaces::refinementSurfaces
 
     List<FixedList<label, 3>> globalGapLevel(surfI);
     List<volumeType> globalGapMode(surfI);
+    boolList globalGapSelf(surfI);
 
     scalarField globalAngle(surfI, -GREAT);
     PtrList<dictionary> globalPatchInfo(surfI);
@@ -213,6 +214,7 @@ Foam::refinementSurfaces::refinementSurfaces
     List<Map<label>> regionLevelIncr(surfI);
     List<Map<FixedList<label, 3>>> regionGapLevel(surfI);
     List<Map<volumeType>> regionGapMode(surfI);
+    List<Map<bool>> regionGapSelf(surfI);
     List<Map<scalar>> regionAngle(surfI);
     List<Map<autoPtr<dictionary>>> regionPatchInfo(surfI);
     List<Map<label>> regionBlockLevel(surfI);
@@ -300,6 +302,8 @@ Foam::refinementSurfaces::refinementSurfaces
                     << exit(FatalIOError);
             }
 
+            globalGapSelf[surfI] =
+                dict.getOrDefault<bool>("gapSelf", true);
 
             const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
 
@@ -418,7 +422,15 @@ Foam::refinementSurfaces::refinementSurfaces
                                 << " gapMode:" << gapModeSpec.str()
                                 << exit(FatalIOError);
                         }
-
+                        regionGapSelf[surfI].insert
+                        (
+                            regionI,
+                            regionDict.getOrDefault<bool>
+                            (
+                                "gapSelf",
+                                true
+                            )
+                        );
 
                         if (regionDict.found("perpendicularAngle"))
                         {
@@ -484,6 +496,8 @@ Foam::refinementSurfaces::refinementSurfaces
     extendedGapLevel_ = nullGapLevel;
     extendedGapMode_.setSize(nRegions);
     extendedGapMode_ = volumeType::UNKNOWN;
+    selfProximity_.setSize(nRegions);
+    selfProximity_ = true;
     perpendicularAngle_.setSize(nRegions);
     perpendicularAngle_ = -GREAT;
     patchInfo_.setSize(nRegions);
@@ -507,6 +521,7 @@ Foam::refinementSurfaces::refinementSurfaces
               + globalLevelIncr[surfI];
             extendedGapLevel_[globalRegionI] = globalGapLevel[surfI];
             extendedGapMode_[globalRegionI] = globalGapMode[surfI];
+            selfProximity_[globalRegionI] = globalGapSelf[surfI];
             perpendicularAngle_[globalRegionI] = globalAngle[surfI];
             if (globalPatchInfo.set(surfI))
             {
@@ -533,6 +548,8 @@ Foam::refinementSurfaces::refinementSurfaces
                 regionGapLevel[surfI][iter.key()];
             extendedGapMode_[globalRegionI] =
                 regionGapMode[surfI][iter.key()];
+            selfProximity_[globalRegionI] =
+                regionGapSelf[surfI][iter.key()];
         }
         forAllConstIters(regionAngle[surfI], iter)
         {

src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H

@@ -102,6 +102,10 @@ class refinementSurfaces
         //- From global region number to side of surface to detect
         List<volumeType> extendedGapMode_;
 
+        //- From global region number to whether to allow selfProximity
+        //  (in gap refinement)
+        boolList selfProximity_;
+
         //- From global region number to perpendicular angle
         scalarField perpendicularAngle_;
 
@@ -241,6 +245,13 @@ public:
                 return extendedGapMode_;
             }
 
+            //- From global region number to whether to detect gaps to same
+            //  surface (in gap refinement)
+            const boolList& gapSelf() const
+            {
+                return selfProximity_;
+            }
+
             //- From global region number to perpendicular angle
             const scalarField& perpendicularAngle() const
             {

src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C

@@ -610,6 +610,7 @@ Foam::shellSurfaces::shellSurfaces
 
     extendedGapLevel_.setSize(shellI);
     extendedGapMode_.setSize(shellI);
+    selfProximity_.setSize(shellI);
 
     FixedList<label, 3> nullGapLevel;
     nullGapLevel[0] = 0;
@@ -724,6 +725,13 @@ Foam::shellSurfaces::shellSurfaces
             extendedGapMode_[shellI] =
                 volumeType("gapMode", dict, volumeType::MIXED);
 
+            // Detect self-intersections
+            selfProximity_[shellI].setSize
+            (
+                regionNames.size(),
+                dict.getOrDefault<bool>("gapSelf", true)
+            );
+
 
             // Override on a per-region basis?
 
@@ -756,6 +764,13 @@ Foam::shellSurfaces::shellSurfaces
                                 regionDict,
                                 volumeType::MIXED
                             );
+
+                        selfProximity_[shellI][regionI] =
+                            regionDict.getOrDefault<bool>
+                            (
+                                "gapSelf",
+                                true
+                            );
                     }
                 }
             }

src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2015 OpenCFD Ltd.
+    Copyright (C) 2015-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -116,6 +116,9 @@ private:
             //- Per shell, per region the small-gap level specification
             List<List<volumeType>> extendedGapMode_;
 
+            //- Per shell, per region whether to allow selfProximity refinement
+            boolListList selfProximity_;
+
 
     // Private data
 
@@ -206,6 +209,13 @@ public:
                 return extendedGapMode_;
             }
 
+            //- Per shell, per region whether to test for gap with same surface
+            const boolListList& gapSelf() const
+            {
+                return selfProximity_;
+            }
+
+
         // Query
 
             //- Highest shell level