ENH: Update collapesEdges to work with any supplied face zone

applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C

@@ -72,10 +72,11 @@ int main(int argc, char *argv[])
         "Collapse small and sliver faces as well as small edges"
     );
 
-    argList::addBoolOption
+    argList::addOption
     (
-        "collapseIndirectPatchFaces",
+        "collapseFaceZone",
-        "Collapse faces that are in the face zone indirectPatchFaces"
+        "zoneName",
+        "Collapse faces that are in the supplied face zone"
     );
 
 #   include "addOverwriteOption.H"
@@ -92,8 +93,7 @@ int main(int argc, char *argv[])
     const bool overwrite = args.optionFound("overwrite");
 
     const bool collapseFaces = args.optionFound("collapseFaces");
-    const bool collapseIndirectPatchFaces =
+    const bool collapseFaceZone = args.optionFound("collapseFaceZone");
-        args.optionFound("collapseIndirectPatchFaces");
 
     forAll(timeDirs, timeI)
     {
@@ -115,11 +115,15 @@ int main(int argc, char *argv[])
             meshMod.changeMesh(mesh, false);
         }
 
-        if (collapseIndirectPatchFaces)
+        if (collapseFaceZone)
         {
+            const word faceZoneName = args.optionRead<word>("collapseFaceZone");
+
+            const faceZone& fZone = mesh.faceZones()[faceZoneName];
+
             // Filter faces. Pass in the number of bad faces that are present
             // from the previous edge filtering to use as a stopping criterion.
-            meshFilter.filterIndirectPatchFaces();
+            meshFilter.filterFaceZone(fZone);
             {
                 polyTopoChange meshMod(newMesh);
 

src/dynamicMesh/polyMeshFilter/polyMeshFilter.C

@@ -32,6 +32,7 @@ License
 #include "polyTopoChange.H"
 #include "globalIndex.H"
 #include "PackedBoolList.H"
+#include "faceZone.H"
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
@@ -1134,83 +1135,330 @@ Foam::label Foam::polyMeshFilter::filterEdges
 }
 
 
-Foam::label Foam::polyMeshFilter::filterIndirectPatchFaces()
+Foam::label Foam::polyMeshFilter::filterFaceZone(const faceZone& fZone)
 {
-    newMeshPtr_ = copyMesh(mesh_);
+    const label nOriginalBadFaces = 0;
-    fvMesh& newMesh = newMeshPtr_();
 
-    label nIterations = 0;
+    label nBadFaces = labelMax;
-    label nBadFaces = 0;
+    label nOuterIterations = 0;
 
-    while (true)
+    minEdgeLen_.resize(mesh_.nEdges(), minLen_);
+    faceFilterFactor_.resize(mesh_.nFaces(), initialFaceLengthFactor_);
+
+    // Maintain the number of times a point has been part of a bad face
+    labelList pointErrorCount(mesh_.nPoints(), 0);
+
+    // Main loop
+    // ~~~~~~~~~
+    // It tries and do some collapses, checks the resulting mesh and
+    // 'freezes' some edges (by marking in minEdgeLen) and tries again.
+    // This will iterate ultimately to the situation where every edge is
+    // frozen and nothing gets collapsed.
+    while
+    (
+        nOuterIterations < maxIterations_
+     && nBadFaces > nOriginalBadFaces
+    )
     {
-        Info<< nl << indent << "Iteration = "
+        Info<< nl << "Outer Iteration = " << nOuterIterations++ << nl
-            << nIterations++ << nl << incrIndent << endl;
+            << endl;
 
-        // Per edge collapse status
+        printScalarFieldStats("Edge Filter Factor", minEdgeLen_);
-        PackedBoolList collapseEdge(newMesh.nEdges());
+        printScalarFieldStats("Face Filter Factor", faceFilterFactor_);
 
-        Map<point> collapsePointToLocation(newMesh.nPoints());
+        // Reset the new mesh to the old mesh
+        newMeshPtr_ = copyMesh(mesh_);
+        fvMesh& newMesh = newMeshPtr_();
 
-        labelList boundaryPoint(newMesh.nPoints());
+        scalarField newMeshFaceFilterFactor = faceFilterFactor_;
 
-        edgeCollapser collapser(newMesh, collapseFacesCoeffDict_);
+        labelList origToCurrentPointMap(identity(newMesh.nPoints()));
-
-        collapser.markIndirectPatchFaces
-        (
-            collapseEdge,
-            collapsePointToLocation
-        );
-
-        // Merge edge collapses into consistent collapse-network.
-        // Make sure no cells get collapsed.
-        List<pointEdgeCollapse> allPointInfo;
-        const globalIndex globalPoints(newMesh.nPoints());
-
-        collapser.consistentCollapse
-        (
-            globalPoints,
-            boundaryPoint,
-            collapsePointToLocation,
-            collapseEdge,
-            allPointInfo
-        );
-
-        label nLocalCollapse = collapseEdge.count();
-
-        reduce(nLocalCollapse, sumOp<label>());
-        Info<< nl << indent << "Collapsing " << nLocalCollapse
-            << " edges after synchronisation and PointEdgeWave" << endl;
-
-        if (nLocalCollapse == 0)
         {
-            break;
+
+            label nInnerIterations = 0;
+            label nPrevLocalCollapse = labelMax;
+
+            Info<< incrIndent;
+
+            while (true)
+            {
+                Info<< nl << indent << "Inner iteration = "
+                    << nInnerIterations++ << nl << incrIndent << endl;
+
+                // Per edge collapse status
+                PackedBoolList collapseEdge(newMesh.nEdges());
+
+                Map<point> collapsePointToLocation(newMesh.nPoints());
+
+                // Mark points on boundary
+                const labelList boundaryPoint = findBoundaryPoints(newMesh);
+
+                edgeCollapser collapser(newMesh, collapseFacesCoeffDict_);
+
+                {
+                    // Collapse faces
+                    labelPair nCollapsedPtEdge = collapser.markFaceZoneEdges
+                    (
+                        fZone,
+                        newMeshFaceFilterFactor,
+                        boundaryPoint,
+                        collapseEdge,
+                        collapsePointToLocation
+                    );
+
+                    label nCollapsed = 0;
+                    forAll(nCollapsedPtEdge, collapseTypeI)
+                    {
+                        nCollapsed += nCollapsedPtEdge[collapseTypeI];
+                    }
+
+                    reduce(nCollapsed, sumOp<label>());
+
+                    Info<< indent
+                        << "Collapsing " << nCollapsed << " faces"
+                        << " (to point = "
+                        << returnReduce
+                           (
+                               nCollapsedPtEdge.first(),
+                               sumOp<label>()
+                           )
+                        << ", to edge = "
+                        << returnReduce
+                           (
+                               nCollapsedPtEdge.second(),
+                               sumOp<label>()
+                           )
+                        << ")" << endl;
+
+                    if (nCollapsed == 0)
+                    {
+                        Info<< decrIndent;
+                        Info<< decrIndent;
+                        break;
+                    }
+                }
+
+                // Merge edge collapses into consistent collapse-network.
+                // Make sure no cells get collapsed.
+                List<pointEdgeCollapse> allPointInfo;
+                const globalIndex globalPoints(newMesh.nPoints());
+
+                collapser.consistentCollapse
+                (
+                    globalPoints,
+                    boundaryPoint,
+                    collapsePointToLocation,
+                    collapseEdge,
+                    allPointInfo
+                );
+
+                label nLocalCollapse = collapseEdge.count();
+
+                reduce(nLocalCollapse, sumOp<label>());
+                Info<< nl << indent << "Collapsing " << nLocalCollapse
+                    << " edges after synchronisation and PointEdgeWave" << endl;
+
+                if (nLocalCollapse >= nPrevLocalCollapse)
+                {
+                    Info<< decrIndent;
+                    Info<< decrIndent;
+                    break;
+                }
+                else
+                {
+                    nPrevLocalCollapse = nLocalCollapse;
+                }
+
+                {
+                    // Apply collapses to current mesh
+                    polyTopoChange newMeshMod(newMesh);
+
+                    // Insert mesh refinement into polyTopoChange.
+                    collapser.setRefinement(allPointInfo, newMeshMod);
+
+                    Info<< indent << "Apply changes to the current mesh"
+                        << decrIndent << endl;
+
+                    // Apply changes to current mesh
+                    autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
+                    (
+                        newMesh,
+                        false
+                    );
+                    const mapPolyMesh& newMap = newMapPtr();
+
+                    // Update fields
+                    newMesh.updateMesh(newMap);
+                    if (newMap.hasMotionPoints())
+                    {
+                        newMesh.movePoints(newMap.preMotionPoints());
+                    }
+
+                    // Relabel the boundary points
+        //            labelList newBoundaryPoints(newMesh.nPoints(), -1);
+        //            forAll(newBoundaryPoints, pI)
+        //            {
+        //                const label newToOldptI= map.pointMap()[pI];
+        //                newBoundaryPoints[pI] = boundaryIOPts[newToOldptI];
+        //            }
+        //            boundaryIOPts = newBoundaryPoints;
+
+
+                    mapOldMeshFaceFieldToNewMesh
+                    (
+                        newMesh,
+                        newMap.faceMap(),
+                        newMeshFaceFilterFactor
+                    );
+
+                    updateOldToNewPointMap
+                    (
+                        newMap.reversePointMap(),
+                        origToCurrentPointMap
+                    );
+                }
+            }
         }
 
-        // Apply collapses to current mesh
-        polyTopoChange newMeshMod(newMesh);
 
-        // Insert mesh refinement into polyTopoChange.
+        scalarField newMeshMinEdgeLen = minEdgeLen_;
-        collapser.setRefinement(allPointInfo, newMeshMod);
 
-        Info<< indent << "Apply changes to the current mesh"
+        label nInnerIterations = 0;
-            << decrIndent << endl;
+        label nPrevLocalCollapse = labelMax;
 
-        // Apply changes to current mesh
+        while (true)
-        autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
-        (
-            newMesh,
-            false
-        );
-        const mapPolyMesh& newMap = newMapPtr();
-
-        // Update fields
-        newMesh.updateMesh(newMap);
-        if (newMap.hasMotionPoints())
         {
-            newMesh.movePoints(newMap.preMotionPoints());
+            Info<< nl << indent << "Inner iteration = "
+                << nInnerIterations++ << nl << incrIndent << endl;
+
+            // Per edge collapse status
+            PackedBoolList collapseEdge(newMesh.nEdges());
+
+            Map<point> collapsePointToLocation(newMesh.nPoints());
+
+            // Mark points on boundary
+            const labelList boundaryPoint = findBoundaryPoints
+            (
+                newMesh//,
+//                    boundaryIOPts
+            );
+
+            edgeCollapser collapser(newMesh, collapseFacesCoeffDict_);
+
+            // Work out which edges to collapse
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            // This is by looking at minEdgeLen (to avoid frozen edges)
+            // and marking in collapseEdge.
+            label nSmallCollapsed = collapser.markSmallEdges
+            (
+                newMeshMinEdgeLen,
+                boundaryPoint,
+                collapseEdge,
+                collapsePointToLocation
+            );
+
+            reduce(nSmallCollapsed, sumOp<label>());
+            Info<< indent << "Collapsing " << nSmallCollapsed
+                << " small edges" << endl;
+
+            // Merge inline edges
+            label nMerged = collapser.markMergeEdges
+            (
+                maxCos_,
+                boundaryPoint,
+                collapseEdge,
+                collapsePointToLocation
+            );
+
+            reduce(nMerged, sumOp<label>());
+            Info<< indent << "Collapsing " << nMerged << " in line edges"
+                << endl;
+
+            if (nMerged + nSmallCollapsed == 0)
+            {
+                Info<< decrIndent;
+                break;
+            }
+
+            // Merge edge collapses into consistent collapse-network.
+            // Make sure no cells get collapsed.
+            List<pointEdgeCollapse> allPointInfo;
+            const globalIndex globalPoints(newMesh.nPoints());
+
+            collapser.consistentCollapse
+            (
+                globalPoints,
+                boundaryPoint,
+                collapsePointToLocation,
+                collapseEdge,
+                allPointInfo
+            );
+
+            label nLocalCollapse = collapseEdge.count();
+
+            reduce(nLocalCollapse, sumOp<label>());
+            Info<< nl << indent << "Collapsing " << nLocalCollapse
+                << " edges after synchronisation and PointEdgeWave" << endl;
+
+            if (nLocalCollapse >= nPrevLocalCollapse)
+            {
+                Info<< decrIndent;
+                break;
+            }
+            else
+            {
+                nPrevLocalCollapse = nLocalCollapse;
+            }
+
+            // Apply collapses to current mesh
+            polyTopoChange newMeshMod(newMesh);
+
+            // Insert mesh refinement into polyTopoChange.
+            collapser.setRefinement(allPointInfo, newMeshMod);
+
+            Info<< indent << "Apply changes to the current mesh"
+                << decrIndent << endl;
+
+            // Apply changes to current mesh
+            autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
+            (
+                newMesh,
+                false
+            );
+            const mapPolyMesh& newMap = newMapPtr();
+
+            // Update fields
+            newMesh.updateMesh(newMap);
+            if (newMap.hasMotionPoints())
+            {
+                newMesh.movePoints(newMap.preMotionPoints());
+            }
+
+            // Relabel the boundary points
+//            labelList newBoundaryPoints(newMesh.nPoints(), -1);
+//            forAll(newBoundaryPoints, pI)
+//            {
+//                const label newToOldptI= map.pointMap()[pI];
+//                newBoundaryPoints[pI] = boundaryIOPts[newToOldptI];
+//            }
+//            boundaryIOPts = newBoundaryPoints;
+
+            // Synchronise the factors
+            mapOldMeshEdgeFieldToNewMesh
+            (
+                newMesh,
+                newMap.pointMap(),
+                newMeshMinEdgeLen
+            );
+
+            updateOldToNewPointMap
+            (
+                newMap.reversePointMap(),
+                origToCurrentPointMap
+            );
         }
 
+
         // Mesh check
         // ~~~~~~~~~~~~~~~~~~
         // Do not allow collapses in regions of error.
@@ -1230,6 +1478,33 @@ Foam::label Foam::polyMeshFilter::filterIndirectPatchFaces()
                 << "    Number of marked points : "
                 << returnReduce(isErrorPoint.count(), sumOp<unsigned int>())
                 << endl;
+
+            updatePointErrorCount
+            (
+                isErrorPoint,
+                origToCurrentPointMap,
+                pointErrorCount
+            );
+
+            checkMeshEdgesAndRelaxEdges
+            (
+                newMesh,
+                origToCurrentPointMap,
+                isErrorPoint,
+                pointErrorCount
+            );
+
+            checkMeshFacesAndRelaxEdges
+            (
+                newMesh,
+                origToCurrentPointMap,
+                isErrorPoint,
+                pointErrorCount
+            );
+        }
+        else
+        {
+            return -1;
         }
     }
 

src/dynamicMesh/polyMeshFilter/polyMeshFilter.H

@@ -52,6 +52,7 @@ namespace Foam
 class polyMesh;
 class fvMesh;
 class PackedBoolList;
+class faceZone;
 
 /*---------------------------------------------------------------------------*\
                        Class polyMeshFilter Declaration
@@ -238,7 +239,7 @@ public:
             label filterEdges(const label nOriginalBadFaces);
 
             //- Filter all faces that are in the face zone indirectPatchFaces
-            label filterIndirectPatchFaces();
+            label filterFaceZone(const faceZone& fZone);
 };
 
 

src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C

@@ -32,6 +32,7 @@ License
 #include "globalIndex.H"
 #include "removePoints.H"
 #include "motionSmoother.H"
+#include "OFstream.H"
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
@@ -2029,94 +2030,151 @@ Foam::labelPair Foam::edgeCollapser::markSmallSliverFaces
 }
 
 
-void Foam::edgeCollapser::markIndirectPatchFaces
+Foam::labelPair Foam::edgeCollapser::markFaceZoneEdges
 (
+    const faceZone& fZone,
+    const scalarField& faceFilterFactor,
+    const labelList& pointPriority,
     PackedBoolList& collapseEdge,
     Map<point>& collapsePointToLocation
 ) const
 {
-    const faceZone& indirectFaceZone = mesh_.faceZones()["indirectPatchFaces"];
+    const faceList& faces = mesh_.faces();
 
-    const edgeList& edges = mesh_.edges();
+    const scalarField targetFaceSizes = calcTargetFaceSizes();
-    const pointField& points = mesh_.points();
-    const labelListList& edgeFaces = mesh_.edgeFaces();
-    const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
 
-    forAll(edges, eI)
+    // Calculate number of faces that will be collapsed to a point or an edge
+    label nCollapseToPoint = 0;
+    label nCollapseToEdge = 0;
+
+    forAll(faces, fI)
     {
-        const edge& e = edges[eI];
+        if (fZone.whichFace(fI) == -1)
-
-        const labelList& eFaces = edgeFaces[eI];
-
-        bool keepEdge = false;
-
-        label nInternalFaces = 0;
-        label nPatchFaces = 0;
-        label nIndirectFaces = 0;
-
-        bool coupled = false;
-
-        forAll(eFaces, eFaceI)
         {
-            const label eFaceIndex = eFaces[eFaceI];
+            continue;
-
-            if (mesh_.isInternalFace(eFaceIndex))
-            {
-                nInternalFaces++;
-            }
-            else
-            {
-                const label patchIndex = bMesh.whichPatch(eFaceIndex);
-                const polyPatch& pPatch = bMesh[patchIndex];
-
-                if (pPatch.coupled())
-                {
-                    coupled = true;
-                    nInternalFaces++;
-                }
-                else
-                {
-                    // Keep the edge if an attached face is not in the face zone
-                    if (indirectFaceZone.whichFace(eFaceIndex) == -1)
-                    {
-                        nPatchFaces++;
-                    }
-                    else
-                    {
-                        nIndirectFaces++;
-                    }
-                }
-            }
         }
 
-        if (eFaces.size() != nInternalFaces + nPatchFaces + nIndirectFaces)
+        const face& f = faces[fI];
+
+        if (faceFilterFactor[fI] == 0)
         {
-            Pout<< eFaces.size() << " ("
+            continue;
-                << nInternalFaces << "/" << nPatchFaces << "/" << nIndirectFaces
-                << ")" << endl;
         }
 
-        if
+        collapseType flagCollapseFace = collapseFace
         (
-            eFaces.size() == nInternalFaces
+            pointPriority,
-         || nIndirectFaces < (coupled ? 1 : 2)
+            f,
-        )
+            fI,
+            targetFaceSizes[fI],
+            collapseEdge,
+            collapsePointToLocation,
+            faceFilterFactor
+        );
+
+        if (flagCollapseFace == noCollapse)
         {
-            keepEdge = true;
+            continue;
         }
-
+        else if (flagCollapseFace == toPoint)
-        if (!keepEdge)
         {
-            collapseEdge[eI] = true;
+            nCollapseToPoint++;
-
+        }
-            const Foam::point collapsePoint =
+        else if (flagCollapseFace == toEdge)
-                0.5*(points[e.end()] + points[e.start()]);
+        {
-
+            nCollapseToEdge++;
-            collapsePointToLocation.insert(e.start(), collapsePoint);
+        }
-            collapsePointToLocation.insert(e.end(), collapsePoint);
+        else
+        {
+            FatalErrorIn("collapseFaces(const polyMesh&, List<labelPair>&)")
+                << "Face is marked to be collapsed to " << flagCollapseFace
+                << ". Currently can only collapse to point/edge."
+                << abort(FatalError);
         }
     }
 
+    return labelPair(nCollapseToPoint, nCollapseToEdge);
+
+//    const edgeList& edges = mesh_.edges();
+//    const pointField& points = mesh_.points();
+//    const labelListList& edgeFaces = mesh_.edgeFaces();
+//    const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
+//
+//    forAll(edges, eI)
+//    {
+//        const edge& e = edges[eI];
+//
+//        const labelList& eFaces = edgeFaces[eI];
+//
+//        bool keepEdge = false;
+//
+//        label nInternalFaces = 0;
+//        label nPatchFaces = 0;
+//        label nIndirectFaces = 0;
+//
+//        bool coupled = false;
+//
+//        forAll(eFaces, eFaceI)
+//        {
+//            const label eFaceIndex = eFaces[eFaceI];
+//
+//            if (mesh_.isInternalFace(eFaceIndex))
+//            {
+//                nInternalFaces++;
+//            }
+//            else
+//            {
+//                const label patchIndex = bMesh.whichPatch(eFaceIndex);
+//                const polyPatch& pPatch = bMesh[patchIndex];
+//
+//                if (pPatch.coupled())
+//                {
+//                    coupled = true;
+//                    nInternalFaces++;
+//                }
+//                else
+//                {
+//                    // Keep the edge if an attached face is not in the zone
+//                    if (fZone.whichFace(eFaceIndex) == -1)
+//                    {
+//                        nPatchFaces++;
+//                    }
+//                    else
+//                    {
+//                        nIndirectFaces++;
+//                    }
+//                }
+//            }
+//        }
+//
+//        if (eFaces.size() != nInternalFaces + nPatchFaces + nIndirectFaces)
+//        {
+//            Pout<< eFaces.size() << " ("
+//                << nInternalFaces << "/" << nPatchFaces << "/"
+//                << nIndirectFaces << ")" << endl;
+//        }
+//
+//        if
+//        (
+//            eFaces.size() == nInternalFaces
+//         || nIndirectFaces < (coupled ? 1 : 2)
+//        )
+//        {
+//            keepEdge = true;
+//        }
+//
+//        if (!keepEdge)
+//        {
+//            collapseEdge[eI] = true;
+//
+//            const Foam::point collapsePoint =
+//                0.5*(points[e.end()] + points[e.start()]);
+//
+//            collapsePointToLocation.insert(e.start(), collapsePoint);
+//            collapsePointToLocation.insert(e.end(), collapsePoint);
+//        }
+//    }
+
 //    OFstream str
 //    (
 //        mesh_.time().path()

src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H

@@ -339,8 +339,11 @@ public:
             ) const;
 
             //- Marks edges in the faceZone indirectPatchFaces for collapse
-            void markIndirectPatchFaces
+            labelPair markFaceZoneEdges
             (
+                const faceZone& fZone,
+                const scalarField& faceFilterFactor,
+                const labelList& pointPriority,
                 PackedBoolList& collapseEdge,
                 Map<point>& collapsePointToLocation
             ) const;