ENH: autoHexMesh: added geometric problem cell removal

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C

@@ -35,6 +35,7 @@ License
 #include "shellSurfaces.H"
 #include "mapDistributePolyMesh.H"
 #include "unitConversion.H"
+#include "snapParameters.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -610,9 +611,16 @@ Foam::label Foam::autoRefineDriver::danglingCellRefine
             << " cells (out of " << mesh.globalData().nTotalCells()
             << ')' << endl;
 
-        // Stop when no cells to refine. No checking of minRefineCells since
+        // Stop when no cells to refine. After a few iterations check if too
-        // too few cells
+        // few cells
-        if (nCellsToRefine == 0)
+        if
+        (
+            nCellsToRefine == 0
+         || (
+                iter >= 1
+             && nCellsToRefine <= refineParams.minRefineCells()
+            )
+        )
         {
             Info<< "Stopping refining since too few cells selected."
                 << nl << endl;
@@ -870,6 +878,7 @@ Foam::label Foam::autoRefineDriver::shellRefine
 void Foam::autoRefineDriver::baffleAndSplitMesh
 (
     const refinementParameters& refineParams,
+    const snapParameters& snapParams,
     const bool handleSnapProblems,
     const dictionary& motionDict
 )
@@ -887,10 +896,17 @@ void Foam::autoRefineDriver::baffleAndSplitMesh
     meshRefiner_.baffleAndSplitMesh
     (
         handleSnapProblems,             // detect&remove potential snap problem
+
+        // Snap problem cell detection
+        snapParams,
+        refineParams.useTopologicalSnapDetection(),
         false,                          // perpendicular edge connected cells
         scalarField(0),                 // per region perpendicular angle
+
+        // Free standing baffles
         !handleSnapProblems,            // merge free standing baffles?
         refineParams.planarAngle(),
+
         motionDict,
         const_cast<Time&>(mesh.time()),
         globalToMasterPatch_,
@@ -951,6 +967,7 @@ void Foam::autoRefineDriver::zonify
 void Foam::autoRefineDriver::splitAndMergeBaffles
 (
     const refinementParameters& refineParams,
+    const snapParameters& snapParams,
     const bool handleSnapProblems,
     const dictionary& motionDict
 )
@@ -973,10 +990,17 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
     meshRefiner_.baffleAndSplitMesh
     (
         handleSnapProblems,
+
+        // Snap problem cell detection
+        snapParams,
+        refineParams.useTopologicalSnapDetection(),
         handleSnapProblems,                 // remove perp edge connected cells
         perpAngle,                          // perp angle
+
+        // Free standing baffles
         true,                               // merge free standing baffles?
         refineParams.planarAngle(),         // planar angle
+
         motionDict,
         const_cast<Time&>(mesh.time()),
         globalToMasterPatch_,
@@ -1081,6 +1105,7 @@ void Foam::autoRefineDriver::doRefine
 (
     const dictionary& refineDict,
     const refinementParameters& refineParams,
+    const snapParameters& snapParams,
     const bool prepareForSnapping,
     const dictionary& motionDict
 )
@@ -1146,13 +1171,25 @@ void Foam::autoRefineDriver::doRefine
 
     // Introduce baffles at surface intersections. Remove sections unreachable
     // from keepPoint.
-    baffleAndSplitMesh(refineParams, prepareForSnapping, motionDict);
+    baffleAndSplitMesh
+    (
+        refineParams,
+        snapParams,
+        prepareForSnapping,
+        motionDict
+    );
 
     // Mesh is at its finest. Do optional zoning.
     zonify(refineParams);
 
     // Pull baffles apart
-    splitAndMergeBaffles(refineParams, prepareForSnapping, motionDict);
+    splitAndMergeBaffles
+    (
+        refineParams,
+        snapParams,
+        prepareForSnapping,
+        motionDict
+    );
 
     // Do something about cells with refined faces on the boundary
     if (prepareForSnapping)

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H

@@ -43,6 +43,8 @@ namespace Foam
 
 // Forward declaration of classes
 class refinementParameters;
+class snapParameters;
+
 class meshRefinement;
 class decompositionMethod;
 class fvMeshDistribute;
@@ -121,6 +123,7 @@ class autoRefineDriver
         void baffleAndSplitMesh
         (
             const refinementParameters& refineParams,
+            const snapParameters& snapParams,
             const bool handleSnapProblems,
             const dictionary& motionDict
         );
@@ -131,6 +134,7 @@ class autoRefineDriver
         void splitAndMergeBaffles
         (
             const refinementParameters& refineParams,
+            const snapParameters& snapParams,
             const bool handleSnapProblems,
             const dictionary& motionDict
         );
@@ -176,6 +180,7 @@ public:
         (
             const dictionary& refineDict,
             const refinementParameters& refineParams,
+            const snapParameters& snapParams,
             const bool prepareForSnapping,
             const dictionary& motionDict
         );

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C

@@ -121,7 +121,7 @@ Foam::pointField Foam::autoSnapDriver::smoothPatchDisplacement
 (
     const motionSmoother& meshMover,
     const List<labelPair>& baffles
-) const
+)
 {
     const indirectPrimitivePatch& pp = meshMover.patch();
 
@@ -615,14 +615,14 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::mergeZoneBaffles
 
 Foam::scalarField Foam::autoSnapDriver::calcSnapDistance
 (
+    const fvMesh& mesh,
     const snapParameters& snapParams,
     const indirectPrimitivePatch& pp
-) const
+)
 {
     const edgeList& edges = pp.edges();
     const labelListList& pointEdges = pp.pointEdges();
     const pointField& localPoints = pp.localPoints();
-    const fvMesh& mesh = meshRefiner_.mesh();
 
     scalarField maxEdgeLen(localPoints.size(), -GREAT);
 
@@ -655,13 +655,14 @@ Foam::scalarField Foam::autoSnapDriver::calcSnapDistance
 
 void Foam::autoSnapDriver::preSmoothPatch
 (
+    const meshRefinement& meshRefiner,
     const snapParameters& snapParams,
     const label nInitErrors,
     const List<labelPair>& baffles,
     motionSmoother& meshMover
-) const
+)
 {
-    const fvMesh& mesh = meshRefiner_.mesh();
+    const fvMesh& mesh = meshRefiner.mesh();
 
     labelList checkFaces;
 
@@ -724,11 +725,11 @@ void Foam::autoSnapDriver::preSmoothPatch
     {
         const_cast<Time&>(mesh.time())++;
         Info<< "Writing patch smoothed mesh to time "
-            << meshRefiner_.timeName() << '.' << endl;
+            << meshRefiner.timeName() << '.' << endl;
-        meshRefiner_.write
+        meshRefiner.write
         (
             debug,
-            mesh.time().path()/meshRefiner_.timeName()
+            mesh.time().path()/meshRefiner.timeName()
         );
         Info<< "Dumped mesh in = "
             << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
@@ -742,12 +743,11 @@ void Foam::autoSnapDriver::preSmoothPatch
 // Get (pp-local) indices of points that are both on zone and on patched surface
 Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
 (
+    const fvMesh& mesh,
     const indirectPrimitivePatch& pp,
     const word& zoneName
-) const
+)
 {
-    const fvMesh& mesh = meshRefiner_.mesh();
-
     label zoneI = mesh.faceZones().findZoneID(zoneName);
 
     if (zoneI == -1)
@@ -755,7 +755,7 @@ Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
         FatalErrorIn
         (
             "autoSnapDriver::getZoneSurfacePoints"
-            "(const indirectPrimitivePatch&, const word&)"
+            "(const fvMesh&, const indirectPrimitivePatch&, const word&)"
         )   << "Cannot find zone " << zoneName
             << exit(FatalError);
     }
@@ -793,17 +793,17 @@ Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
 
 Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
 (
+    const meshRefinement& meshRefiner,
     const scalarField& snapDist,
-    motionSmoother& meshMover
+    const indirectPrimitivePatch& pp
-) const
+)
 {
     Info<< "Calculating patchDisplacement as distance to nearest surface"
         << " point ..." << endl;
 
-    const indirectPrimitivePatch& pp = meshMover.patch();
     const pointField& localPoints = pp.localPoints();
-    const refinementSurfaces& surfaces = meshRefiner_.surfaces();
+    const refinementSurfaces& surfaces = meshRefiner.surfaces();
-    const fvMesh& mesh = meshRefiner_.mesh();
+    const fvMesh& mesh = meshRefiner.mesh();
 
     // Displacement per patch point
     vectorField patchDisp(localPoints.size(), vector::zero);
@@ -815,9 +815,9 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
 
         // Divide surfaces into zoned and unzoned
         labelList zonedSurfaces =
-            meshRefiner_.surfaces().getNamedSurfaces();
+            meshRefiner.surfaces().getNamedSurfaces();
         labelList unzonedSurfaces =
-            meshRefiner_.surfaces().getUnnamedSurfaces();
+            meshRefiner.surfaces().getUnnamedSurfaces();
 
 
         // 1. All points to non-interface surfaces
@@ -870,6 +870,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
             (
                 getZoneSurfacePoints
                 (
+                    mesh,
                     pp,
                     faceZoneNames[zoneSurfI]
                 )
@@ -966,6 +967,254 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
 }
 
 
+////XXXXXXXXX
+//// Get (pp-local) indices of points that are on both patches
+//Foam::labelList Foam::autoSnapDriver::getPatchSurfacePoints
+//(
+//    const fvMesh& mesh,
+//    const indirectPrimitivePatch& allPp,
+//    const polyPatch& pp
+//)
+//{
+//    // Could use PrimitivePatch & localFaces to extract points but might just
+//    // as well do it ourselves.
+//
+//    boolList pointOnZone(allPp.nPoints(), false);
+//
+//    forAll(pp, i)
+//    {
+//        const face& f = pp[i];
+//
+//        forAll(f, fp)
+//        {
+//            label meshPointI = f[fp];
+//
+//            Map<label>::const_iterator iter =
+//                allPp.meshPointMap().find(meshPointI);
+//
+//            if (iter != allPp.meshPointMap().end())
+//            {
+//                label pointI = iter();
+//                pointOnZone[pointI] = true;
+//            }
+//        }
+//    }
+//
+//    return findIndices(pointOnZone, true);
+//}
+//Foam::vectorField Foam::autoSnapDriver::calcNearestLocalSurface
+//(
+//    const meshRefinement& meshRefiner,
+//    const scalarField& snapDist,
+//    const indirectPrimitivePatch& pp
+//)
+//{
+//    Info<< "Calculating patchDisplacement as distance to nearest"
+//        << " local surface point ..." << endl;
+//
+//    const pointField& localPoints = pp.localPoints();
+//    const refinementSurfaces& surfaces = meshRefiner.surfaces();
+//    const fvMesh& mesh = meshRefiner.mesh();
+//    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+//
+//
+////    // Assume that all patch-internal points get attracted to their surface
+////    // only. So we want to know if point is on multiple regions
+////
+////    labelList minPatch(mesh.nPoints(), labelMax);
+////    labelList maxPatch(mesh.nPoints(), labelMin);
+////
+////    forAll(meshMover.adaptPatchIDs(), i)
+////    {
+////        label patchI = meshMover.adaptPatchIDs()[i];
+////        const labelList& meshPoints = pbm[patchI].meshPoints();
+////
+////        forAll(meshPoints, meshPointI)
+////        {
+////            label meshPointI = meshPoints[meshPointI];
+////            minPatch[meshPointI] = min(minPatch[meshPointI], patchI);
+////            maxPatch[meshPointI] = max(maxPatch[meshPointI], patchI);
+////        }
+////    }
+////
+////    syncTools::syncPointList
+////    (
+////        mesh,
+////        minPatch,
+////        minEqOp<label>(),   // combine op
+////        labelMax            // null value
+////    );
+////    syncTools::syncPointList
+////    (
+////        mesh,
+////        maxPatch,
+////        maxEqOp<label>(),   // combine op
+////        labelMin            // null value
+////    );
+//
+//    // Now all points with minPatch != maxPatch will be on the outside of
+//    // the patch.
+//
+//    // Displacement per patch point
+//    vectorField patchDisp(localPoints.size(), vector::zero);
+//    // Current best snap distance
+//    scalarField minSnapDist(snapDist);
+//    // Current surface snapped to
+//    labelList snapSurf(localPoints.size(), -1);
+//
+//    const labelList& surfaceGeometry = surfaces.surfaces();
+//    forAll(surfaceGeometry, surfI)
+//    {
+//        label geomI = surfaceGeometry[surfI];
+//        const wordList& regNames = allGeometry.regionNames()[geomI];
+//        forAll(regNames, regionI)
+//        {
+//            label globalRegionI = surfaces.globalRegion(surfI, regionI);
+//            // Collect master patch points
+//            label masterPatchI = globalToMasterPatch_[globalRegionI];
+//            label slavePatchI = globalToSlavePatch_[globalRegionI];
+//
+//            labelList patchPointIndices
+//            (
+//                getPatchSurfacePoints
+//                (
+//                    mesh,
+//                    pp,
+//                    pbm[masterPatchI]
+//                )
+//            );
+//
+//            // Find nearest for points both on faceZone and pp.
+//            List<pointIndexHit> hitInfo;
+//            surfaces.findNearest
+//            (
+//                surfI,
+//                regionI,
+//                pointField(localPoints, patchPointIndices),
+//                sqr(scalarField(minSnapDist, patchPointIndices)),
+//                hitInfo
+//            );
+//
+//            forAll(hitInfo, i)
+//            {
+//                label pointI = patchPointIndices[i];
+//
+//                if (hitInfo[i].hit())
+//                {
+//                    const point& pt = hitInfo[i].hitPoint();
+//                    patchDisp[pointI] = pt-localPoints[pointI];
+//                    minSnapDist[pointI] = min
+//                    (
+//                        minSnapDist[pointI],
+//                        mag(patchDisp[pointI])
+//                    );
+//                    snapSurf[pointI] = surfI;
+//                }
+//            }
+//
+//            // Slave patch
+//            if (slavePatchI != masterPatchI)
+//            {
+//                labelList patchPointIndices
+//                (
+//                    getPatchSurfacePoints
+//                    (
+//                        mesh,
+//                        pp,
+//                        pbm[slavePatchI]
+//                    )
+//                );
+//
+//                // Find nearest for points both on faceZone and pp.
+//                List<pointIndexHit> hitInfo;
+//                surfaces.findNearest
+//                (
+//                    surfI,
+//                    regionI,
+//                    pointField(localPoints, patchPointIndices),
+//                    sqr(scalarField(minSnapDist, patchPointIndices)),
+//                    hitInfo
+//                );
+//
+//                forAll(hitInfo, i)
+//                {
+//                    label pointI = patchPointIndices[i];
+//
+//                    if (hitInfo[i].hit())
+//                    {
+//                        const point& pt = hitInfo[i].hitPoint();
+//                        patchDisp[pointI] = pt-localPoints[pointI];
+//                        minSnapDist[pointI] = min
+//                        (
+//                            minSnapDist[pointI],
+//                            mag(patchDisp[pointI])
+//                        );
+//                        snapSurf[pointI] = surfI;
+//                    }
+//                }
+//            }
+//        }
+//    }
+//
+//
+//    // Check if all points are being snapped
+//    forAll(snapSurf, pointI)
+//    {
+//        if (snapSurf[pointI] == -1)
+//        {
+//            WarningIn("autoSnapDriver::calcNearestLocalSurface(..)")
+//                << "For point:" << pointI
+//                << " coordinate:" << localPoints[pointI]
+//                << " did not find any surface within:"
+//                << minSnapDist[pointI]
+//                << " metre." << endl;
+//        }
+//    }
+//
+//    {
+//        scalarField magDisp(mag(patchDisp));
+//
+//        Info<< "Wanted displacement : average:"
+//            << gSum(magDisp)/returnReduce(patchDisp.size(), sumOp<label>())
+//            << " min:" << gMin(magDisp)
+//            << " max:" << gMax(magDisp) << endl;
+//    }
+//
+//
+//    Info<< "Calculated surface displacement in = "
+//        << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
+//
+//
+//    // Limit amount of movement.
+//    forAll(patchDisp, patchPointI)
+//    {
+//        scalar magDisp = mag(patchDisp[patchPointI]);
+//
+//        if (magDisp > snapDist[patchPointI])
+//        {
+//            patchDisp[patchPointI] *= snapDist[patchPointI] / magDisp;
+//
+//            Pout<< "Limiting displacement for " << patchPointI
+//                << " from " << magDisp << " to " << snapDist[patchPointI]
+//                << endl;
+//        }
+//    }
+//
+//    // Points on zones in one domain but only present as point on other
+//    // will not do condition 2 on all. Sync explicitly.
+//    syncTools::syncPointList
+//    (
+//        mesh,
+//        pp.meshPoints(),
+//        patchDisp,
+//        minMagSqrEqOp<point>(),    // combine op
+//        vector(GREAT, GREAT, GREAT)// null value (note: cannot use VGREAT)
+//    );
+//
+//    return patchDisp;
+//}
+////XXXXXXXXX
+
 void Foam::autoSnapDriver::smoothDisplacement
 (
     const snapParameters& snapParams,
@@ -1153,7 +1402,15 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
         scalarField faceSnapDist(pp.size(), -GREAT);
         {
             // Distance to attract to nearest feature on surface
-            const scalarField snapDist(calcSnapDistance(snapParams, pp));
+            const scalarField snapDist
+            (
+                calcSnapDistance
+                (
+                    mesh,
+                    snapParams,
+                    pp
+                )
+            );
 
             const faceList& localFaces = pp.localFaces();
 
@@ -1429,7 +1686,7 @@ void Foam::autoSnapDriver::doSnap
         indirectPrimitivePatch& pp = ppPtr();
 
         // Distance to attract to nearest feature on surface
-        const scalarField snapDist(calcSnapDistance(snapParams, pp));
+        const scalarField snapDist(calcSnapDistance(mesh, snapParams, pp));
 
 
         // Construct iterative mesh mover.
@@ -1467,7 +1724,14 @@ void Foam::autoSnapDriver::doSnap
             << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
 
         // Pre-smooth patch vertices (so before determining nearest)
-        preSmoothPatch(snapParams, nInitErrors, baffles, meshMover);
+        preSmoothPatch
+        (
+            meshRefiner_,
+            snapParams,
+            nInitErrors,
+            baffles,
+            meshMover
+        );
 
 
         for (label iter = 0; iter < nFeatIter; iter++)
@@ -1478,7 +1742,12 @@ void Foam::autoSnapDriver::doSnap
 
             // Calculate displacement at every patch point. Insert into
             // meshMover.
-            vectorField disp = calcNearestSurface(snapDist, meshMover);
+            vectorField disp = calcNearestSurface
+            (
+                meshRefiner_,
+                snapDist,
+                meshMover.patch()
+            );
 
             // Override displacement with feature edge attempt
             if (doFeatures)

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,11 +81,11 @@ class autoSnapDriver
 
             //- Calculate displacement per patch point to smooth out patch.
             //  Quite complicated in determining which points to move where.
-            pointField smoothPatchDisplacement
+            static pointField smoothPatchDisplacement
             (
                 const motionSmoother&,
                 const List<labelPair>&
-            ) const;
+            );
 
             //- Check that face zones are synced
             void checkCoupledFaceZones() const;
@@ -406,37 +406,51 @@ public:
             autoPtr<mapPolyMesh> mergeZoneBaffles(const List<labelPair>&);
 
             //- Calculate edge length per patch point.
-            scalarField calcSnapDistance
+            static scalarField calcSnapDistance
             (
+                const fvMesh& mesh,
                 const snapParameters& snapParams,
                 const indirectPrimitivePatch&
-            ) const;
+            );
 
             //- Smooth the mesh (patch and internal) to increase visibility
             //  of surface points (on castellated mesh) w.r.t. surface.
-            void preSmoothPatch
+            static void preSmoothPatch
             (
+                const meshRefinement& meshRefiner,
                 const snapParameters& snapParams,
                 const label nInitErrors,
                 const List<labelPair>& baffles,
                 motionSmoother&
-            ) const;
+            );
 
             //- Get points both on patch and facezone.
-            labelList getZoneSurfacePoints
+            static labelList getZoneSurfacePoints
             (
+                const fvMesh& mesh,
                 const indirectPrimitivePatch&,
                 const word& zoneName
-            ) const;
+            );
 
             //- Per patch point calculate point on nearest surface. Set as
             //  boundary conditions of motionSmoother displacement field. Return
             //  displacement of patch points.
-            vectorField calcNearestSurface
+            static vectorField calcNearestSurface
             (
+                const meshRefinement& meshRefiner,
                 const scalarField& snapDist,
-                motionSmoother& meshMover
+                const indirectPrimitivePatch&
-            ) const;
+            );
+
+            ////- Per patch point calculate point on nearest surface. Set as
+            ////  boundary conditions of motionSmoother displacement field.
+            ////  Return displacement of patch points.
+            //static vectorField calcNearestLocalSurface
+            //(
+            //    const meshRefinement& meshRefiner,
+            //    const scalarField& snapDist,
+            //    const indirectPrimitivePatch&
+            //);
 
             //- Smooth the displacement field to the internal.
             void smoothDisplacement

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C

@@ -45,7 +45,11 @@ Foam::refinementParameters::refinementParameters
     nBufferLayers_(readLabel(dict.lookup("nBufferLayers"))),
     keepPoints_(dict.lookup("keepPoints")),
     allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")),
-    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0))
+    useTopologicalSnapDetection_
+    (
+        dict.lookupOrDefault<bool>("useTopologicalSnapDetection", true)
+    ),
+    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance", 0))
 {}
 
 
@@ -65,7 +69,11 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
     nBufferLayers_(readLabel(dict.lookup("nCellsBetweenLevels"))),
     keepPoints_(pointField(1, dict.lookup("locationInMesh"))),
     allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")),
-    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0))
+    useTopologicalSnapDetection_
+    (
+        dict.lookupOrDefault<bool>("useTopologicalSnapDetection", true)
+    ),
+    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance", 0))
 {
     scalar featAngle(readScalar(dict.lookup("resolveFeatureAngle")));
 
@@ -83,7 +91,7 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::labelList Foam::refinementParameters::findCells(const polyMesh& mesh)
- const
+const
 {
     // Force calculation of tet-diag decomposition (for use in findCell)
     (void)mesh.tetBasePtIs();

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H

@@ -80,6 +80,10 @@ class refinementParameters
         //  cellZone?
         Switch allowFreeStandingZoneFaces_;
 
+        //- Use old topology based problem-cell removal (cells with 8 points
+        //  on surface)
+        Switch useTopologicalSnapDetection_;
+
         //- Allowed load unbalance
         scalar maxLoadUnbalance_;
 
@@ -157,6 +161,13 @@ public:
                 return allowFreeStandingZoneFaces_;
             }
 
+            //- Use old topology based problem-cell removal
+            //  (cells with 8 points on surface)
+            bool useTopologicalSnapDetection() const
+            {
+                return useTopologicalSnapDetection_;
+            }
+
             //- Allowed load unbalance
             scalar maxLoadUnbalance() const
             {

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

@@ -72,6 +72,8 @@ class globalIndex;
 class removePoints;
 class localPointRegion;
 
+class snapParameters;
+
 /*---------------------------------------------------------------------------*\
                            Class meshRefinement Declaration
 \*---------------------------------------------------------------------------*/
@@ -470,6 +472,14 @@ private:
                 const labelList& globalToMasterPatch
             ) const;
 
+            //- Returns list with for every internal face -1 or the patch
+            //  they should be baffled into.
+            labelList markFacesOnProblemCellsGeometric
+            (
+                const snapParameters& snapParams,
+                const dictionary& motionDict
+            ) const;
+
 
         // Baffle merging
 
@@ -536,6 +546,8 @@ private:
             //- Remove any loose standing cells
             void handleSnapProblems
             (
+                const snapParameters& snapParams,
+                const bool useTopologicalSnapDetection,
                 const bool removeEdgeConnectedCells,
                 const scalarField& perpendicularAngle,
                 const dictionary& motionDict,
@@ -760,10 +772,17 @@ public:
             void baffleAndSplitMesh
             (
                 const bool handleSnapProblems,
+
+                // How to remove problem snaps
+                const snapParameters& snapParams,
+                const bool useTopologicalSnapDetection,
                 const bool removeEdgeConnectedCells,
                 const scalarField& perpendicularAngle,
+
+                // How to handle free-standing baffles
                 const bool mergeFreeStanding,
                 const scalar freeStandingAngle,
+
                 const dictionary& motionDict,
                 Time& runTime,
                 const labelList& globalToMasterPatch,

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

@@ -1872,6 +1872,8 @@ void Foam::meshRefinement::makeConsistentFaceIndex
 
 void Foam::meshRefinement::handleSnapProblems
 (
+    const snapParameters& snapParams,
+    const bool useTopologicalSnapDetection,
     const bool removeEdgeConnectedCells,
     const scalarField& perpendicularAngle,
     const dictionary& motionDict,
@@ -1885,44 +1887,39 @@ void Foam::meshRefinement::handleSnapProblems
         << "----------------------------------------------" << nl
         << endl;
 
-    labelList facePatch
+    labelList facePatch;
-    (
+    if (useTopologicalSnapDetection)
-        markFacesOnProblemCells
+    {
+        facePatch = markFacesOnProblemCells
         (
             motionDict,
             removeEdgeConnectedCells,
             perpendicularAngle,
             globalToMasterPatch
-        )
+        );
-    );
+    }
+    else
+    {
+        facePatch = markFacesOnProblemCellsGeometric(snapParams, motionDict);
+    }
     Info<< "Analyzed problem cells in = "
         << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
 
     if (debug&meshRefinement::MESH)
     {
-        faceSet problemTopo(mesh_, "problemFacesTopo", 100);
+        faceSet problemFaces(mesh_, "problemFaces", 100);
 
-        const labelList facePatchTopo
+        forAll(facePatch, faceI)
-        (
-            markFacesOnProblemCells
-            (
-                motionDict,
-                removeEdgeConnectedCells,
-                perpendicularAngle,
-                globalToMasterPatch
-            )
-        );
-        forAll(facePatchTopo, faceI)
         {
-            if (facePatchTopo[faceI] != -1)
+            if (facePatch[faceI] != -1)
             {
-                problemTopo.insert(faceI);
+                problemFaces.insert(faceI);
             }
         }
-        problemTopo.instance() = timeName();
+        problemFaces.instance() = timeName();
-        Pout<< "Dumping " << problemTopo.size()
+        Pout<< "Dumping " << problemFaces.size()
-            << " problem faces to " << problemTopo.objectPath() << endl;
+            << " problem faces to " << problemFaces.objectPath() << endl;
-        problemTopo.write();
+        problemFaces.write();
     }
 
     Info<< "Introducing baffles to delete problem cells." << nl << endl;
@@ -1962,6 +1959,8 @@ void Foam::meshRefinement::handleSnapProblems
 void Foam::meshRefinement::baffleAndSplitMesh
 (
     const bool doHandleSnapProblems,
+    const snapParameters& snapParams,
+    const bool useTopologicalSnapDetection,
     const bool removeEdgeConnectedCells,
     const scalarField& perpendicularAngle,
     const bool mergeFreeStanding,
@@ -2029,83 +2028,10 @@ void Foam::meshRefinement::baffleAndSplitMesh
 
     if (doHandleSnapProblems)
     {
-        //Info<< nl
-        //    << "Introducing baffles to block off problem cells" << nl
-        //    << "----------------------------------------------" << nl
-        //    << endl;
-        //
-        //labelList facePatch
-        //(
-        //    markFacesOnProblemCells
-        //    (
-        //        motionDict,
-        //        removeEdgeConnectedCells,
-        //        perpendicularAngle,
-        //        globalToMasterPatch
-        //    )
-        //    //markFacesOnProblemCellsGeometric(motionDict)
-        //);
-        //Info<< "Analyzed problem cells in = "
-        //    << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
-        //
-        //if (debug&meshRefinement::MESH)
-        //{
-        //    faceSet problemTopo(mesh_, "problemFacesTopo", 100);
-        //
-        //    const labelList facePatchTopo
-        //    (
-        //        markFacesOnProblemCells
-        //        (
-        //            motionDict,
-        //            removeEdgeConnectedCells,
-        //            perpendicularAngle,
-        //            globalToMasterPatch
-        //        )
-        //    );
-        //    forAll(facePatchTopo, faceI)
-        //    {
-        //        if (facePatchTopo[faceI] != -1)
-        //        {
-        //            problemTopo.insert(faceI);
-        //        }
-        //    }
-        //    problemTopo.instance() = timeName();
-        //    Pout<< "Dumping " << problemTopo.size()
-        //        << " problem faces to " << problemTopo.objectPath() << endl;
-        //    problemTopo.write();
-        //}
-        //
-        //Info<< "Introducing baffles to delete problem cells." << nl << endl;
-        //
-        //if (debug)
-        //{
-        //    runTime++;
-        //}
-        //
-        //// Create baffles with same owner and neighbour for now.
-        //createBaffles(facePatch, facePatch);
-        //
-        //if (debug)
-        //{
-        //    // Debug:test all is still synced across proc patches
-        //    checkData();
-        //}
-        //Info<< "Created baffles in = "
-        //    << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
-        //
-        //printMeshInfo(debug, "After introducing baffles");
-        //
-        //if (debug&meshRefinement::MESH)
-        //{
-        //    Pout<< "Writing extra baffled mesh to time "
-        //        << timeName() << endl;
-        //    write(debug, runTime.path()/"extraBaffles");
-        //    Pout<< "Dumped debug data in = "
-        //        << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
-        //}
-
         handleSnapProblems
         (
+            snapParams,
+            useTopologicalSnapDetection,
             removeEdgeConnectedCells,
             perpendicularAngle,
             motionDict,
@@ -2192,6 +2118,8 @@ void Foam::meshRefinement::baffleAndSplitMesh
             // and delete them
             handleSnapProblems
             (
+                snapParams,
+                useTopologicalSnapDetection,
                 removeEdgeConnectedCells,
                 perpendicularAngle,
                 motionDict,

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

@@ -36,6 +36,10 @@ License
 #include "polyMeshGeometry.H"
 #include "IOmanip.H"
 #include "unitConversion.H"
+#include "autoSnapDriver.H"
+
+#include "snapParameters.H"
+#include "motionSmoother.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -968,164 +972,241 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 }
 
 
-//// Mark faces to be baffled to prevent snapping problems. Does
+// Mark faces to be baffled to prevent snapping problems. Does
-//// test to find nearest surface and checks which faces would get squashed.
+// test to find nearest surface and checks which faces would get squashed.
-//Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
+Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
-//(
+(
-//    const dictionary& motionDict
+    const snapParameters& snapParams,
-//) const
+    const dictionary& motionDict
-//{
+) const
-//    // Construct addressing engine.
+{
-//    autoPtr<indirectPrimitivePatch> ppPtr
+    pointField oldPoints(mesh_.points());
-//    (
+
-//        meshRefinement::makePatch
+    // Repeat (most of) autoSnapDriver::doSnap
-//        (
+    {
-//            mesh_,
+        labelList adaptPatchIDs(meshedPatches());
-//            meshedPatches()
+
-//        )
+        // Construct addressing engine.
-//    );
+        autoPtr<indirectPrimitivePatch> ppPtr
-//    const indirectPrimitivePatch& pp = ppPtr();
+        (
-//    const pointField& localPoints = pp.localPoints();
+            meshRefinement::makePatch
-//    const labelList& meshPoints = pp.meshPoints();
+            (
-//
+                mesh_,
-//    // Find nearest (non-baffle) surface
+                adaptPatchIDs
-//    pointField newPoints(mesh_.points());
+            )
-//    {
+        );
-//        List<pointIndexHit> hitInfo;
+        indirectPrimitivePatch& pp = ppPtr();
-//        labelList hitSurface;
+
-//        surfaces_.findNearest
+        // Distance to attract to nearest feature on surface
-//        (
+        const scalarField snapDist
-//            surfaces_.getUnnamedSurfaces(),
+        (
-//            localPoints,
+            autoSnapDriver::calcSnapDistance(mesh_, snapParams, pp)
-//            scalarField(localPoints.size(), sqr(GREAT)),// sqr of attraction
+        );
-//            hitSurface,
+
-//            hitInfo
+
-//        );
+        // Construct iterative mesh mover.
-//
+        Info<< "Constructing mesh displacer ..." << endl;
-//        forAll(hitInfo, i)
+        Info<< "Using mesh parameters " << motionDict << nl << endl;
-//        {
+
-//            if (hitInfo[i].hit())
+        const pointMesh& pMesh = pointMesh::New(mesh_);
-//            {
+
-//                //label pointI = meshPoints[i];
+        motionSmoother meshMover
-//                //Pout<< "   " << pointI << " moved from "
+        (
-//                //    << mesh_.points()[pointI] << " by "
+            mesh_,
-//                //    << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
+            pp,
-//                //    << endl;
+            adaptPatchIDs,
-//                newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
+            meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs)(),
-//            }
+            motionDict
-//        }
+        );
-//    }
+
-//
+
-//    // Per face (internal or coupled!) the patch that the
+        // Check initial mesh
-//    // baffle should get (or -1).
+        Info<< "Checking initial mesh ..." << endl;
-//    labelList facePatch(mesh_.nFaces(), -1);
+        labelHashSet wrongFaces(mesh_.nFaces()/100);
-//    // Count of baffled faces
+        motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
-//    label nBaffleFaces = 0;
+        const label nInitErrors = returnReduce
-//
+        (
-//    {
+            wrongFaces.size(),
-//        pointField oldPoints(mesh_.points());
+            sumOp<label>()
-//        mesh_.movePoints(newPoints);
+        );
-//        faceSet wrongFaces(mesh_, "wrongFaces", 100);
+
-//        {
+        Info<< "Detected " << nInitErrors << " illegal faces"
-//            //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
+            << " (concave, zero area or negative cell pyramid volume)"
-//
+            << endl;
-//            // Just check the errors from squashing
+
-//            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
-//
+        Info<< "Checked initial mesh in = "
-//            const labelList allFaces(identity(mesh_.nFaces()));
+            << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
-//            label nWrongFaces = 0;
+
-//
+        // Pre-smooth patch vertices (so before determining nearest)
-//            scalar minArea(readScalar(motionDict.lookup("minArea")));
+        autoSnapDriver::preSmoothPatch
-//            if (minArea > -SMALL)
+        (
-//            {
+            *this,
-//                polyMeshGeometry::checkFaceArea
+            snapParams,
-//                (
+            nInitErrors,
-//                    false,
+            List<labelPair>(0), //baffles
-//                    minArea,
+            meshMover
-//                    mesh_,
+        );
-//                    mesh_.faceAreas(),
+
-//                    allFaces,
+        const vectorField disp
-//                    &wrongFaces
+        (
-//                );
+            autoSnapDriver::calcNearestSurface
-//
+            (
-//                label nNewWrongFaces = returnReduce
+                *this,
-//                (
+                snapDist,   // attraction
-//                    wrongFaces.size(),
+                pp
-//                    sumOp<label>()
+            )
-//                );
+        );
-//
+
-//                Info<< "    faces with area < "
+        const labelList& meshPoints = pp.meshPoints();
-//                    << setw(5) << minArea
+
-//                    << " m^2                            : "
+        pointField newPoints(mesh_.points());
-//                    << nNewWrongFaces-nWrongFaces << endl;
+        forAll(meshPoints, i)
-//
+        {
-//                nWrongFaces = nNewWrongFaces;
+            newPoints[meshPoints[i]] += disp[i];
-//            }
+        }
-//
+        mesh_.movePoints(newPoints);
-////            scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
+    }
-//            scalar minDet = 0.01;
+
-//            if (minDet > -1)
+
-//            {
+    // Per face (internal or coupled!) the patch that the
-//                polyMeshGeometry::checkCellDeterminant
+    // baffle should get (or -1).
-//                (
+    labelList facePatch(mesh_.nFaces(), -1);
-//                    false,
+    // Count of baffled faces
-//                    minDet,
+    label nBaffleFaces = 0;
-//                    mesh_,
+
-//                    mesh_.faceAreas(),
+    {
-//                    allFaces,
+        faceSet wrongFaces(mesh_, "wrongFaces", 100);
-//                    polyMeshGeometry::affectedCells(mesh_, allFaces),
+        {
-//                    &wrongFaces
+            //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
-//                );
+
-//
+            // Just check the errors from squashing
-//                label nNewWrongFaces = returnReduce
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-//                (
+
-//                    wrongFaces.size(),
+            const labelList allFaces(identity(mesh_.nFaces()));
-//                    sumOp<label>()
+            label nWrongFaces = 0;
-//                );
+
-//
+            //const scalar minV(readScalar(motionDict.lookup("minVol", true)));
-//                Info<< "    faces on cells with determinant < "
+            //if (minV > -GREAT)
-//                    << setw(5) << minDet << "                : "
+            //{
-//                    << nNewWrongFaces-nWrongFaces << endl;
+            //    polyMeshGeometry::checkFacePyramids
-//
+            //    (
-//                nWrongFaces = nNewWrongFaces;
+            //        false,
-//            }
+            //        minV,
-//        }
+            //        mesh_,
-//
+            //        mesh_.cellCentres(),
-//
+            //        mesh_.points(),
-//        forAllConstIter(faceSet, wrongFaces, iter)
+            //        allFaces,
-//        {
+            //        List<labelPair>(0),
-//            label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
+            //        &wrongFaces
-//
+            //    );
-//            if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
+            //
-//            {
+            //    label nNewWrongFaces = returnReduce
-//                facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
+            //    (
-//                nBaffleFaces++;
+            //        wrongFaces.size(),
-//
+            //        sumOp<label>()
-//                //Pout<< "    " << iter.key()
+            //    );
-//                //    //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
+            //
-//                //    << " is destined for patch " << facePatch[iter.key()]
+            //    Info<< "    faces with pyramid volume < "
-//                //    << endl;
+            //        << setw(5) << minV
-//            }
+            //        << " m^3                  : "
-//        }
+            //        << nNewWrongFaces-nWrongFaces << endl;
-//        // Restore points.
+            //
-//        mesh_.movePoints(oldPoints);
+            //    nWrongFaces = nNewWrongFaces;
-//    }
+            //}
-//
+
-//
+            scalar minArea(readScalar(motionDict.lookup("minArea")));
-//    Info<< "markFacesOnProblemCellsGeometric : marked "
+            if (minArea > -SMALL)
-//        << returnReduce(nBaffleFaces, sumOp<label>())
+            {
-//        << " additional internal and coupled faces"
+                polyMeshGeometry::checkFaceArea
-//        << " to be converted into baffles." << endl;
+                (
-//
+                    false,
-//    syncTools::syncFaceList
+                    minArea,
-//    (
+                    mesh_,
-//        mesh_,
+                    mesh_.faceAreas(),
-//        facePatch,
+                    allFaces,
-//        maxEqOp<label>()
+                    &wrongFaces
-//    );
+                );
-//
+
-//    return facePatch;
+                label nNewWrongFaces = returnReduce
-//}
+                (
+                    wrongFaces.size(),
+                    sumOp<label>()
+                );
+
+                Info<< "    faces with area < "
+                    << setw(5) << minArea
+                    << " m^2                            : "
+                    << nNewWrongFaces-nWrongFaces << endl;
+
+                nWrongFaces = nNewWrongFaces;
+            }
+
+            scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
+            if (minDet > -1)
+            {
+                polyMeshGeometry::checkCellDeterminant
+                (
+                    false,
+                    minDet,
+                    mesh_,
+                    mesh_.faceAreas(),
+                    allFaces,
+                    polyMeshGeometry::affectedCells(mesh_, allFaces),
+                    &wrongFaces
+                );
+
+                label nNewWrongFaces = returnReduce
+                (
+                    wrongFaces.size(),
+                    sumOp<label>()
+                );
+
+                Info<< "    faces on cells with determinant < "
+                    << setw(5) << minDet << "                : "
+                    << nNewWrongFaces-nWrongFaces << endl;
+
+                nWrongFaces = nNewWrongFaces;
+            }
+        }
+
+
+        forAllConstIter(faceSet, wrongFaces, iter)
+        {
+            label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
+
+            if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
+            {
+                facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
+                nBaffleFaces++;
+
+                //Pout<< "    " << iter.key()
+                //    //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
+                //    << " is destined for patch " << facePatch[iter.key()]
+                //    << endl;
+            }
+        }
+    }
+
+
+    // Restore points.
+    mesh_.movePoints(oldPoints);
+
+
+    Info<< "markFacesOnProblemCellsGeometric : marked "
+        << returnReduce(nBaffleFaces, sumOp<label>())
+        << " additional internal and coupled faces"
+        << " to be converted into baffles." << endl;
+
+    syncTools::syncFaceList
+    (
+        mesh_,
+        facePatch,
+        maxEqOp<label>()
+    );
+
+    return facePatch;
+}
 
 
 // ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -77,7 +77,8 @@ const Foam::NamedEnum<Foam::refinementSurfaces::faceZoneType, 3>
 Foam::refinementSurfaces::refinementSurfaces
 (
     const searchableSurfaces& allGeometry,
-    const dictionary& surfacesDict
+    const dictionary& surfacesDict,
+    const label gapLevelIncrement
 )
 :
     allGeometry_(allGeometry),
@@ -143,7 +144,7 @@ Foam::refinementSurfaces::refinementSurfaces
             globalLevelIncr[surfI] = dict.lookupOrDefault
             (
                 "gapLevelIncrement",
-                0
+                gapLevelIncrement
             );
 
             if
@@ -274,7 +275,7 @@ Foam::refinementSurfaces::refinementSurfaces
                         label levelIncr = regionDict.lookupOrDefault
                         (
                             "gapLevelIncrement",
-                            0
+                            gapLevelIncrement
                         );
                         regionLevelIncr[surfI].insert(regionI, levelIncr);
 

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H

@@ -149,7 +149,8 @@ public:
         refinementSurfaces
         (
             const searchableSurfaces& allGeometry,
-            const dictionary&
+            const dictionary&,
+            const label gapLevelIncrement
         );