problem cell deletion

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -131,6 +131,12 @@ castellatedMeshControls
                     level (3 3);
                 }
             }
+
+            // Optional angle to detect small-large cell situation perpendicular
+            // to the surface. Is the angle of face w.r.t the local surface
+            // normal. Use on flat(ish) surfaces only. Otherwise
+            // leave out or set to negative number.
+            //perpendicularAngle 10;
         }
     }
 
@@ -317,9 +323,9 @@ meshQualityControls
     minTriangleTwist -1;
 
     //- if >0 : preserve single cells with all points on the surface if the
-    //  resulting volume after snapping is larger than minVolFraction times old
+    //  resulting volume after snapping (approximation) is larger than
-    //  volume. If <0 : delete always.
+    //  minVolFraction times old volume. If <0 : delete always.
-    minVolFraction 0.1;
+    minVolFraction 0.5;
 
 
     // Advanced

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C

@@ -58,35 +58,6 @@ defineTypeNameAndDebug(autoSnapDriver, 0);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::autoSnapDriver::getZonedSurfaces
-(
-    labelList& zonedSurfaces,
-    labelList& unzonedSurfaces
-) const
-{    // Surfaces with zone information
-    const wordList& faceZoneNames = meshRefiner_.surfaces().faceZoneNames();
-
-    zonedSurfaces.setSize(faceZoneNames.size());
-    label zonedI = 0;
-    unzonedSurfaces.setSize(faceZoneNames.size());
-    label unzonedI = 0;
-
-    forAll(faceZoneNames, surfI)
-    {
-        if (faceZoneNames[surfI].size())
-        {
-            zonedSurfaces[zonedI++] = surfI;
-        }
-        else
-        {
-            unzonedSurfaces[unzonedI++] = surfI;
-        }
-    }
-    zonedSurfaces.setSize(zonedI);
-    unzonedSurfaces.setSize(unzonedI);
-}
-
-
 // Get faces to repatch. Returns map from face to patch.
 Foam::Map<Foam::label> Foam::autoSnapDriver::getZoneBafflePatches
 (
@@ -711,9 +682,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::createZoneBaffles
     List<labelPair>& baffles
 )
 {
-    labelList zonedSurfaces;
+    labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
-    labelList unzonedSurfaces;
-    getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
 
     autoPtr<mapPolyMesh> map;
 
@@ -798,9 +767,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::mergeZoneBaffles
     const List<labelPair>& baffles
 )
 {
-    labelList zonedSurfaces;
+    labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
-    labelList unzonedSurfaces;
-    getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
 
     autoPtr<mapPolyMesh> map;
 
@@ -1048,9 +1015,10 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
         labelList snapSurf(localPoints.size(), -1);
 
         // Divide surfaces into zoned and unzoned
-        labelList zonedSurfaces;
+        labelList zonedSurfaces =
-        labelList unzonedSurfaces;
+            meshRefiner_.surfaces().getNamedSurfaces();
-        getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
+        labelList unzonedSurfaces =
+            meshRefiner_.surfaces().getUnnamedSurfaces();
 
 
         // 1. All points to non-interface surfaces
@@ -1368,9 +1336,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
     indirectPrimitivePatch& pp = ppPtr();
 
     // Divide surfaces into zoned and unzoned
-    labelList zonedSurfaces;
+    labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
-    labelList unzonedSurfaces;
+    labelList unzonedSurfaces = meshRefiner_.surfaces().getUnnamedSurfaces();
-    getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
 
 
     // Faces that do not move

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H

@@ -84,9 +84,6 @@ class autoSnapDriver
 
         // Snapping
 
-             //- Split surfaces into non-zoned and zones ones
-            void getZonedSurfaces(labelList&, labelList&) const;
-
             //- Get faces to repatch. Returns map from face to patch.
             Map<label> getZoneBafflePatches(const bool allowBoundary) const;
 

src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H

@@ -374,11 +374,28 @@ private:
                 const labelList&
             ) const;
 
+            bool isCollapsedFace
+            (
+                const pointField&,
+                const pointField& neiCc,
+                const scalar minFaceArea,
+                const scalar maxNonOrtho,
+                const label faceI
+            ) const;
+
+            bool isCollapsedCell
+            (
+                const pointField&,
+                const scalar volFraction,
+                const label cellI
+            ) const;
+
             //- Returns list with for every internal face -1 or the patch
             //  they should be baffled into. If removeEdgeConnectedCells is set
             //  removes cells based on perpendicularAngle.
             labelList markFacesOnProblemCells
             (
+                const dictionary& motionDict,
                 const bool removeEdgeConnectedCells,
                 const scalarField& perpendicularAngle,
                 const labelList& globalToPatch

src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C

@@ -1503,6 +1503,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
         (
             markFacesOnProblemCells
             (
+                motionDict,
                 removeEdgeConnectedCells,
                 perpendicularAngle,
                 globalToPatch
@@ -1548,6 +1549,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
             (
                 markFacesOnProblemCells
                 (
+                    motionDict,
                     removeEdgeConnectedCells,
                     perpendicularAngle,
                     globalToPatch

src/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C

@@ -278,13 +278,105 @@ Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
 }
 
 
+// Check if moving face to new points causes a 'collapsed' face.
+// Uses new point position only for the face, not the neighbouring
+// cell centres
+bool Foam::meshRefinement::isCollapsedFace
+(
+    const pointField& points,
+    const pointField& neiCc,
+    const scalar minFaceArea,
+    const scalar maxNonOrtho,
+    const label faceI
+) const
+{
+    vector s = mesh_.faces()[faceI].normal(points);
+    scalar magS = mag(s);
+
+    // Check face area
+    if (magS < minFaceArea)
+    {
+        return true;
+    }
+
+    // Check orthogonality
+    const point& ownCc = mesh_.cellCentres()[mesh_.faceOwner()[faceI]];
+
+    if (mesh_.isInternalFace(faceI))
+    {
+        label nei = mesh_.faceNeighbour()[faceI];
+        vector d = ownCc - mesh_.cellCentres()[nei];
+
+        scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
+        
+        if (dDotS < maxNonOrtho)
+        {
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+    else
+    {
+        label patchI = mesh_.boundaryMesh().whichPatch(faceI);
+
+        if (mesh_.boundaryMesh()[patchI].coupled())
+        {
+            vector d = ownCc - neiCc[faceI-mesh_.nInternalFaces()];
+
+            scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
+
+            if (dDotS < maxNonOrtho)
+            {
+                return true;
+            }
+            else
+            {
+                return false;
+            }
+        }
+        else
+        {
+            // Collapsing normal boundary face does not cause problems with
+            // non-orthogonality
+            return true;
+        }
+    }
+}
+
+
+// Check if moving cell to new points causes it to collapse.
+bool Foam::meshRefinement::isCollapsedCell
+(
+    const pointField& points,
+    const scalar volFraction,
+    const label cellI
+) const
+{
+    scalar vol = mesh_.cells()[cellI].mag(points, mesh_.faces());
+
+    if (vol/mesh_.cellVolumes()[cellI] < volFraction)
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 // Returns list with for every internal face -1 or the patch they should
-// be baffled into. Gets run after createBaffles so all the surface
+// be baffled into. Gets run after createBaffles so all the unzoned surface
-// intersections have already been turned into baffles. Used to remove cells
+// intersections have already been turned into baffles. (Note: zoned surfaces
-// by baffling all their faces and have the splitMeshRegions chuck away non
+// are not baffled at this stage)
-// used regions.
+// Used to remove cells by baffling all their faces and have the
+// splitMeshRegions chuck away non used regions.
 Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 (
+    const dictionary& motionDict,
     const bool removeEdgeConnectedCells,
     const scalarField& perpendicularAngle,
     const labelList& globalToPatch
@@ -294,6 +386,10 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     const labelList& pointLevel = meshCutter_.pointLevel();
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
+    // Swap neighbouring cell centres and cell level
+    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
+    pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
+    calcNeighbourData(neiLevel, neiCc);
 
     // Per internal face (boundary faces not used) the patch that the
     // baffle should get (or -1)
@@ -334,6 +430,9 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     // Count of faces marked for baffling
     label nBaffleFaces = 0;
 
+    // Count of faces not baffled since would not cause a collapse
+    label nPrevented = 0;
+
     if (removeEdgeConnectedCells && max(perpendicularAngle) >= 0)
     {
         Info<< "markFacesOnProblemCells :"
@@ -420,6 +519,74 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     );
 
 
+    // See if checking for collapse
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Collapse checking parameters
+    scalar volFraction = -1;
+    if (motionDict.found("minVolFraction"))
+    {
+        volFraction = readScalar(motionDict.lookup("minVolFraction"));
+    }
+    const bool checkCollapse = (volFraction > 0);
+    scalar minArea = -1;
+    scalar maxNonOrtho = -1;
+
+
+    // Find nearest (non-baffle) surface
+    pointField newPoints;
+
+    if (checkCollapse)
+    {
+        minArea = readScalar(motionDict.lookup("minArea"));
+        maxNonOrtho = readScalar(motionDict.lookup("maxNonOrtho"));
+
+        Info<< "markFacesOnProblemCells :"
+            << " Deleting all-anchor surface cells only if"
+            << "snapping them violates mesh quality constraints:" << nl
+            << "    snapped/original cell volume < " << volFraction << nl
+            << "    face area                    < " << minArea << nl
+            << "    non-orthogonality            > " << maxNonOrtho << nl
+            << endl;
+
+        // Construct addressing engine.
+        autoPtr<indirectPrimitivePatch> ppPtr
+        (
+            meshRefinement::makePatch
+            (
+                mesh_,
+                adaptPatchIDs
+            )
+        );
+        const indirectPrimitivePatch& pp = ppPtr();
+        const pointField& localPoints = pp.localPoints();
+        const labelList& meshPoints = pp.meshPoints();
+
+        List<pointIndexHit> hitInfo;
+        labelList hitSurface;
+        surfaces_.findNearest
+        (
+            surfaces_.getUnnamedSurfaces(),
+            localPoints,
+            scalarField(localPoints.size(), sqr(GREAT)),    // sqr of attraction
+            hitSurface,
+            hitInfo
+        );
+
+        // Start of from current points
+        newPoints = mesh_.points();
+
+        forAll(hitInfo, i)
+        {
+            if (hitInfo[i].hit())
+            {
+                newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
+            }
+        }
+    }
+
+
+
     // For each cell count the number of anchor points that are on
     // the boundary:
     // 8 : check the number of (baffle) boundary faces. If 3 or more block
@@ -430,6 +597,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     //     3 or more of these cells. (note that on a flat surface a boundary
     //     point will only have 4 cells connected to it)
 
+
     // Does cell have exactly 7 of its 8 anchor points on the boundary?
     PackedBoolList hasSevenBoundaryAnchorPoints(mesh_.nCells());
     // If so what is the remaining non-boundary anchor point?
@@ -494,23 +662,45 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
             // Eugene: delete cell no matter what.
             //if (nBfaces > 1)
             {
-                forAll(cFaces, cf)
+                if
+                (
+                    checkCollapse
+                && !isCollapsedCell(newPoints, volFraction, cellI)
+                )
                 {
-                    label faceI = cFaces[cf];
+                    nPrevented++;
-
+                    //Pout<< "Preventing collapse of 8 anchor point cell "
-                    if (facePatch[faceI] == -1 && mesh_.isInternalFace(faceI))
+                    //    << cellI << " at " << mesh_.cellCentres()[cellI]
+                    //    << " since new volume "
+                    //    << mesh_.cells()[cellI].mag(newPoints, mesh_.faces())
+                    //    << " old volume " << mesh_.cellVolumes()[cellI]
+                    //    << endl;
+                }
+                else
+                {
+                    // Block all faces of cell
+                    forAll(cFaces, cf)
                     {
-                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                        label faceI = cFaces[cf];
-                        nBaffleFaces++;
 
-                        // Mark face as a 'boundary'
+                        if
-                        markBoundaryFace
                         (
-                            faceI,
+                            facePatch[faceI] == -1
-                            isBoundaryFace,
+                         && mesh_.isInternalFace(faceI)
-                            isBoundaryEdge,
+                        )
-                            isBoundaryPoint
+                        {
-                        );
+                            facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                            nBaffleFaces++;
+
+                            // Mark face as a 'boundary'
+                            markBoundaryFace
+                            (
+                                faceI,
+                                isBoundaryFace,
+                                isBoundaryEdge,
+                                isBoundaryPoint
+                            );
+                        }
                     }
                 }
             }
@@ -556,29 +746,52 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 
                 if (hasSevenBoundaryAnchorPoints.get(cellI) == 1u)
                 {
-                    const cell& cFaces = mesh_.cells()[cellI];
+                    if
-
+                    (
-                    forAll(cFaces, cf)
+                        checkCollapse
+                    && !isCollapsedCell(newPoints, volFraction, cellI)
+                    )
                     {
-                        label faceI = cFaces[cf];
+                        nPrevented++;
+                        //Pout<< "Preventing collapse of 7 anchor cell "
+                        //    << cellI
+                        //    << " at " << mesh_.cellCentres()[cellI]
+                        //    << " since new volume "
+                        //    << mesh_.cells()[cellI].mag
+                        //        (newPoints, mesh_.faces())
+                        //    << " old volume " << mesh_.cellVolumes()[cellI]
+                        //    << endl;
+                    }
+                    else
+                    {
+                        const cell& cFaces = mesh_.cells()[cellI];
 
-                        if
+                        forAll(cFaces, cf)
-                        (
-                            facePatch[faceI] == -1
-                         && mesh_.isInternalFace(faceI)
-                        )
                         {
-                            facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                            label faceI = cFaces[cf];
-                            nBaffleFaces++;
 
-                            // Mark face as a 'boundary'
+                            if
-                            markBoundaryFace
                             (
-                                faceI,
+                                facePatch[faceI] == -1
-                                isBoundaryFace,
+                             && mesh_.isInternalFace(faceI)
-                                isBoundaryEdge,
+                            )
-                                isBoundaryPoint
+                            {
-                            );
+                                facePatch[faceI] = getBafflePatch
+                                (
+                                    facePatch,
+                                    faceI
+                                );
+                                nBaffleFaces++;
+
+                                // Mark face as a 'boundary'
+                                markBoundaryFace
+                                (
+                                    faceI,
+                                    isBoundaryFace,
+                                    isBoundaryEdge,
+                                    isBoundaryPoint
+                                );
+                            }
                         }
                     }
                 }
@@ -635,11 +848,33 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 
             if (nFaceBoundaryEdges == fEdges.size())
             {
-                facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                if
-                nBaffleFaces++;
+                (
+                    checkCollapse
+                && !isCollapsedFace
+                    (
+                        newPoints,
+                        neiCc,
+                        minArea,
+                        maxNonOrtho,
+                        faceI
+                    )
+                )
+                {
+                    nPrevented++;
+                    //Pout<< "Preventing collapse of face " << faceI
+                    //    << " with all boundary edges "
+                    //    << " at " << mesh_.faceCentres()[faceI]
+                    //    << endl;
+                }
+                else
+                {
+                    facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                    nBaffleFaces++;
 
-                // Do NOT update boundary data since this would grow blocked
+                    // Do NOT update boundary data since this would grow blocked
-                // faces across gaps.
+                    // faces across gaps.
+                }
             }
         }
     }
@@ -669,11 +904,34 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 
                     if (nFaceBoundaryEdges == fEdges.size())
                     {
-                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                        if
-                        nBaffleFaces++;
+                        (
+                            checkCollapse
+                        && !isCollapsedFace
+                            (
+                                newPoints,
+                                neiCc,
+                                minArea,
+                                maxNonOrtho,
+                                faceI
+                            )
+                        )
+                        {
+                            nPrevented++;
+                            //Pout<< "Preventing collapse of coupled face "
+                            //    << faceI
+                            //    << " with all boundary edges "
+                            //    << " at " << mesh_.faceCentres()[faceI]
+                            //    << endl;
+                        }
+                        else
+                        {
+                            facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                            nBaffleFaces++;
 
-                        // Do NOT update boundary data since this would grow
+                            // Do NOT update boundary data since this would grow
-                        // blocked faces across gaps.
+                            // blocked faces across gaps.
+                        }
                     }
                 }
 
@@ -687,266 +945,181 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
         << " additional internal faces to be converted into baffles."
         << endl;
 
-    return facePatch;
+    if (checkCollapse)
-}
-//XXXXXXXXXXXXXX
-// Mark faces to be baffled to prevent snapping problems. Does
-// test to find nearest surface and checks which faces would get squashed.
-Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
-(
-    const dictionary& motionDict,
-    const labelList& globalToPatch
-) const
-{
-    // Get the labels of added patches.
-    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
-
-    // Construct addressing engine.
-    autoPtr<indirectPrimitivePatch> ppPtr
-    (
-        meshRefinement::makePatch
-        (
-            mesh_,
-            adaptPatchIDs
-        )
-    );
-    const indirectPrimitivePatch& pp = ppPtr();
-    const pointField& localPoints = pp.localPoints();
-    const labelList& meshPoints = pp.meshPoints();
-
-    // Find nearest (non-baffle) surface
-    pointField newPoints(mesh_.points());
     {
-        List<pointIndexHit> hitInfo;
+        Info<< "markFacesOnProblemCells : prevented "
-        labelList hitSurface;
+            << returnReduce(nPrevented, sumOp<label>())
-        surfaces_.findNearest
+            << " internal faces fom getting converted into baffles."
-        (
+            << endl;
-            surfaces_.getUnnamedSurfaces(),
-            localPoints,
-            scalarField(localPoints.size(), sqr(GREAT)),    // sqr of attraction
-            hitSurface,
-            hitInfo
-        );
-
-        forAll(hitInfo, i)
-        {
-            if (hitInfo[i].hit())
-            {
-                //label pointI = meshPoints[i];
-                //Pout<< "   " << pointI << " moved from "
-                //    << mesh_.points()[pointI] << " by "
-                //    << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
-                //    << endl;
-                newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
-            }
-        }
     }
 
-    // Per face (internal or coupled!) the patch that the
+    return facePatch;
-    // baffle should get (or -1).
+}
-    labelList facePatch(mesh_.nFaces(), -1);
-    // Count of baffled faces
-    label nBaffleFaces = 0;
 
 
-//    // Sync position? Or not since same face on both side so just sync
+//// Mark faces to be baffled to prevent snapping problems. Does
-//    // result of baffle.
+//// test to find nearest surface and checks which faces would get squashed.
+//Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
+//(
+//    const dictionary& motionDict,
+//    const labelList& globalToPatch
+//) const
+//{
+//    // Get the labels of added patches.
+//    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
 //
-//    const scalar minArea(readScalar(motionDict.lookup("minArea")));
+//    // Construct addressing engine.
+//    autoPtr<indirectPrimitivePatch> ppPtr
+//    (
+//        meshRefinement::makePatch
+//        (
+//            mesh_,
+//            adaptPatchIDs
+//        )
+//    );
+//    const indirectPrimitivePatch& pp = ppPtr();
+//    const pointField& localPoints = pp.localPoints();
+//    const labelList& meshPoints = pp.meshPoints();
 //
-//    Pout<< "markFacesOnProblemCellsGeometric : Comparing to minArea:"
+//    // Find nearest (non-baffle) surface
-//        << minArea << endl;
+//    pointField newPoints(mesh_.points());
-//
-//    pointField facePoints;
-//    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
 //    {
-//        const face& f = mesh_.faces()[faceI];
+//        List<pointIndexHit> hitInfo;
+//        labelList hitSurface;
+//        surfaces_.findNearest
+//        (
+//            surfaces_.getUnnamedSurfaces(),
+//            localPoints,
+//            scalarField(localPoints.size(), sqr(GREAT)),// sqr of attraction
+//            hitSurface,
+//            hitInfo
+//        );
 //
-//        bool usesPatchPoint = false;
+//        forAll(hitInfo, i)
-//
-//        facePoints.setSize(f.size());
-//        forAll(f, fp)
 //        {
-//            Map<label>::const_iterator iter = pp.meshPointMap().find(f[fp]);
+//            if (hitInfo[i].hit())
-//
-//            if (iter != pp.meshPointMap().end())
 //            {
-//                facePoints[fp] = newPosition[iter()];
+//                //label pointI = meshPoints[i];
-//                usesPatchPoint = true;
+//                //Pout<< "   " << pointI << " moved from "
+//                //    << mesh_.points()[pointI] << " by "
+//                //    << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
+//                //    << endl;
+//                newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
 //            }
-//            else
+//        }
+//    }
+//
+//    // Per face (internal or coupled!) the patch that the
+//    // baffle should get (or -1).
+//    labelList facePatch(mesh_.nFaces(), -1);
+//    // Count of baffled faces
+//    label nBaffleFaces = 0;
+//
+//    {
+//        pointField oldPoints(mesh_.points());
+//        mesh_.movePoints(newPoints);
+//        faceSet wrongFaces(mesh_, "wrongFaces", 100);
+//        {
+//            //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
+//
+//            // Just check the errors from squashing
+//            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+//
+//            const labelList allFaces(identity(mesh_.nFaces()));
+//            label nWrongFaces = 0;
+//
+//            scalar minArea(readScalar(motionDict.lookup("minArea")));
+//            if (minArea > -SMALL)
 //            {
-//                facePoints[fp] = mesh_.points()[f[fp]];
+//                polyMeshGeometry::checkFaceArea
+//                (
+//                    false,
+//                    minArea,
+//                    mesh_,
+//                    mesh_.faceAreas(),
+//                    allFaces,
+//                    &wrongFaces
+//                );
+//
+//                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")));
+//            scalar minDet = 0.01;
+//            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;
 //            }
 //        }
 //
-//        if (usesPatchPoint)
-//        {
-//            // Check area of face wrt original area
-//            face identFace(identity(f.size()));
 //
-//            if (identFace.mag(facePoints) < minArea)
+//        forAllConstIter(faceSet, wrongFaces, iter)
+//        {
+//            label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
+//
+//            if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
 //            {
-//                facePatch[faceI] = getBafflePatch(facePatch, faceI);
+//                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);
 //    }
 //
 //
-//    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+//    Info<< "markFacesOnProblemCellsGeometric : marked "
-//    forAll(patches, patchI)
+//        << returnReduce(nBaffleFaces, sumOp<label>())
-//    {
+//        << " additional internal and coupled faces"
-//        const polyPatch& pp = patches[patchI];
+//        << " to be converted into baffles." << endl;
 //
-//        if (pp.coupled())
+//    syncTools::syncFaceList
-//        {
+//    (
-//            forAll(pp, i)
+//        mesh_,
-//            {
+//        facePatch,
-//                label faceI = pp.start()+i;
+//        maxEqOp<label>(),
+//        false               // no separation
+//    );
 //
-//                const face& f = mesh_.faces()[faceI];
+//    return facePatch;
-//
+//}
-//                bool usesPatchPoint = false;
-//
-//                facePoints.setSize(f.size());
-//                forAll(f, fp)
-//                {
-//                    Map<label>::const_iterator iter =
-//                        pp.meshPointMap().find(f[fp]);
-//
-//                    if (iter != pp.meshPointMap().end())
-//                    {
-//                        facePoints[fp] = newPosition[iter()];
-//                        usesPatchPoint = true;
-//                    }
-//                    else
-//                    {
-//                        facePoints[fp] = mesh_.points()[f[fp]];
-//                    }
-//                }
-//
-//                if (usesPatchPoint)
-//                {
-//                    // Check area of face wrt original area
-//                    face identFace(identity(f.size()));
-//
-//                    if (identFace.mag(facePoints) < minArea)
-//                    {
-//                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
-//                        nBaffleFaces++;
-//                    }
-//                }
-//            }
-//        }
-//    }
-
-    {
-        pointField oldPoints(mesh_.points());
-        mesh_.movePoints(newPoints);
-        faceSet wrongFaces(mesh_, "wrongFaces", 100);
-        {
-            //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
-
-            // Just check the errors from squashing
-            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-            const labelList allFaces(identity(mesh_.nFaces()));
-            label nWrongFaces = 0;
-
-            scalar minArea(readScalar(motionDict.lookup("minArea")));
-            if (minArea > -SMALL)
-            {
-                polyMeshGeometry::checkFaceArea
-                (
-                    false,
-                    minArea,
-                    mesh_,
-                    mesh_.faceAreas(),
-                    allFaces,
-                    &wrongFaces
-                );
-
-                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")));
-            scalar minDet = 0.01;
-            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>(),
-        false               // no separation
-    );
-
-    return facePatch;
-}
-//XXXXXXXX
 
 
 // ************************************************************************* //