ENH: cvMesh: speed up clipLineToProc.

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

@@ -178,8 +178,7 @@ Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment
             "("
                 "const Foam::point& pt"
             ") const"
-        )
+        )   << "No secondary surface hit found in spoke search "
-            << "No secondary surface hit found in spoke search "
             << "using " << s
             << " spokes, try increasing alignmentSearchSpokes."
             << endl;

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

@@ -124,8 +124,15 @@ public:
 
 private:
 
+    // Static data
+
         static const scalar tolParallel;
 
+        static const scalar searchConeAngle;
+
+        static const scalar searchAngleOppositeSurface;
+
+
     // Private data
 
         //- The time registry of the application
@@ -420,6 +427,8 @@ private:
         //- Determine and insert point groups at the feature points
         void insertFeaturePoints();
 
+        bool edgesShareNormal(const label e1, const label e2) const;
+
         //- Create point groups at convex feature points
         void createConvexFeaturePoints
         (
@@ -470,6 +479,11 @@ private:
         //- Store the locations of all of the features to be conformed to
         void constructFeaturePointLocations();
 
+        List<pointIndexHit> findSurfacePtLocationsNearFeaturePoint
+        (
+            const Foam::point& featurePoint
+        ) const;
+
         //- Reinsert stored feature point defining points
         void reinsertFeaturePoints(bool distribute = false);
 
@@ -582,6 +596,7 @@ private:
         //  intersect.
         bool clipLineToProc
         (
+            const Foam::point& pt,
             Foam::point& a,
             Foam::point& b
         ) const;
@@ -699,9 +714,9 @@ private:
         ) const;
 
         //- Check if a point is near any feature edge points.
-        bool pointIsNearFeatureEdge(const Foam::point& pt) const;
+        bool pointIsNearFeatureEdgeLocation(const Foam::point& pt) const;
 
-        bool pointIsNearFeatureEdge
+        bool pointIsNearFeatureEdgeLocation
         (
             const Foam::point& pt,
             pointIndexHit& info

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

@@ -29,6 +29,12 @@ License
 
 using namespace Foam::vectorTools;
 
+const Foam::scalar Foam::conformalVoronoiMesh::searchConeAngle
+    = Foam::cos(degToRad(30));
+
+const Foam::scalar Foam::conformalVoronoiMesh::searchAngleOppositeSurface
+    = Foam::cos(degToRad(150));
+
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
 void Foam::conformalVoronoiMesh::conformToSurface()
@@ -277,79 +283,6 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
         timeCheck("After initial conformation");
 
         initialTotalHits = nSurfHits + nFeatEdHits;
-
-        // Filter small edges at the boundary by inserting surface point pairs
-//        for
-//        (
-//            Delaunay::Finite_cells_iterator cit = finite_cells_begin();
-//            cit != finite_cells_end();
-//            cit++
-//        )
-//        {
-//            const Foam::point dualVertex = topoint(dual(cit));
-//
-//
-//        }
-//        for
-//        (
-//            Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-//            vit != finite_vertices_end();
-//            vit++
-//        )
-//        {
-//            if (vit->nearBoundary())
-//            {
-//                std::list<Facet> faces;
-//                incident_facets(vit, std::back_inserter(faces));
-//
-//                List<scalar> edgeLengthList(faces.size());
-//                scalar totalLength = 0;
-//                label count = 0;
-//
-//                for
-//                (
-//                    std::list<Facet>::iterator fit=faces.begin();
-//                    fit != faces.end();
-//                    ++fit
-//                )
-//                {
-//                    if
-//                    (
-//                        is_infinite(fit->first)
-//                     || is_infinite(fit->first->neighbor(fit->second))
-//                    )
-//                    {
-//                        continue;
-//                    }
-//
-//                    const Point dE0 = dual(fit->first);
-//                    const Point dE1 = dual(fit->first->neighbor(fit->second));
-//
-//                    const Segment s(dE0, dE1);
-//
-//                    const scalar length = Foam::sqrt(s.squared_length());
-//
-//                    edgeLengthList[count++] = length;
-//
-//                    totalLength += length;
-//
-//                    //Info<< length << " / " << totalLength << endl;
-//                }
-//
-//                const scalar averageLength = totalLength/edgeLengthList.size();
-//
-//                forAll(edgeLengthList, eI)
-//                {
-//                    const scalar& el = edgeLengthList[eI];
-//
-//                    if (el < 0.1*averageLength)
-//                    {
-//                        //Info<< "SMALL" << endl;
-//                    }
-//                }
-//            }
-//        }
-
     }
 
     // Remember which vertices were referred to each processor so only updates
@@ -598,7 +531,7 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection
             Foam::point& a = dE0;
             Foam::point& b = dE1;
 
-            bool inProc = clipLineToProc(a, b);
+            bool inProc = clipLineToProc(topoint(vit->point()), a, b);
 
             // Check for the edge passing through a surface
             if
@@ -642,7 +575,7 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
 
     for
     (
-        std::list<Facet>::iterator fit=facets.begin();
+        std::list<Facet>::iterator fit = facets.begin();
         fit != facets.end();
         ++fit
     )
@@ -666,7 +599,7 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
 
         if (Pstream::parRun())
         {
-            bool inProc = clipLineToProc(dE0, dE1);
+            bool inProc = clipLineToProc(topoint(vit->point()), dE0, dE1);
 
             if (!inProc)
             {
@@ -684,8 +617,6 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
 
         if (infoIntersection.hit())
         {
-            flagIntersection = true;
-
             vectorField norm(1);
 
             allGeometry_[hitSurfaceIntersection].getNormal
@@ -709,10 +640,10 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
             pointIndexHit info;
             label hitSurface = -1;
 
-            geometryToConformTo_.findSurfaceNearestIntersection
+            geometryToConformTo_.findSurfaceNearest
             (
-                vertex,
+                newPoint,
-                newPoint + SMALL*n,
+                surfaceSearchDistanceSqr(newPoint),
                 info,
                 hitSurface
             );
@@ -738,12 +669,12 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
                             = mag(p - info.hitPoint());
 
                         const scalar minSepDist
-                            = cvMeshControls().removalDistCoeff()
+                            = sqr(cvMeshControls().removalDistCoeff()
-                             *targetCellSize(p);
+                             *targetCellSize(p));
 
                         // Reject the point if it is too close to another
                         // surface point.
-                        // Could merge?
+                        // Could merge the points?
                         if (separationDistance < minSepDist)
                         {
                             rejectPoint = true;
@@ -757,6 +688,7 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
             // because another surface may be in the way.
             if (!rejectPoint && info.hit())
             {
+                flagIntersection = true;
                 infoList.append(info);
                 hitSurfaceList.append(hitSurface);
             }
@@ -769,72 +701,45 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
 
 bool Foam::conformalVoronoiMesh::clipLineToProc
 (
+    const Foam::point& pt,
     Foam::point& a,
     Foam::point& b
 ) const
 {
-    bool intersects = false;
+    bool inProc = false;
 
-    if (decomposition_().positionOnThisProcessor(a))
+    pointIndexHit findAnyIntersection = decomposition_().findLine(a, b);
+
+    if (!findAnyIntersection.hit())
     {
-        // a is inside this processor
+        pointIndexHit info = decomposition_().findLine(a, pt);
 
-        if (decomposition_().positionOnThisProcessor(b))
+        if (!info.hit())
         {
-            // both a and b inside, no clip required
+            inProc = true;
-            intersects = true;
         }
         else
         {
-            // b is outside, clip b to bounding box.
+            inProc = false;
-
-            pointIndexHit info = decomposition_().findLine(b, a);
-
-            if (info.hit())
-            {
-                intersects = true;
-                b = info.hitPoint();
-            }
         }
     }
     else
     {
-        // a is outside this processor
+        pointIndexHit info = decomposition_().findLine(a, pt);
 
-        if (decomposition_().positionOnThisProcessor(b))
+        if (!info.hit())
         {
-            // b is inside this processor, clip a to processor
+            inProc = true;
-
+            b = findAnyIntersection.hitPoint();
-            pointIndexHit info = decomposition_().findLine(a, b);
-
-            if (info.hit())
-            {
-                intersects = true;
-                a = info.hitPoint();
-            }
         }
         else
         {
-            // both a and b outside, but they can still intersect the processor
+            inProc = true;
-
+            a = findAnyIntersection.hitPoint();
-            pointIndexHit info = decomposition_().findLine(a, b);
-
-            if (info.hit())
-            {
-                intersects = true;
-                a = info.hitPoint();
-
-                info = decomposition_().findLine(b, a);
-
-                if (info.hit())
-                {
-                    b = info.hitPoint();
-                }
-            }
         }
     }
 
-    return intersects;
+    return inProc;
 }
 
 
@@ -1877,7 +1782,7 @@ Foam::scalar Foam::conformalVoronoiMesh::angleBetweenSurfacePoints
 
     const vector nB = norm[0];
 
-    return vectorTools::degAngleBetween(nA, nB);
+    return vectorTools::cosPhi(nA, nB);
 }
 
 
@@ -1896,11 +1801,11 @@ bool Foam::conformalVoronoiMesh::nearSurfacePoint
 
     if (closeToSurfacePt)
     {
-        const scalar angle
+        const scalar cosAngle
             = angleBetweenSurfacePoints(pt, closePoint.hitPoint());
 
         // @todo make this tolerance run-time selectable?
-        if (angle > 170.0)
+        if (cosAngle < searchAngleOppositeSurface)
         {
             pointIndexHit pCloseHit;
             label pCloseSurfaceHit = -1;
@@ -1999,7 +1904,7 @@ bool Foam::conformalVoronoiMesh::appendToEdgeLocationTree
 }
 
 
-bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdge
+bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdgeLocation
 (
     const Foam::point& pt
 ) const
@@ -2013,7 +1918,7 @@ bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdge
 }
 
 
-bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdge
+bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdgeLocation
 (
     const Foam::point& pt,
     pointIndexHit& info
@@ -2066,7 +1971,7 @@ bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
 
     pointIndexHit info;
 
-    bool closeToFeatureEdge = pointIsNearFeatureEdge(pt, info);
+    bool closeToFeatureEdge = pointIsNearFeatureEdgeLocation(pt, info);
 
     if (!closeToFeatureEdge)
     {
@@ -2075,7 +1980,7 @@ bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
     else
     {
         // Check if the edge location that the new edge location is near to
-        // might be on a different edge. If so, add it anyway.
+        // "might" be on a different edge. If so, add it anyway.
         pointIndexHit edgeHit;
         label featureHit = -1;
 
@@ -2094,14 +1999,17 @@ bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
 
         const vector lineBetweenPoints = pt - info.hitPoint();
 
-        const scalar angle = degAngleBetween(edgeDir, lineBetweenPoints);
+        const scalar cosAngle = vectorTools::cosPhi(edgeDir, lineBetweenPoints);
 
         // Allow the point to be added if it is almost at right angles to the
         // other point. Also check it is not the same point.
+//        Info<< cosAngle<< " "
+//            << radToDeg(acos(cosAngle)) << " "
+//            << searchConeAngle << " "
+//            << radToDeg(acos(searchConeAngle)) << endl;
         if
         (
-            angle < 100.0
+            mag(cosAngle) < searchConeAngle
-         && angle > 80.0
          && mag(lineBetweenPoints) > SMALL
         )
         {
@@ -2274,7 +2182,16 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
             continue;
         }
 
-        if (nearFeaturePt(surfHitI.hitPoint()))
+        bool isNearFeaturePt = nearFeaturePt(surfHitI.hitPoint());
+
+        bool isNearSurfacePt = nearSurfacePoint
+        (
+            surfHitI,
+            hitSurfaceI,
+            existingSurfacePtLocations
+        );
+
+        if (isNearFeaturePt || isNearSurfacePt)
         {
             keepSurfacePoint = false;
 
@@ -2289,19 +2206,6 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
             }
         }
 
-        if
-        (
-            nearSurfacePoint
-            (
-                surfHitI,
-                hitSurfaceI,
-                existingSurfacePtLocations
-            )
-        )
-        {
-            keepSurfacePoint = false;
-        }
-
         List<List<pointIndexHit> > edHitsByFeature;
 
         labelList featuresHit;
@@ -2332,7 +2236,7 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
             {
                 pointIndexHit& edHit = edHits[eHitI];
 
-                if (!nearFeaturePt(edHit.hitPoint()))
+                if (!nearFeaturePt(edHit.hitPoint()) && keepSurfacePoint)
                 {
                     if
                     (
@@ -2345,6 +2249,8 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
                         // this will prevent "pits" forming.
 
                         keepSurfacePoint = false;
+
+                        // NEED TO REMOVE FROM THE SURFACE TREE...
                     }
 
                     if

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

@@ -84,12 +84,12 @@ void Foam::conformalVoronoiMesh::createEdgePointGroup
         }
         case extendedFeatureEdgeMesh::OPEN:
         {
-            createOpenEdgePointGroup(feMesh, edHit, pts, indices, types);
+            //createOpenEdgePointGroup(feMesh, edHit, pts, indices, types);
             break;
         }
         case extendedFeatureEdgeMesh::MULTIPLE:
         {
-            createMultipleEdgePointGroup(feMesh, edHit, pts, indices, types);
+            //createMultipleEdgePointGroup(feMesh, edHit, pts, indices, types);
             break;
         }
         case extendedFeatureEdgeMesh::NONE:
@@ -124,7 +124,6 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
     {
         // The normals are nearly parallel, so this is too sharp a feature to
         // conform to.
-
         return;
     }
 
@@ -448,6 +447,40 @@ void Foam::conformalVoronoiMesh::reinsertFeaturePoints(bool distribute)
 }
 
 
+bool Foam::conformalVoronoiMesh::edgesShareNormal
+(
+    const label e1,
+    const label e2
+) const
+{
+    const PtrList<extendedFeatureEdgeMesh>& feMeshes
+    (
+        geometryToConformTo_.features()
+    );
+
+    forAll(feMeshes, i)
+    {
+        const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
+
+        const vectorField& e1normals = feMesh.edgeNormals(e1);
+        const vectorField& e2normals = feMesh.edgeNormals(e2);
+
+        forAll(e1normals, nI1)
+        {
+            forAll(e2normals, nI2)
+            {
+                if (degAngleBetween(e1normals[nI1], e2normals[nI2]) < 1)
+                {
+                    return true;
+                }
+            }
+        }
+    }
+
+    return false;
+}
+
+
 void Foam::conformalVoronoiMesh::createConvexFeaturePoints
 (
     DynamicList<Foam::point>& pts,
@@ -479,6 +512,7 @@ void Foam::conformalVoronoiMesh::createConvexFeaturePoints
             }
 
             const vectorField& featPtNormals = feMesh.featurePointNormals(ptI);
+
             const scalar ppDist = - pointPairDistance(featPt);
 
             vector cornerNormal = sum(featPtNormals);
@@ -706,6 +740,61 @@ void Foam::conformalVoronoiMesh::createMixedFeaturePoints
     }
 }
 
+//
+//void Foam::conformalVoronoiMesh::createFeaturePoints
+//(
+//    DynamicList<Foam::point>& pts,
+//    DynamicList<label>& indices,
+//    DynamicList<label>& types
+//)
+//{
+//    const PtrList<extendedFeatureEdgeMesh>& feMeshes
+//    (
+//        geometryToConformTo_.features()
+//    );
+//
+//    forAll(feMeshes, i)
+//    {
+//        const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
+//
+//        for
+//        (
+//            label ptI = feMesh.convexStart();
+//            ptI < feMesh.mixedStart();
+//            ++ptI
+//        )
+//        {
+//            const Foam::point& featPt = feMesh.points()[ptI];
+//
+//            if (!positionOnThisProc(featPt))
+//            {
+//                continue;
+//            }
+//
+//            const scalar searchRadiusSqr = 5.0*targetCellSize(featPt);
+//
+//            labelList indices =
+//                surfacePtLocationTreePtr_().findSphere(featPt, searchRadiusSqr);
+//
+//            pointField nearestSurfacePoints(indices.size());
+//
+//            forAll(indices, pI)
+//            {
+//                nearestSurfacePoints[pI] =
+//                    surfaceConformationVertices_[indices[pI]];
+//            }
+//
+//            forAll()
+//
+//            // Now find the nearest points within the edge cones.
+//
+//            // Calculate preliminary surface point locations
+//
+//
+//        }
+//    }
+//}
+
 
 void Foam::conformalVoronoiMesh::insertFeaturePoints()
 {
@@ -721,6 +810,8 @@ void Foam::conformalVoronoiMesh::insertFeaturePoints()
 
     createConcaveFeaturePoints(pts, indices, types);
 
+//    createFeaturePoints(pts, indices, types);
+
     createMixedFeaturePoints(pts, indices, types);
 
     // Insert the created points, distributing to the appropriate processor
@@ -798,3 +889,35 @@ void Foam::conformalVoronoiMesh::constructFeaturePointLocations()
 
     featurePointLocations_.transfer(ftPtLocs);
 }
+
+
+Foam::List<Foam::pointIndexHit>
+Foam::conformalVoronoiMesh::findSurfacePtLocationsNearFeaturePoint
+(
+    const Foam::point& featurePoint
+) const
+{
+    DynamicList<pointIndexHit> dynPointList;
+
+    const scalar searchRadiusSqr = 3*targetCellSize(featurePoint);
+
+    labelList surfacePtList = surfacePtLocationTreePtr_().findSphere
+    (
+       featurePoint,
+       searchRadiusSqr
+    );
+
+    forAll(surfacePtList, elemI)
+    {
+       label index = surfacePtList[elemI];
+
+       const Foam::point& p
+           = surfacePtLocationTreePtr_().shapes().shapePoints()[index];
+
+       pointIndexHit nearHit(true, p, index);
+
+       dynPointList.append(nearHit);
+    }
+
+    return dynPointList.shrink();
+}