ENH: Improve overlap testing.

Only test circumspheres that overlap the current processor boundary.

Sphere overlap test is just a findNearest test, no need for findSphere.

Function for cell: internalOrBoundaryDualVertex.

src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C

@@ -1159,6 +1159,16 @@ bool Foam::backgroundMeshDecomposition::overlapsThisProcessor
 }
 
 
+bool Foam::backgroundMeshDecomposition::overlapsThisProcessor
+(
+    const point& centre,
+    scalar radiusSqr
+) const
+{
+    return bFTreePtr_().findNearest(centre, radiusSqr).hit();
+}
+
+
 Foam::pointIndexHit Foam::backgroundMeshDecomposition::findLine
 (
     const point& start,
@@ -1169,6 +1179,16 @@ Foam::pointIndexHit Foam::backgroundMeshDecomposition::findLine
 }
 
 
+Foam::pointIndexHit Foam::backgroundMeshDecomposition::findLineAny
+(
+    const point& start,
+    const point& end
+) const
+{
+    return bFTreePtr_().findLineAny(start, end);
+}
+
+
 Foam::labelList Foam::backgroundMeshDecomposition::processorPosition
 (
     const List<point>& pts
@@ -1601,8 +1621,8 @@ Foam::labelListList Foam::backgroundMeshDecomposition::overlapsProcessors
         const point& c = testCentres[sI];
         scalar rSqr = testRadiusSqrs[sI];
 
-        // If the sphere finds some elements of the patch, then it overlaps
-        sphereOverlapsCandidate[sI] = !bFTreePtr_().findSphere(c, rSqr).empty();
+        // If the sphere finds a nearest element of the patch, then it overlaps
+        sphereOverlapsCandidate[sI] = bFTreePtr_().findNearest(c, rSqr).hit();
     }
 
     map().reverseDistribute

src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H

@@ -229,6 +229,14 @@ public:
         //  processor
         bool overlapsThisProcessor(const treeBoundBox& box) const;
 
+        //- Does the given sphere overlap the faces of the bounday of this
+        //  processor
+        bool overlapsThisProcessor
+        (
+            const point& centre,
+            scalar radiusSqr
+        ) const;
+
         //- Find nearest intersection of line between start and end, (exposing
         //  underlying indexedOctree)
         pointIndexHit findLine
@@ -237,6 +245,14 @@ public:
             const point& end
         ) const;
 
+        //- Find any intersection of line between start and end, (exposing
+        //  underlying indexedOctree)
+        pointIndexHit findLineAny
+        (
+            const point& start,
+            const point& end
+        ) const;
+
         //- What processor is the given position on?
         labelList processorPosition(const List<point>& pts) const;
 

src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C

@@ -1964,13 +1964,7 @@ void Foam::conformalVoronoiMesh::move()
     {
         cit->cellIndex() = -1;
 
-        if
-        (
-            cit->vertex(0)->anyInternalOrBoundaryPoint()
-         || cit->vertex(1)->anyInternalOrBoundaryPoint()
-         || cit->vertex(2)->anyInternalOrBoundaryPoint()
-         || cit->vertex(3)->anyInternalOrBoundaryPoint()
-        )
+        if (cit->anyInternalOrBoundaryDualVertex())
         {
             cit->cellIndex() = dualVertI;
 
@@ -2336,7 +2330,7 @@ bool Foam::conformalVoronoiMesh::positionOnThisProc
 
 Foam::boolList Foam::conformalVoronoiMesh::positionOnThisProc
 (
-    const Foam::pointField& pts
+    const Foam::List<Foam::point>& pts
 ) const
 {
     if (Pstream::parRun())
@@ -2350,7 +2344,7 @@ Foam::boolList Foam::conformalVoronoiMesh::positionOnThisProc
 
 Foam::labelList Foam::conformalVoronoiMesh::positionProc
 (
-    const Foam::pointField& pts
+    const Foam::List<Foam::point>& pts
 ) const
 {
     if (!Pstream::parRun())
@@ -2365,8 +2359,8 @@ Foam::labelList Foam::conformalVoronoiMesh::positionProc
 Foam::List<Foam::List<Foam::pointIndexHit> >
 Foam::conformalVoronoiMesh::intersectsProc
 (
-    const pointField& starts,
-    const pointField& ends
+    const List<Foam::point>& starts,
+    const List<Foam::point>& ends
 ) const
 {
     if (!Pstream::parRun())
@@ -2380,8 +2374,8 @@ Foam::conformalVoronoiMesh::intersectsProc
 
 Foam::labelListList Foam::conformalVoronoiMesh::overlapsProc
 (
-    const pointField& centres,
-    const scalarField& radiusSqrs
+    const List<Foam::point>& centres,
+    const List<scalar>& radiusSqrs
 ) const
 {
     if (!Pstream::parRun())

src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H

@@ -702,13 +702,6 @@ private:
             const Delaunay::Finite_edges_iterator& eit
         ) const;
 
-        //- Determines if the dual vertex constructed by the Delaunay
-        //  cell is a boundary vertex
-        inline bool isBoundaryDualVertex
-        (
-            const Delaunay::Finite_cells_iterator& cit
-        ) const;
-
         //- Merge vertices that are very close together
         void mergeCloseDualVertices
         (
@@ -895,23 +888,23 @@ public:
         bool positionOnThisProc(const Foam::point& pt) const;
 
         //- Check if the point is in the domain handled by this processor
-        boolList positionOnThisProc(const Foam::pointField& pts) const;
+        boolList positionOnThisProc(const Foam::List<Foam::point>& pts) const;
 
         //- Which processor's domain handles this point
-        labelList positionProc(const Foam::pointField& pts) const;
+        labelList positionProc(const Foam::List<Foam::point>& pts) const;
 
         //- Which other processors does each line segment intersect
         List<List<pointIndexHit> > intersectsProc
         (
-            const pointField& starts,
-            const pointField& ends
+            const List<Foam::point>& starts,
+            const List<Foam::point>& ends
         ) const;
 
         //- Which other processors does each sphere overlap
         labelListList overlapsProc
         (
-            const pointField& centres,
-            const scalarField& radiusSqrs
+            const List<Foam::point>& centres,
+            const List<scalar>& radiusSqrs
         ) const;
 
         //- Conversion functions between point (FOAM::) and Point (CGAL)

src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C

@@ -254,13 +254,7 @@ void Foam::conformalVoronoiMesh::calcDualMesh
 //         ++cit
 //     )
 //     {
-//         if
-//         (
-//             cit->vertex(0)->internalOrBoundaryPoint()
-//          || cit->vertex(1)->internalOrBoundaryPoint()
-//          || cit->vertex(2)->internalOrBoundaryPoint()
-//          || cit->vertex(3)->internalOrBoundaryPoint()
-//         )
+//         if (cit->internalOrBoundaryDualVertex())
 //         {
 //              cit->cellIndex() = cellI++;
 //         }
@@ -420,6 +414,7 @@ void Foam::conformalVoronoiMesh::mergeCloseDualVertices
     // Assess close points to be merged
 
     label nPtsMerged = 0;
+    label nPtsMergedSum = 0;
 
     do
     {
@@ -427,14 +422,20 @@ void Foam::conformalVoronoiMesh::mergeCloseDualVertices
 
         nPtsMerged = mergeCloseDualVertices(pts, boundaryPts, dualPtIndexMap);
 
-        if (nPtsMerged > 0)
-        {
-            Pout<< "    Merged " << nPtsMerged << " points " << endl;
-        }
+        // if (nPtsMerged > 0)
+        // {
+        //     Pout<< "    Merged " << nPtsMerged << " points " << endl;
+        // }
+
+        nPtsMergedSum += nPtsMerged;
 
         reindexDualVertices(dualPtIndexMap);
 
     } while (nPtsMerged > 0);
+
+    Info<< "    Merged "
+        << returnReduce(nPtsMergedSum, sumOp<label>())
+        << " points " << endl;
 }
 
 
@@ -1733,13 +1734,7 @@ void Foam::conformalVoronoiMesh::indexDualVertices
         ++cit
     )
     {
-        if
-        (
-            cit->vertex(0)->internalOrBoundaryPoint()
-         || cit->vertex(1)->internalOrBoundaryPoint()
-         || cit->vertex(2)->internalOrBoundaryPoint()
-         || cit->vertex(3)->internalOrBoundaryPoint()
-        )
+        if (cit->internalOrBoundaryDualVertex())
         {
             cit->cellIndex() = dualVertI;
 

src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C

@@ -740,9 +740,7 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceIntersection
     )
     {
         const Cell_handle c1(fit->first);
-        // oppositeVertex number
-        const int oV = fit->second;
-        const Cell_handle c2(c1->neighbor(oV));
+        const Cell_handle c2(c1->neighbor(fit->second));
 
         // If either Delaunay cell at the end of the Dual edge is infinite,
         // skip.
@@ -755,18 +753,8 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceIntersection
         // i.e. all vertices form part of the internal or boundary definition
         if
         (
-            (
-            c1->vertex(oV)->internalOrBoundaryPoint()
-         || c1->vertex(vertex_triple_index(oV, 0))->internalOrBoundaryPoint()
-         || c1->vertex(vertex_triple_index(oV, 1))->internalOrBoundaryPoint()
-         || c1->vertex(vertex_triple_index(oV, 2))->internalOrBoundaryPoint()
-            )
-         && (
-            c2->vertex(oV)->internalOrBoundaryPoint()
-         || c2->vertex(vertex_triple_index(oV, 0))->internalOrBoundaryPoint()
-         || c2->vertex(vertex_triple_index(oV, 1))->internalOrBoundaryPoint()
-         || c2->vertex(vertex_triple_index(oV, 2))->internalOrBoundaryPoint()
-            )
+            c1->internalOrBoundaryDualVertex()
+         && c2->internalOrBoundaryDualVertex()
         )
         {
             dE0[fI] = topoint(dual(c1));
@@ -799,9 +787,7 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceIntersection
     )
     {
         const Cell_handle c1(fit->first);
-        // oppositeVertex number
-        const int oV = fit->second;
-        const Cell_handle c2(c1->neighbor(oV));
+        const Cell_handle c2(c1->neighbor(fit->second));
 
         // If either Delaunay cell at the end of the Dual edge is infinite,
         // skip.
@@ -814,18 +800,8 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceIntersection
         // i.e. all vertices form part of the internal or boundary definition
         if
         (
-            (
-            c1->vertex(oV)->internalOrBoundaryPoint()
-         || c1->vertex(vertex_triple_index(oV, 0))->internalOrBoundaryPoint()
-         || c1->vertex(vertex_triple_index(oV, 1))->internalOrBoundaryPoint()
-         || c1->vertex(vertex_triple_index(oV, 2))->internalOrBoundaryPoint()
-            )
-         && (
-            c2->vertex(oV)->internalOrBoundaryPoint()
-         || c2->vertex(vertex_triple_index(oV, 0))->internalOrBoundaryPoint()
-         || c2->vertex(vertex_triple_index(oV, 1))->internalOrBoundaryPoint()
-         || c2->vertex(vertex_triple_index(oV, 2))->internalOrBoundaryPoint()
-            )
+            c1->internalOrBoundaryDualVertex()
+         && c2->internalOrBoundaryDualVertex()
         )
         {
             const Foam::point a = dE0[fI];
@@ -957,8 +933,6 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceIntersection
         }
     }
 
-    timeCheck("buildParallelInterfaceIntersection before referVertices");
-
     referVertices
     (
         targetProcessor,
@@ -983,10 +957,15 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
     DynamicList<label> parallelInfluenceIndices;
 
     // Some of these values will not be used, i.e. for non-real cells
-    pointField circumcentre(number_of_cells(), vector::one*GREAT);
-    scalarField circumradiusSqr(number_of_cells(), scalar(0.0));
+    DynamicList<Foam::point> circumcentre;
+    DynamicList<scalar> circumradiusSqr;
+
+    PackedBoolList testCellInfluence(number_of_cells(), false);
 
-    label cI = 0;
+    // Index outer (all) Delaunauy cells for whether they are potential
+    // overlaps, index (inner) the list of tests an results.
+    label cIInner = 0;
+    label cIOuter = 0;
 
     for
     (
@@ -1002,27 +981,28 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
 
         // The Delaunay cells to assess have to be real, i.e. all vertices form
         // part of the internal or boundary definition
-        if
-        (
-               cit->vertex(0)->internalOrBoundaryPoint()
-            || cit->vertex(1)->internalOrBoundaryPoint()
-            || cit->vertex(2)->internalOrBoundaryPoint()
-            || cit->vertex(3)->internalOrBoundaryPoint()
-        )
+        if (cit->internalOrBoundaryDualVertex())
         {
-            circumcentre[cI] = topoint(dual(cit));
+            Foam::point cc(topoint(dual(cit)));
 
-            circumradiusSqr[cI] = magSqr
+            scalar crSqr
             (
-                circumcentre[cI] - topoint(cit->vertex(0)->point())
+                magSqr(cc - topoint(cit->vertex(0)->point()))
             );
 
-            cI++;
+            // Only if the circumsphere overlaps the boundary of this processor
+            // is there a chance of it overlapping others
+            if (decomposition_().overlapsThisProcessor(cc, crSqr))
+            {
+                circumcentre.append(cc);
+                circumradiusSqr.append(crSqr);
+
+                testCellInfluence[cIOuter] = true;
+            }
         }
-    }
 
-    circumcentre.setSize(cI);
-    circumradiusSqr.setSize(cI);
+        cIOuter++;
+    }
 
     // Increasing the circumspheres to increase the overlaps and compensate for
     // floating point errors missing some referrals
@@ -1031,7 +1011,9 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
         overlapsProc(circumcentre, sqr(1.01)*circumradiusSqr)
     );
 
-    cI = 0;
+    // Reset counters
+    cIInner = 0;
+    cIOuter = 0;
 
     // Relying on the order of iteration of cells being the same as before
     for
@@ -1041,17 +1023,10 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
         ++cit
     )
     {
-        // The Delaunay cells to assess have to be real, i.e. all vertices form
-        // part of the internal or boundary definition
-        if
-        (
-            cit->vertex(0)->internalOrBoundaryPoint()
-         || cit->vertex(1)->internalOrBoundaryPoint()
-         || cit->vertex(2)->internalOrBoundaryPoint()
-         || cit->vertex(3)->internalOrBoundaryPoint()
-        )
+        // Pre-tested circumsphere potential influence
+        if (testCellInfluence[cIOuter])
         {
-            const labelList& citOverlaps = circumsphereOverlaps[cI];
+            const labelList& citOverlaps = circumsphereOverlaps[cIInner];
 
             forAll(citOverlaps, cOI)
             {
@@ -1090,8 +1065,10 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
                 }
             }
 
-            cI++;
+            cIInner++;
         }
+
+        cIOuter++;
     }
 
     referVertices
@@ -1116,6 +1093,8 @@ void Foam::conformalVoronoiMesh::referVertices
     const word& outputName
 )
 {
+    timeCheck("Start of referVertices " + stageName);
+
     if (cvMeshControls().objOutput())
     {
         writePoints
@@ -1183,6 +1162,8 @@ void Foam::conformalVoronoiMesh::referVertices
     }
 
     Info<< "    Total " << stageName << " vertices " << totalVertices << endl;
+
+    timeCheck("End of referVertices " + stageName);
 }
 
 

src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H

@@ -357,32 +357,6 @@ inline bool Foam::conformalVoronoiMesh::isBoundaryDualFace
 }
 
 
-inline bool Foam::conformalVoronoiMesh::isBoundaryDualVertex
-(
-    const Delaunay::Finite_cells_iterator& cit
-) const
-{
-    // A dual vertex on the boundary will result from a Delaunay cell
-    // with least one Delaunay vertex outside and at least one inside
-
-    return
-    (
-        (
-            cit->vertex(0)->internalOrBoundaryPoint()
-         || cit->vertex(1)->internalOrBoundaryPoint()
-         || cit->vertex(2)->internalOrBoundaryPoint()
-         || cit->vertex(3)->internalOrBoundaryPoint()
-        )
-     && (
-            !cit->vertex(0)->internalOrBoundaryPoint()
-         || !cit->vertex(1)->internalOrBoundaryPoint()
-         || !cit->vertex(2)->internalOrBoundaryPoint()
-         || !cit->vertex(3)->internalOrBoundaryPoint()
-        )
-    );
-}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 #ifdef CGAL_INEXACT

src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H

@@ -137,6 +137,56 @@ public:
     {
         return filterCount_;
     }
+
+
+    // Is the Delaunay cell real, i.e. all vertices form part of the internal
+    // or boundary definition
+    inline int internalOrBoundaryDualVertex() const
+    {
+        return
+        (
+            this->vertex(0)->internalOrBoundaryPoint()
+         || this->vertex(1)->internalOrBoundaryPoint()
+         || this->vertex(2)->internalOrBoundaryPoint()
+         || this->vertex(3)->internalOrBoundaryPoint()
+        );
+    }
+
+
+    // Is the Delaunay cell real or referred (or mixed), i.e. all vertices form
+    // part of the real or referred internal or boundary definition
+    inline int anyInternalOrBoundaryDualVertex() const
+    {
+        return
+        (
+            this->vertex(0)->anyInternalOrBoundaryPoint()
+         || this->vertex(1)->anyInternalOrBoundaryPoint()
+         || this->vertex(2)->anyInternalOrBoundaryPoint()
+         || this->vertex(3)->anyInternalOrBoundaryPoint()
+        );
+    }
+
+
+    // A dual vertex on the boundary will result from a Delaunay cell with
+    // least one Delaunay vertex outside and at least one inside
+    inline bool boundaryDualVertex() const
+    {
+        return
+        (
+            (
+               this->vertex(0)->internalOrBoundaryPoint()
+            || this->vertex(1)->internalOrBoundaryPoint()
+            || this->vertex(2)->internalOrBoundaryPoint()
+            || this->vertex(3)->internalOrBoundaryPoint()
+            )
+         && (
+               !this->vertex(0)->internalOrBoundaryPoint()
+            || !this->vertex(1)->internalOrBoundaryPoint()
+            || !this->vertex(2)->internalOrBoundaryPoint()
+            || !this->vertex(3)->internalOrBoundaryPoint()
+            )
+        );
+    }
 };