ENH: conformalMeshCalcDualMesh: protect merging processor points

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

@@ -1294,11 +1294,12 @@ Foam::labelList Foam::backgroundMeshDecomposition::processorPosition
                     "Foam::labelList"
                     "Foam::backgroundMeshDecomposition::processorPosition"
                     "("
-                        "const List<point>& pts"
+                        "const List<point>&"
                     ") const"
                 )
                     << "The position " << pts[pI]
-                    << " is not in any part of the background mesh.  "
+                    << " is not in any part of the background mesh "
+                    << globalBackgroundBounds_ << endl
                     << "A background mesh with a wider margin around "
                     << "the geometry may help."
                     << exit(FatalError);
@@ -1311,11 +1312,11 @@ Foam::labelList Foam::backgroundMeshDecomposition::processorPosition
                     "Foam::labelList"
                     "Foam::backgroundMeshDecomposition::processorPosition"
                     "("
-                        "const List<point>& pts"
+                        "const List<point>&"
                     ") const"
-                )
+                )   << "The position " << pts[pI]
-                    << "The position " << pts[pI]
+                    << " was not found in the background mesh "
-                    << " was not found in the background mesh, finding nearest."
+                    << globalBackgroundBounds_ << ", finding nearest."
                     << endl;
             }
 

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

@@ -30,6 +30,16 @@ License
 #include "backgroundMeshDecomposition.H"
 #include "meshSearch.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(conformalVoronoiMesh, 0);
+
+}
+
+
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
 Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment

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

@@ -722,6 +722,13 @@ private:
             const Delaunay::Finite_facets_iterator& fit
         ) const;
 
+        //- Determines if the edge constructed from the face is on
+        //  a processor patch
+        inline bool isParallelDualEdge
+        (
+            const Delaunay::Finite_facets_iterator& fit
+        ) const;
+
         //- Merge vertices that are very close together
         void mergeCloseDualVertices
         (
@@ -771,7 +778,7 @@ private:
         ) const;
 
         //- Collapse a face to an edge, updating the point and point
-        // map.  Returns the collapse mode that was applied.
+        //  map. Returns the collapse mode that was applied.
         faceCollapseMode collapseFace
         (
             const face& f,
@@ -783,7 +790,7 @@ private:
             label maxFC
         ) const;
 
-        // Identify the index of the longest edge on the face
+        //- Identify the index of the longest edge on the face
         label longestEdge(const face& f, const pointField& pts) const;
 
         //- Identify the face labels of the deferred collapse faces
@@ -886,6 +893,10 @@ private:
 
 public:
 
+    //- Runtime type information
+    ClassName("conformalVoronoiMesh");
+
+
     // Constructors
 
         //- Construct from Time and cvMeshDict

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

@@ -543,8 +543,21 @@ Foam::label Foam::conformalVoronoiMesh::mergeCloseDualVertices
 {
     label nPtsMerged = 0;
 
+    label nIdentical = 0;
+    label nProcEdge = 0;
+
+    // Relative distance for points to be merged
     scalar closenessTolerance = cvMeshControls().mergeClosenessCoeff();
 
+    // Absolute distance for points to be considered coincident. Bit adhoc
+    // but points were seen with distSqr ~ 1E-30 which is SMALL^2. Add a few
+    // digits to account for truncation errors.
+    scalar coincidentDistanceSqr = sqr
+    (
+        SMALL*1E2*geometryToConformTo_.globalBounds().mag()
+    );
+
+
     for
     (
         Delaunay::Finite_facets_iterator fit = finite_facets_begin();
@@ -567,14 +580,14 @@ Foam::label Foam::conformalVoronoiMesh::mergeCloseDualVertices
             continue;
         }
 
-        if (!c1->farCell() && !c2->farCell() && (c1I != c2I))
+        if ((c1I != c2I) && !c1->farCell() && !c2->farCell())
         {
-            if
+            scalar distSqr = magSqr(pts[c1I] - pts[c2I]);
-            (
+
-                magSqr(pts[c1I] - pts[c2I])
+            if (pts[c1I] == pts[c2I] || distSqr < coincidentDistanceSqr)
-              < sqr(averageAnyCellSize(fit)*closenessTolerance)
-            )
             {
+                nIdentical++;
+
                 if (boundaryPts[c2I] == true)
                 {
                     // If c2I is a boundary point, then it is kept.
@@ -583,18 +596,59 @@ Foam::label Foam::conformalVoronoiMesh::mergeCloseDualVertices
 
                     dualPtIndexMap.insert(c1I, c2I);
                     dualPtIndexMap.insert(c2I, c2I);
+                    nPtsMerged++;
                 }
                 else
                 {
                     dualPtIndexMap.insert(c1I, c1I);
                     dualPtIndexMap.insert(c2I, c1I);
+                    nPtsMerged++;
                 }
 
-                nPtsMerged++;
+            }
+            else if (distSqr < sqr(averageAnyCellSize(fit)*closenessTolerance))
+            {
+                if (c1->parallelDualVertex() || c2->parallelDualVertex())
+                //if (isParallelDualEdge(fit))
+                {
+                    // Skip if face uses any edge that becomes a processor
+                    // dual face.
+                    // Note: the real test should be whether the Delaunay edge
+                    //  will form a processor patch.
+                    nProcEdge++;
+                }
+                else if (boundaryPts[c2I] == true)
+                {
+                    // If c2I is a boundary point, then it is kept.
+                    // If both are boundary points then c2I is chosen
+                    // arbitrarily to be kept.
+
+                    dualPtIndexMap.insert(c1I, c2I);
+                    dualPtIndexMap.insert(c2I, c2I);
+                    nPtsMerged++;
+                }
+                else
+                {
+                    dualPtIndexMap.insert(c1I, c1I);
+                    dualPtIndexMap.insert(c2I, c1I);
+                    nPtsMerged++;
+                }
             }
         }
     }
 
+    if (debug)
+    {
+        Info<< "mergeCloseDualVertices:"
+            << " coincident distance:" << coincidentDistanceSqr
+            << " closenessTolerance:" << closenessTolerance << endl
+            << "    identical points : "
+            << returnReduce(nIdentical, sumOp<label>()) << endl
+            << "    processor edges  : "
+            << returnReduce(nProcEdge, sumOp<label>()) << endl
+            << endl;
+    }
+
     return nPtsMerged;
 }
 

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

@@ -384,6 +384,23 @@ inline Foam::List<Foam::label> Foam::conformalVoronoiMesh::processorsAttached
 }
 
 
+inline bool Foam::conformalVoronoiMesh::isParallelDualEdge
+(
+    const Delaunay::Finite_facets_iterator& fit
+) const
+{
+    const Cell_handle c1(fit->first);
+    const int oppositeVertex = fit->second;
+
+    return
+    (
+            c1->vertex(vertex_triple_index(oppositeVertex, 0))->referred()
+         || c1->vertex(vertex_triple_index(oppositeVertex, 1))->referred()
+         || c1->vertex(vertex_triple_index(oppositeVertex, 2))->referred()
+    );
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 #ifdef CGAL_INEXACT

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

@@ -208,7 +208,6 @@ public:
         );
     }
 
-
     //- Does the Dual vertex form part of a processor patch
     inline Foam::label dualVertexMasterProc() const
     {