ENH: foamyHexMesh: Allow surface point insertion on close surfaces

ENH: foamyHexMesh: do not point motion to cross surfaces
ENH: foamyHexMesh: do not allow perturbation to cross surface
ENH: foamyHexMesh: change inside/outside test to check for nearest intersection
BUG: foamyHexMesh: With surfaces of same priority do not overwrite cell size
ENH: foamyHexMesh: insert surface points on close parallel surfaces
BUG: foamyHexMesh: write to own processor file in parallel
ENH: foamyHexMesh: Insert surface points if internal points are exposed
ENH: foamyHexMesh: add hasBoundaryPoint() to indexedCell
REVERT: foamyHexMesh: do not check inside/outside for surface points
REVERT: foamyHexMesh: correct function argument list
REVERT: foamyHexMesh: correct function argument list

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/controlMeshRefinement/controlMeshRefinement.C

@@ -295,14 +295,14 @@ void Foam::controlMeshRefinement::initialMeshPopulation
                     shapeController_.minimumCellSize()
                 );
             }
-            else if (maxPriority == controlFunction.maxPriority())
+//            else if (maxPriority == controlFunction.maxPriority())
-            {
+//            {
-                vertices[vI].targetCellSize() = max
+//                vertices[vI].targetCellSize() = max
-                (
+//                (
-                    min(sizes[vI], size),
+//                    min(sizes[vI], size),
-                    shapeController_.minimumCellSize()
+//                    shapeController_.minimumCellSize()
-                );
+//                );
-            }
+//            }
             else
             {
                 vertices[vI].targetCellSize() = max

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

@@ -1487,8 +1487,13 @@ void Foam::conformalVoronoiMesh::move()
     // Save displacements to file.
     if (foamyHexMeshControls().objOutput() && time().outputTime())
     {
-        Pout<< "Writing point displacement vectors to file." << endl;
+        Info<< "Writing point displacement vectors to file." << endl;
-        OFstream str("displacements_" + runTime_.timeName() + ".obj");
+        OFstream str
+        (
+            time().path()
+          + "displacements_" + runTime_.timeName()
+          + ".obj"
+        );
 
         for
         (

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

@@ -550,53 +550,130 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
             }
         }
 
+//        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
+//            (
+//                vA->referred()
+//             || vB->referred()
+//            )
+//            {
+//                continue;
+//            }
+//
+//            if
+//            (
+//                (vA->internalPoint() && vA->referred())
+//             || (vB->internalPoint() && vB->referred())
+//            )
+//            {
+//                continue;
+//            }
+//
+//            if
+//            (
+//                (vA->internalPoint() && vB->externalBoundaryPoint())
+//             || (vB->internalPoint() && vA->externalBoundaryPoint())
+//             || (vA->internalBoundaryPoint() && vB->internalBoundaryPoint())
+//            )
+//            {
+//                pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
+//                pointIndexHit surfHit;
+//                label hitSurface;
+//
+//                geometryToConformTo_.findSurfaceNearest
+//                (
+//                    (
+//                        vA->internalPoint()
+//                      ? topoint(vA->point())
+//                      : topoint(vB->point())
+//                    ),
+//                    magSqr(topoint(vA->point()) - topoint(vB->point())),
+//                    surfHit,
+//                    hitSurface
+//                );
+//
+//                if (surfHit.hit())
+//                {
+//                    surfaceIntersections.append
+//                    (
+//                        pointIndexHitAndFeature(surfHit, hitSurface)
+//                    );
+//
+//                    addSurfaceAndEdgeHits
+//                    (
+//                        (
+//                            vA->internalPoint()
+//                          ? topoint(vA->point())
+//                          : topoint(vB->point())
+//                        ),
+//                        surfaceIntersections,
+//                        surfacePtReplaceDistCoeffSqr,
+//                        edgeSearchDistCoeffSqr,
+//                        surfaceHits,
+//                        featureEdgeHits,
+//                        surfaceToTreeShape,
+//                        edgeToTreeShape,
+//                        surfacePtToEdgePtDist,
+//                        false
+//                    );
+//                }
+//            }
+//        }
+
         for
         (
-            Delaunay::Finite_edges_iterator eit = finite_edges_begin();
+            Delaunay::Finite_cells_iterator cit = finite_cells_begin();
-            eit != finite_edges_end();
+            cit != finite_cells_end();
-            ++eit
+            ++cit
         )
         {
-            Cell_handle c = eit->first;
+            label nInternal = 0;
-            Vertex_handle vA = c->vertex(eit->second);
+            for (label vI = 0; vI < 4; vI++)
-            Vertex_handle vB = c->vertex(eit->third);
-
-            if
-            (
-                vA->referred()
-             || vB->referred()
-            )
             {
-                continue;
+                if (cit->vertex(vI)->internalPoint())
+                {
+                    nInternal++;
+                }
             }
 
-            if
+            if (nInternal == 1 && cit->hasBoundaryPoint())
-            (
+            //if (cit->boundaryDualVertex() && !cit->hasReferredPoint())
-                (vA->internalPoint() && vA->referred())
-             || (vB->internalPoint() && vB->referred())
-            )
             {
-                continue;
+                const Foam::point& pt = cit->dual();
-            }
+
+                const Foam::point cellCentre =
+                    topoint
+                    (
+                        CGAL::centroid
+                        (
+                            CGAL::Tetrahedron_3<baseK>
+                            (
+                                cit->vertex(0)->point(),
+                                cit->vertex(1)->point(),
+                                cit->vertex(2)->point(),
+                                cit->vertex(3)->point()
+                            )
+                        )
+                    );
 
-            if
-            (
-                (vA->internalPoint() && vB->externalBoundaryPoint())
-             || (vB->internalPoint() && vA->externalBoundaryPoint())
-            )
-            {
                 pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
                 pointIndexHit surfHit;
                 label hitSurface;
 
-                geometryToConformTo_.findSurfaceNearest
+                geometryToConformTo_.findSurfaceNearestIntersection
                 (
-                    (
+                    cellCentre,
-                        vA->internalPoint()
+                    pt,
-                      ? topoint(vA->point())
-                      : topoint(vB->point())
-                    ),
-                    magSqr(topoint(vA->point()) - topoint(vB->point())),
                     surfHit,
                     hitSurface
                 );
@@ -610,11 +687,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
 
                     addSurfaceAndEdgeHits
                     (
-                        (
+                        pt,
-                            vA->internalPoint()
-                          ? topoint(vA->point())
-                          : topoint(vB->point())
-                        ),
                         surfaceIntersections,
                         surfacePtReplaceDistCoeffSqr,
                         edgeSearchDistCoeffSqr,
@@ -1613,6 +1686,7 @@ void Foam::conformalVoronoiMesh::limitDisplacement
     // Do not allow infinite recursion
     if (callCount > 7)
     {
+        displacement = vector::zero;
         return;
     }
 
@@ -1661,6 +1735,7 @@ void Foam::conformalVoronoiMesh::limitDisplacement
             if (magSqr(pt - surfHit.hitPoint()) <= searchDistanceSqr)
             {
                 // Cannot limit displacement, point closer than tolerance
+                displacement = vector::zero;
                 return;
             }
         }
@@ -1748,17 +1823,19 @@ bool Foam::conformalVoronoiMesh::nearSurfacePoint
     const Foam::point& pt = pHit.first().hitPoint();
 
     pointIndexHit closePoint;
-
     const bool closeToSurfacePt = pointIsNearSurfaceLocation(pt, closePoint);
 
     if
     (
         closeToSurfacePt
-     && mag(pt - closePoint.hitPoint()) > pointPairDistance(pt)
+     && (
+            magSqr(pt - closePoint.hitPoint())
+          > sqr(pointPairDistance(pt))
+        )
     )
     {
-        const scalar cosAngle
+        const scalar cosAngle =
-            = angleBetweenSurfacePoints(pt, closePoint.hitPoint());
+            angleBetweenSurfacePoints(pt, closePoint.hitPoint());
 
         // @todo make this tolerance run-time selectable?
         if (cosAngle < searchAngleOppositeSurface)
@@ -1768,11 +1845,6 @@ bool Foam::conformalVoronoiMesh::nearSurfacePoint
 
             const scalar searchDist = targetCellSize(closePoint.hitPoint());
 
-            if (searchDist < SMALL)
-            {
-                Pout<< "WARNING: SMALL CELL SIZE" << endl;
-            }
-
             geometryToConformTo_.findSurfaceNearest
             (
                 closePoint.hitPoint(),
@@ -1789,15 +1861,15 @@ bool Foam::conformalVoronoiMesh::nearSurfacePoint
                 norm
             );
 
-            const vector nA = norm[0];
+            const vector& nA = norm[0];
 
             pointIndexHit oppositeHit;
             label oppositeSurfaceHit = -1;
 
             geometryToConformTo_.findSurfaceNearestIntersection
             (
-                closePoint.hitPoint() + SMALL*nA,
+                closePoint.hitPoint() + 0.5*pointPairDistance(pt)*nA,
-                closePoint.hitPoint() + mag(pt - closePoint.hitPoint())*nA,
+                closePoint.hitPoint() + 5*targetCellSize(pt)*nA,
                 oppositeHit,
                 oppositeSurfaceHit
             );
@@ -1808,21 +1880,10 @@ bool Foam::conformalVoronoiMesh::nearSurfacePoint
                 pHit.first() = oppositeHit;
                 pHit.second() = oppositeSurfaceHit;
 
-//                appendToSurfacePtTree(pHit.first().hitPoint());
-//                surfaceToTreeShape.append
-//                (
-//                    existingSurfacePtLocations_.size() - 1
-//                );
-
                 return !closeToSurfacePt;
             }
         }
     }
-    else
-    {
-//        appendToSurfacePtTree(pt);
-//        surfaceToTreeShape.append(existingSurfacePtLocations_.size() - 1);
-    }
 
     return closeToSurfacePt;
 }

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

@@ -230,13 +230,13 @@ inline void Foam::conformalVoronoiMesh::createPointPair
 {
     vector ppDistn = ppDist*n;
 
-    const Foam::point internalPt = surfPt - ppDistn;
+//    const Foam::point internalPt = surfPt - ppDistn;
-    const Foam::point externalPt = surfPt + ppDistn;
+//    const Foam::point externalPt = surfPt + ppDistn;
 
-    bool internalInside = geometryToConformTo_.inside(internalPt);
+//    bool internalInside = geometryToConformTo_.inside(internalPt);
-    bool externalOutside = geometryToConformTo_.outside(externalPt);
+//    bool externalOutside = geometryToConformTo_.outside(externalPt);
 
-    if (internalInside && externalOutside)
+//    if (internalInside && externalOutside)
     {
         pts.append
         (
@@ -259,15 +259,24 @@ inline void Foam::conformalVoronoiMesh::createPointPair
                 Pstream::myProcNo()
             )
         );
+
+//        if (ptPair)
+//        {
+//            ptPairs_.addPointPair
+//            (
+//                pts[pts.size() - 2].index(),
+//                pts[pts.size() - 1].index() // external 0 -> slave
+//            );
+//        }
     }
-    else
+//    else
-    {
+//    {
-        Info<< "Warning: point pair not inside/outside" << nl
+//        Info<< "Warning: point pair not inside/outside" << nl
-            << "    surfPt = " << surfPt << nl
+//            << "    surfPt = " << surfPt << nl
-            << "    internal = " << internalPt << " " << internalInside << nl
+//            << "    internal = " << internalPt << " " << internalInside << nl
-            << "    external = " << externalPt << " " << externalOutside
+//            << "    external = " << externalPt << " " << externalOutside
-            << endl;
+//            << endl;
-    }
+//    }
 }
 
 

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

@@ -183,6 +183,7 @@ public:
         inline bool hasSeedPoint() const;
 
         inline bool hasInternalPoint() const;
+        inline bool hasBoundaryPoint() const;
 
         inline bool hasConstrainedPoint() const;
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H

@@ -228,6 +228,19 @@ inline bool CGAL::indexedCell<Gt, Cb>::hasInternalPoint() const
 }
 
 
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::hasBoundaryPoint() const
+{
+    return
+    (
+        this->vertex(0)->boundaryPoint()
+     || this->vertex(1)->boundaryPoint()
+     || this->vertex(2)->boundaryPoint()
+     || this->vertex(3)->boundaryPoint()
+    );
+}
+
+
 template<class Gt, class Cb>
 inline bool CGAL::indexedCell<Gt, Cb>::hasConstrainedPoint() const
 {

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C

@@ -683,14 +683,18 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide
                 vector hitDir = info[0].rawPoint() - samplePts[i];
                 hitDir /= mag(hitDir) + SMALL;
 
-                if
+                pointIndexHit surfHit;
+                label hitSurface;
+
+                findSurfaceNearestIntersection
                 (
-                    findSurfaceAnyIntersection
+                    samplePts[i],
-                    (
+                    info[0].rawPoint(),
-                        samplePts[i],
+                    surfHit,
-                        info[0].rawPoint() - 1e-3*mag(hitDir)*(hitDir)
+                    hitSurface
-                    )
+                );
-                )
+
+                if (surfHit.hit() && hitSurface != surfaces_[s])
                 {
                     continue;
                 }

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C

@@ -37,6 +37,90 @@ namespace Foam
 defineTypeNameAndDebug(rayShooting, 0);
 addToRunTimeSelectionTable(initialPointsMethod, rayShooting, dictionary);
 
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void rayShooting::splitLine
+(
+    const line<point, point>& l,
+    const scalar& pert,
+    DynamicList<Vb::Point>& initialPoints
+) const
+{
+    Foam::point midPoint(l.centre());
+    const scalar localCellSize(cellShapeControls().cellSize(midPoint));
+    const scalar lineLength(l.mag());
+
+    const scalar minDistFromSurfaceSqr
+    (
+        minimumSurfaceDistanceCoeffSqr_
+       *sqr(localCellSize)
+    );
+
+    if
+    (
+        magSqr(midPoint - l.start()) > minDistFromSurfaceSqr
+     && magSqr(midPoint - l.end()) > minDistFromSurfaceSqr
+    )
+    {
+        // Add extra points if line length is much bigger than local cell size
+//        if (lineLength > 4.0*localCellSize)
+//        {
+//            splitLine
+//            (
+//                line<point, point>(l.start(), midPoint),
+//                pert,
+//                initialPoints
+//            );
+//
+//            splitLine
+//            (
+//                line<point, point>(midPoint, l.end()),
+//                pert,
+//                initialPoints
+//            );
+//        }
+
+        if (randomiseInitialGrid_)
+        {
+            Foam::point newPt
+            (
+                midPoint.x() + pert*(rndGen().scalar01() - 0.5),
+                midPoint.y() + pert*(rndGen().scalar01() - 0.5),
+                midPoint.z() + pert*(rndGen().scalar01() - 0.5)
+            );
+
+            if
+            (
+                !geometryToConformTo().findSurfaceAnyIntersection
+                (
+                    midPoint,
+                    newPt
+                )
+            )
+            {
+                midPoint = newPt;
+            }
+            else
+            {
+                WarningIn
+                (
+                    "rayShooting::splitLine"
+                    "("
+                    "   const line<point,point>&,"
+                    "   const scalar&,"
+                    "   DynamicList<Vb::Point>&"
+                    ")"
+                )   << "Point perturbation crosses a surface. Not inserting."
+                    << endl;
+            }
+        }
+
+        initialPoints.append(toPoint(midPoint));
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 rayShooting::rayShooting
@@ -75,7 +159,7 @@ List<Vb::Point> rayShooting::initialPoints() const
     const searchableSurfaces& surfaces = geometryToConformTo().geometry();
     const labelList& surfacesToConformTo = geometryToConformTo().surfaces();
 
-    const scalar maxRayLength = surfaces.bounds().mag();
+    const scalar maxRayLength(surfaces.bounds().mag());
 
     // Initialise points list
     label initialPointsSize = 0;
@@ -94,8 +178,7 @@ List<Vb::Point> rayShooting::initialPoints() const
         const pointField& faceCentres = faceCentresTmp();
 
         Info<< "    Shoot rays from " << s.name() << nl
-            << "    nRays = " << faceCentres.size() << endl;
+            << "        nRays = " << faceCentres.size() << endl;
-
 
         forAll(faceCentres, fcI)
         {
@@ -110,9 +193,11 @@ List<Vb::Point> rayShooting::initialPoints() const
                 continue;
             }
 
-            const scalar pert =
+            const scalar pert
+            (
                 randomPerturbationCoeff_
-               *cellShapeControls().cellSize(fC);
+               *cellShapeControls().cellSize(fC)
+            );
 
             pointIndexHit surfHitStart;
             label hitSurfaceStart;
@@ -181,27 +266,12 @@ List<Vb::Point> rayShooting::initialPoints() const
                         }
                     }
 
-                    Foam::point midPoint(l.centre());
+                    splitLine
-
-                    const scalar minDistFromSurfaceSqr =
-                        minimumSurfaceDistanceCoeffSqr_
-                       *sqr(cellShapeControls().cellSize(midPoint));
-
-                    if (randomiseInitialGrid_)
-                    {
-                        midPoint.x() += pert*(rndGen().scalar01() - 0.5);
-                        midPoint.y() += pert*(rndGen().scalar01() - 0.5);
-                        midPoint.z() += pert*(rndGen().scalar01() - 0.5);
-                    }
-
-                    if
                     (
-                        magSqr(midPoint - l.start()) > minDistFromSurfaceSqr
+                        l,
-                     && magSqr(midPoint - l.end()) > minDistFromSurfaceSqr
+                        pert,
-                    )
+                        initialPoints
-                    {
+                    );
-                        initialPoints.append(toPoint(midPoint));
-                    }
                 }
             }
         }

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H

@@ -61,6 +61,16 @@ private:
         scalar randomPerturbationCoeff_;
 
 
+    // Private Member Functions
+
+        void splitLine
+        (
+            const line<point, point>& l,
+            const scalar& pert,
+            DynamicList<Vb::Point>& initialPoints
+        ) const;
+
+
 public:
 
     //- Runtime type information