ENH: cvMesh: add filter function to remove nearly parallel edges

The edges are connected by a single point only.

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H

@@ -536,6 +536,11 @@ private:
             const Delaunay::Finite_edges_iterator& eit
         ) const;
 
+        boolList dualFaceBoundaryPoints
+        (
+            const Delaunay::Finite_edges_iterator& eit
+        ) const;
+
         //- Finds the maximum filterCount of the dual vertices
         //  (Delaunay cells) that form the dual face produced by the
         //  supplied edge
@@ -842,6 +847,13 @@ private:
             const Delaunay::Finite_facets_iterator& fit
         ) const;
 
+        //- Merge adjacent edges that are not attached to other faces
+        label mergeNearlyParallelEdges
+        (
+            const pointField& pts,
+            const scalar maxCosAngle
+        );
+
         //- Merge vertices that are very close together
         void mergeCloseDualVertices
         (

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C

@@ -212,6 +212,20 @@ void Foam::conformalVoronoiMesh::calcDualMesh
                     Info<< nl << "Collapsing unnecessary faces" << endl;
 
                     collapseFaces(points, boundaryPts, deferredCollapseFaces);
+
+                    const scalar maxCosAngle
+                        = cos(degToRad(cvMeshControls().edgeMergeAngle()));
+
+                    Info<< nl << "Merging adjacent edges which have an angle "
+                        << "of greater than "
+                        << cvMeshControls().edgeMergeAngle() << ": " << endl;
+
+                    label nRemovedEdges =
+                        mergeNearlyParallelEdges(points, maxCosAngle);
+
+                    reduce(nRemovedEdges, sumOp<label>());
+
+                    Info<< "    Merged " << nRemovedEdges << " edges" << endl;
                 }
 
                 labelHashSet wrongFaces = checkPolyMeshQuality(points);
@@ -652,6 +666,113 @@ Foam::label Foam::conformalVoronoiMesh::mergeCloseDualVertices
 }
 
 
+Foam::label Foam::conformalVoronoiMesh::mergeNearlyParallelEdges
+(
+    const pointField& pts,
+    const scalar maxCosAngle
+)
+{
+    List<HashSet<label> > pointFaceCount(number_of_cells());
+    labelList pointNeighbour(number_of_cells(), -1);
+
+    Map<label> dualPtIndexMap;
+
+    for
+    (
+        Delaunay::Finite_edges_iterator eit = finite_edges_begin();
+        eit != finite_edges_end();
+        ++eit
+    )
+    {
+        Cell_handle c = eit->first;
+        Vertex_handle vA = c->vertex(eit->second);
+        Vertex_handle vB = c->vertex(eit->third);
+
+        if (isBoundaryDualFace(eit))
+        {
+            const face f = buildDualFace(eit);
+
+            forAll(f, pI)
+            {
+                const label pIndex = f[pI];
+
+                const label prevPointI = f.prevLabel(pI);
+                const label nextPointI = f.nextLabel(pI);
+
+                pointFaceCount[pIndex].insert(prevPointI);
+                pointFaceCount[pIndex].insert(nextPointI);
+                pointNeighbour[pIndex] = nextPointI;
+            }
+        }
+        else if
+        (
+            vA->internalOrBoundaryPoint()
+         || vB->internalOrBoundaryPoint()
+        )
+        {
+            const face f = buildDualFace(eit);
+            const boolList faceBoundaryPoints = dualFaceBoundaryPoints(eit);
+
+            forAll(f, pI)
+            {
+                const label pIndex = f[pI];
+
+                const label prevPointI = f.prevLabel(pI);
+                const label nextPointI = f.nextLabel(pI);
+
+                pointFaceCount[pIndex].insert(prevPointI);
+                pointFaceCount[pIndex].insert(nextPointI);
+
+                if (faceBoundaryPoints[pI] == false)
+                {
+                    pointNeighbour[pIndex] = nextPointI;
+                }
+                else
+                {
+                    if (faceBoundaryPoints[prevPointI] == true)
+                    {
+                        pointNeighbour[pIndex] = prevPointI;
+                    }
+                    else if (faceBoundaryPoints[nextPointI] == true)
+                    {
+                        pointNeighbour[pIndex] = nextPointI;
+                    }
+                    else
+                    {
+                        pointNeighbour[pIndex] = pIndex;
+                    }
+                }
+            }
+        }
+    }
+
+    forAll(pointFaceCount, pI)
+    {
+        if (pointFaceCount[pI].size() == 2)
+        {
+            List<vector> edges(2, vector(0, 0, 0));
+
+            label count = 0;
+            forAllConstIter(HashSet<label>, pointFaceCount[pI], iter)
+            {
+                edges[count] = pts[pI] - pts[iter.key()];
+                edges[count] /= mag(edges[count]) + VSMALL;
+                count++;
+            }
+
+            if (mag(edges[0] & edges[1]) > maxCosAngle)
+            {
+                dualPtIndexMap.insert(pI, pointNeighbour[pI]);
+            }
+        }
+    }
+
+    reindexDualVertices(dualPtIndexMap);
+
+    return dualPtIndexMap.size();
+}
+
+
 void Foam::conformalVoronoiMesh::smoothSurface
 (
     pointField& pts,
@@ -696,7 +817,6 @@ void Foam::conformalVoronoiMesh::smoothSurface
     } while (nCollapsedFaces > 0);
 
     // Force all points of boundary faces to be on the surface
-
     for
     (
         Delaunay::Finite_cells_iterator cit = finite_cells_begin();
@@ -917,6 +1037,7 @@ void Foam::conformalVoronoiMesh::collapseFaces
 
         mergeCloseDualVertices(pts, boundaryPts);
 
+
         if (nCollapsedFaces > 0)
         {
             Info<< "    Collapsed " << nCollapsedFaces << " faces" << endl;
@@ -932,6 +1053,7 @@ void Foam::conformalVoronoiMesh::collapseFaces
         }
 
     } while (nCollapsedFaces > 0);
+
 }
 
 
@@ -1064,6 +1186,25 @@ Foam::conformalVoronoiMesh::collapseFace
 
     bool allowEarlyCollapseToPoint = true;
 
+
+//    // Quick exit
+//    label smallEdges = 0;
+//    const edgeList& fEdges = f.edges();
+//    forAll(fEdges, eI)
+//    {
+//        const edge& e = fEdges[eI];
+//
+//        if (e.mag(pts) < 0.2*targetFaceSize)
+//        {
+//            smallEdges++;
+//        }
+//    }
+//    if (smallEdges == 0)
+//    {
+//        return fcmNone;
+//    }
+
+
     // if (maxFC > cvMeshControls().filterCountSkipThreshold() - 3)
     // {
     //     limitToQuadsOrTris = true;
@@ -1071,6 +1212,7 @@ Foam::conformalVoronoiMesh::collapseFace
     //     allowEarlyCollapseToPoint = false;
     // }
 
+
     collapseSizeLimitCoeff *= pow
     (
         cvMeshControls().filterErrorReductionCoeff(),
@@ -1158,6 +1300,91 @@ Foam::conformalVoronoiMesh::collapseFace
         }
     }
 
+
+    scalar maxDist = 0;
+    scalar minDist = GREAT;
+//
+//    if (f.size() <= 3)
+//    {
+//        const edgeList& fEdges = f.edges();
+//
+//        forAll(fEdges, eI)
+//        {
+//            const edge& e = fEdges[eI];
+//            const scalar d = e.mag(pts);
+//
+//            if (d > maxDist)
+//            {
+//                maxDist = d;
+//                collapseAxis = e.vec(pts);
+//            }
+//            else if (d < minDist && d != 0)
+//            {
+//                minDist = d;
+//            }
+//        }
+//    }
+//    else
+//    {
+//        forAll(f, pI)
+//        {
+//            for (label i = pI + 1; i < f.size(); ++i)
+//            {
+//                if
+//                (
+//                    f[i] != f.nextLabel(pI)
+//                 && f[i] != f.prevLabel(pI)
+//                )
+//                {
+//                    scalar d = mag(pts[f[pI]] - pts[f[i]]);
+//
+//                    if (d > maxDist)
+//                    {
+//                        maxDist = d;
+//                        collapseAxis = pts[f[pI]] - pts[f[i]];
+//                    }
+//                    else if (d < minDist && d != 0)
+//                    {
+//                        minDist = d;
+//                    }
+//                }
+//            }
+//        }
+//    }
+//
+//    const edgeList& fEdges = f.edges();
+//
+//    scalar perimeter = 0;
+//
+//    forAll(fEdges, eI)
+//    {
+//        const edge& e = fEdges[eI];
+//        const scalar d = e.mag(pts);
+//
+//        perimeter += d;
+//
+////        collapseAxis += e.vec(pts);
+//
+//        if (d > maxDist)
+//        {
+//            collapseAxis = e.vec(pts);
+//            maxDist = d;
+//        }
+//        else if (d < minDist && d != 0)
+//        {
+//            minDist = d;
+//        }
+//    }
+//
+//    collapseAxis /= mag(collapseAxis);
+//
+////    Info<< f.size() << " " << minDist << " " << maxDist << " | "
+////        << collapseAxis << endl;
+//
+//    aspectRatio = maxDist/minDist;
+//
+//    aspectRatio = min(aspectRatio, sqr(perimeter)/(16.0*fA));
+
     if (magSqr(collapseAxis) < VSMALL)
     {
         WarningIn
@@ -1170,9 +1397,9 @@ Foam::conformalVoronoiMesh::collapseFace
         // Output face and collapse axis for visualisation
 
         Pout<< "# Aspect ratio = " << aspectRatio << nl
-            << "# inertia = " << J << nl
+//            << "# inertia = " << J << nl
-            << "# determinant = " << detJ << nl
+//            << "# determinant = " << detJ << nl
-            << "# eigenvalues = " << eigenValues(J) << nl
+//            << "# eigenvalues = " << eigenValues(J) << nl
             << "# collapseAxis = " << collapseAxis << nl
             << "# facePts = " << facePts << nl
             << endl;
@@ -1280,8 +1507,6 @@ Foam::conformalVoronoiMesh::collapseFace
     {
         scalar guardFraction = cvMeshControls().edgeCollapseGuardFraction();
 
-        cvMeshControls().maxCollapseFaceToPointSideLengthCoeff();
-
         if
         (
             allowEarlyCollapseToPoint
@@ -1337,6 +1562,8 @@ Foam::conformalVoronoiMesh::collapseFace
             Foam::point collapseToPt =
                 collapseAxis*(sum(dNeg)/dNeg.size() - dShift) + fC;
 
+//            DynamicList<label> faceBoundaryPts(f.size());
+
             forAll(facePtsNeg, fPtI)
             {
                 if (boundaryPts[facePtsNeg[fPtI]] == true)
@@ -1351,9 +1578,32 @@ Foam::conformalVoronoiMesh::collapseFace
                     collapseToPt = pts[collapseToPtI];
 
                     break;
+
+//                    faceBoundaryPts.append(facePtsNeg[fPtI]);
                 }
             }
 
+//            if (!faceBoundaryPts.empty())
+//            {
+//                if (faceBoundaryPts.size() == 2)
+//                {
+//                    collapseToPtI = faceBoundaryPts[0];
+//
+//                    collapseToPt =
+//                        0.5*(pts[faceBoundaryPts[0]] + pts[faceBoundaryPts[1]]);
+//                }
+//                else if (faceBoundaryPts.size() < f.size())
+//                {
+//                    face bFace(faceBoundaryPts);
+//
+//                    collapseToPtI = faceBoundaryPts.first();
+//
+//                    collapseToPt = bFace.centre(pts);
+//                }
+//            }
+//
+//            faceBoundaryPts.clear();
+
             // ...otherwise arbitrarily choosing the most distant
             // point as the index to collapse to.
 
@@ -1384,9 +1634,30 @@ Foam::conformalVoronoiMesh::collapseFace
                     collapseToPt = pts[collapseToPtI];
 
                     break;
+
+//                    faceBoundaryPts.append(facePtsNeg[fPtI]);
                 }
             }
 
+//            if (!faceBoundaryPts.empty())
+//            {
+//                if (faceBoundaryPts.size() == 2)
+//                {
+//                    collapseToPtI = faceBoundaryPts[0];
+//
+//                    collapseToPt =
+//                        0.5*(pts[faceBoundaryPts[0]] + pts[faceBoundaryPts[1]]);
+//                }
+//                else if (faceBoundaryPts.size() < f.size())
+//                {
+//                    face bFace(faceBoundaryPts);
+//
+//                    collapseToPtI = faceBoundaryPts.first();
+//
+//                    collapseToPt = bFace.centre(pts);
+//                }
+//            }
+
             // ...otherwise arbitrarily choosing the most distant
             // point as the index to collapse to.
 
@@ -1406,6 +1677,8 @@ Foam::conformalVoronoiMesh::collapseFace
 
             Foam::point collapseToPt = fC;
 
+            DynamicList<label> faceBoundaryPts(f.size());
+
             forAll(facePts, fPtI)
             {
                 if (boundaryPts[facePts[fPtI]] == true)
@@ -1415,11 +1688,32 @@ Foam::conformalVoronoiMesh::collapseFace
                     // use the first boundary point encountered if
                     // there are multiple boundary points.
 
-                    collapseToPtI = facePts[fPtI];
+//                    collapseToPtI = facePts[fPtI];
+//
+//                    collapseToPt = pts[collapseToPtI];
+//
+//                    break;
 
-                    collapseToPt = pts[collapseToPtI];
+                    faceBoundaryPts.append(facePts[fPtI]);
+                }
+            }
 
-                    break;
+            if (!faceBoundaryPts.empty())
+            {
+                if (faceBoundaryPts.size() == 2)
+                {
+                    collapseToPtI = faceBoundaryPts[0];
+
+                    collapseToPt =
+                        0.5*(pts[faceBoundaryPts[0]] + pts[faceBoundaryPts[1]]);
+                }
+                else if (faceBoundaryPts.size() < f.size())
+                {
+                    face bFace(faceBoundaryPts);
+
+                    collapseToPtI = faceBoundaryPts.first();
+
+                    collapseToPt = bFace.centre(pts);
                 }
             }
 
@@ -1858,7 +2152,10 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::findOffsetPatchFaces
         }
     }
 
-    offsetBoundaryCells.write();
+    if (cvMeshControls().objOutput())
+    {
+        offsetBoundaryCells.write();
+    }
 
     return offsetBoundaryCells;
 }
@@ -2092,13 +2389,7 @@ void Foam::conformalVoronoiMesh::indexDualVertices
 
             pts[dualVertI] = cit->dual();
 
-            if
+            if (cit->boundaryDualVertex())
-            (
-                !cit->vertex(0)->internalOrBoundaryPoint()
-             || !cit->vertex(1)->internalOrBoundaryPoint()
-             || !cit->vertex(2)->internalOrBoundaryPoint()
-             || !cit->vertex(3)->internalOrBoundaryPoint()
-            )
             {
                 // This is a boundary dual vertex
                 boundaryPts[dualVertI] = true;

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H

@@ -366,6 +366,49 @@ inline bool Foam::conformalVoronoiMesh::isBoundaryDualFace
 }
 
 
+inline Foam::List<bool> Foam::conformalVoronoiMesh::dualFaceBoundaryPoints
+(
+    const Delaunay::Finite_edges_iterator& eit
+) const
+{
+    Cell_circulator ccStart = incident_cells(*eit);
+    Cell_circulator cc1 = ccStart;
+    Cell_circulator cc2 = cc1;
+
+    // Advance the second circulator so that it always stays on the next
+    // cell around the edge;
+    cc2++;
+
+    DynamicList<bool> tmpFaceBoundaryPoints;
+
+    do
+    {
+        label cc1I = cc1->cellIndex();
+
+        label cc2I = cc2->cellIndex();
+
+        if (cc1I != cc2I)
+        {
+            if (cc1->boundaryDualVertex())
+            {
+                tmpFaceBoundaryPoints.append(true);
+            }
+            else
+            {
+                tmpFaceBoundaryPoints.append(false);
+            }
+        }
+
+        cc1++;
+
+        cc2++;
+
+    } while (cc1 != ccStart);
+
+    return tmpFaceBoundaryPoints;
+}
+
+
 inline Foam::List<Foam::label> Foam::conformalVoronoiMesh::processorsAttached
 (
     const Delaunay::Finite_facets_iterator& fit