problem cell removal

src/autoMesh/Make/files

@@ -15,6 +15,7 @@ $(autoHexMeshDriver)/pointData/pointData.C
 $(autoHexMesh)/meshRefinement/meshRefinementBaffles.C
 $(autoHexMesh)/meshRefinement/meshRefinement.C
 $(autoHexMesh)/meshRefinement/meshRefinementMerge.C
+$(autoHexMesh)/meshRefinement/meshRefinementProblemCells.C
 $(autoHexMesh)/meshRefinement/meshRefinementRefine.C
 $(autoHexMesh)/refinementSurfaces/refinementSurfaces.C
 $(autoHexMesh)/shellSurfaces/shellSurfaces.C

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C

@@ -518,7 +518,9 @@ void Foam::autoRefineDriver::baffleAndSplitMesh
     // be like boundary face from now on so not coupled anymore.
     meshRefiner_.baffleAndSplitMesh
     (
-        handleSnapProblems,
+        handleSnapProblems,             // detect&remove potential snap problem
+        false,                          // perpendicular edge connected cells
+        scalarField(0),                 // per region perpendicular angle
         !handleSnapProblems,            // merge free standing baffles?
         motionDict,
         const_cast<Time&>(mesh.time()),
@@ -592,10 +594,14 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
         const_cast<Time&>(mesh.time())++;
     }
 
+    const scalarField& perpAngle = meshRefiner_.surfaces().perpendicularAngle();
+
     meshRefiner_.baffleAndSplitMesh
     (
         handleSnapProblems,
-        false,                  // merge free standing baffles?
+        handleSnapProblems,                 // remove perp edge connected cells
+        perpAngle,                          // perp angle
+        false,                              // merge free standing baffles?
         motionDict,
         const_cast<Time&>(mesh.time()),
         globalToPatch_,

src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H

@@ -36,6 +36,7 @@ SourceFiles
     meshRefinement.C
     meshRefinementBaffles.C
     meshRefinementMerge.C
+    meshRefinementProblemCells.C
     meshRefinementRefine.C
 
 \*---------------------------------------------------------------------------*/
@@ -51,6 +52,7 @@ SourceFiles
 #include "indirectPrimitivePatch.H"
 #include "pointFieldsFwd.H"
 #include "Tuple2.H"
+#include "pointIndexHit.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -345,6 +347,8 @@ private:
                 polyTopoChange& meshMod
             ) const;
 
+        // Problem cell handling
+
             //- Helper function to mark face as being on 'boundary'. Used by
             //  markFacesOnProblemCells
             void markBoundaryFace
@@ -355,15 +359,32 @@ private:
                 boolList& isBoundaryPoint
             ) const;
 
+            void findNearest
+            (
+                const labelList& meshFaces,
+                List<pointIndexHit>& nearestInfo,
+                labelList& nearestSurface,
+                labelList& nearestRegion,
+                vectorField& nearestNormal
+            ) const;
+
+            Map<label> findEdgeConnectedProblemCells
+            (
+                const scalarField& perpendicularAngle,
+                const labelList&
+            ) const;
+
             //- Returns list with for every internal face -1 or the patch
-            //  they should be baffled into.
+            //  they should be baffled into. If removeEdgeConnectedCells is set
+            //  removes cells based on perpendicularAngle.
             labelList markFacesOnProblemCells
             (
+                const bool removeEdgeConnectedCells,
+                const scalarField& perpendicularAngle,
                 const labelList& globalToPatch
             ) const;
 
-            //- Returns list with for every internal face -1 or the patch
-            //  they should be baffled into.
+            //- Initial test of marking faces using geometric information.
             labelList markFacesOnProblemCellsGeometric
             (
                 const dictionary& motionDict,
@@ -589,6 +610,8 @@ public:
             void baffleAndSplitMesh
             (
                 const bool handleSnapProblems,
+                const bool removeEdgeConnectedCells,
+                const scalarField& perpendicularAngle,
                 const bool mergeFreeStanding,
                 const dictionary& motionDict,
                 Time& runTime,

src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -53,8 +53,10 @@ Foam::refinementSurfaces::refinementSurfaces
 {
     labelList globalMinLevel(surfaceDicts.size(), 0);
     labelList globalMaxLevel(surfaceDicts.size(), 0);
+    scalarField globalAngle(surfaceDicts.size(), -GREAT);
     List<Map<label> > regionMinLevel(surfaceDicts.size());
     List<Map<label> > regionMaxLevel(surfaceDicts.size());
+    List<Map<scalar> > regionAngle(surfaceDicts.size());
 
     //wordList globalPatchType(surfaceDicts.size());
     //List<HashTable<word> > regionPatchType(surfaceDicts.size());
@@ -80,6 +82,12 @@ Foam::refinementSurfaces::refinementSurfaces
             dict.lookup("zoneInside") >> zoneInside_[surfI];
         }
 
+        // Global perpendicular angle
+        if (dict.found("perpendicularAngle"))
+        {
+            globalAngle[surfI] = readScalar(dict.lookup("perpendicularAngle"));
+        }
+
         //// Global patch name per surface
         //if (dict.found("patchType"))
         //{
@@ -130,6 +138,15 @@ Foam::refinementSurfaces::refinementSurfaces
                         << exit(FatalError);
                 }
                 regionMaxLevel[surfI].insert(regionI, max);
+
+                if (regionDict.found("perpendicularAngle"))
+                {
+                    regionAngle[surfI].insert
+                    (
+                        regionI,
+                        readScalar(regionDict.lookup("perpendicularAngle"))
+                    );
+                }
             }
         }
     }
@@ -170,6 +187,8 @@ Foam::refinementSurfaces::refinementSurfaces
     minLevel_ = 0;
     maxLevel_.setSize(nRegions);
     maxLevel_ = 0;
+    perpendicularAngle_.setSize(nRegions);
+    perpendicularAngle_ = -GREAT;
     //patchName_.setSize(nRegions);
     //patchType_.setSize(nRegions);
 
@@ -182,6 +201,7 @@ Foam::refinementSurfaces::refinementSurfaces
         {
             minLevel_[regionOffset_[surfI] + i] = globalMinLevel[surfI];
             maxLevel_[regionOffset_[surfI] + i] = globalMaxLevel[surfI];
+            perpendicularAngle_[regionOffset_[surfI] + i] = globalAngle[surfI];
         }
 
         // Overwrite with region specific information
@@ -191,6 +211,7 @@ Foam::refinementSurfaces::refinementSurfaces
 
             minLevel_[globalRegionI] = iter();
             maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
+            perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
 
             // Check validity
             if
@@ -240,8 +261,10 @@ Foam::refinementSurfaces::refinementSurfaces
 {
     labelList globalMinLevel(surfacesDict.size(), 0);
     labelList globalMaxLevel(surfacesDict.size(), 0);
+    scalarField globalAngle(surfacesDict.size(), -GREAT);
     List<Map<label> > regionMinLevel(surfacesDict.size());
     List<Map<label> > regionMaxLevel(surfacesDict.size());
+    List<Map<scalar> > regionAngle(surfacesDict.size());
 
     label surfI = 0;
     forAllConstIter(dictionary, surfacesDict, iter)
@@ -274,6 +297,12 @@ Foam::refinementSurfaces::refinementSurfaces
             dict.lookup("zoneInside") >> zoneInside_[surfI];
         }
 
+        // Global perpendicular angle
+        if (dict.found("perpendicularAngle"))
+        {
+            globalAngle[surfI] = readScalar(dict.lookup("perpendicularAngle"));
+        }
+
         if (dict.found("regions"))
         {
             const dictionary& regionsDict = dict.subDict("regions");
@@ -306,6 +335,15 @@ Foam::refinementSurfaces::refinementSurfaces
 
                     regionMinLevel[surfI].insert(regionI, refLevel[0]);
                     regionMaxLevel[surfI].insert(regionI, refLevel[1]);
+
+                    if (regionDict.found("perpendicularAngle"))
+                    {
+                        regionAngle[surfI].insert
+                        (
+                            regionI,
+                            readScalar(regionDict.lookup("perpendicularAngle"))
+                        );
+                    }
                 }
             }
         }
@@ -326,6 +364,8 @@ Foam::refinementSurfaces::refinementSurfaces
     minLevel_ = 0;
     maxLevel_.setSize(nRegions);
     maxLevel_ = 0;
+    perpendicularAngle_.setSize(nRegions);
+    perpendicularAngle_ = -GREAT;
 
 
     forAll(globalMinLevel, surfI)
@@ -337,6 +377,7 @@ Foam::refinementSurfaces::refinementSurfaces
         {
             minLevel_[regionOffset_[surfI] + i] = globalMinLevel[surfI];
             maxLevel_[regionOffset_[surfI] + i] = globalMaxLevel[surfI];
+            perpendicularAngle_[regionOffset_[surfI] + i] = globalAngle[surfI];
         }
 
         // Overwrite with region specific information
@@ -346,6 +387,7 @@ Foam::refinementSurfaces::refinementSurfaces
 
             minLevel_[globalRegionI] = iter();
             maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
+            perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
 
             // Check validity
             if
@@ -454,8 +496,6 @@ void Foam::refinementSurfaces::setMinLevelFields
     const shellSurfaces& shells
 )
 {
-    //minLevelFields_.setSize(surfaces_.size());
-
     forAll(surfaces_, surfI)
     {
         const searchableSurface& geom = allGeometry_[surfaces_[surfI]];

src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H

@@ -90,6 +90,9 @@ class refinementSurfaces
         //- From global region number to refinement level
         labelList maxLevel_;
 
+        //- From global region number to perpendicular angle
+        scalarField perpendicularAngle_;
+
 
     // Private Member Functions
 
@@ -178,6 +181,12 @@ public:
                 return maxLevel_;
             }
 
+            //- From global region number to perpendicular angle
+            const scalarField& perpendicularAngle() const
+            {
+                return perpendicularAngle_;
+            }
+
 
         // Helper