Implemented face rotation controller and some associated functions and fixes.

Some test code code commented out at point of displacementAccumulator
addition/clip.

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

@@ -63,9 +63,15 @@ int main(int argc, char *argv[])
 
         mesh.move();
 
+        mesh.conformToSurface();
+
         runTime++;
     }
 
+    mesh.writeMesh();
+
+    mesh.writePoints("allFinalPoints.obj", false);
+
     Info<< nl << "End\n" << endl;
 
     return 0;

src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C

@@ -95,25 +95,13 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
         )
     )
 {
-    timeCheck();
-
     createFeaturePoints();
-    timeCheck();
 
     insertInitialPoints();
-    timeCheck();
 
     conformToSurface();
-    timeCheck();
 
-    writePoints("allPoints.obj", false);
+    writePoints("allInitialPoints.obj", false);
-    timeCheck();
-
-    writeMesh();
-    timeCheck();
-
-    writeTargetCellSize();
-    timeCheck();
 }
 
 
@@ -256,7 +244,7 @@ Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment
 
     // Auxiliary alignment generated by spoke intersection normal.
 
-    allGeometry_[hitSurface].getNormal
+    allGeometry_[closestSpokeSurface].getNormal
     (
         List<pointIndexHit>(1, closestSpokeHit),
         norm
@@ -270,8 +258,13 @@ Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment
     if (mag(ns) <  SMALL)
     {
         FatalErrorIn("conformalVoronoiMesh::requiredAlignment")
-            << "Parallel normals detected in spoke search."
+        << "Parallel normals detected in spoke search." << nl
-            << nl << exit(FatalError);
+            << "point: " << pt << nl
+            << "closest surface point: " << surfHit.hitPoint() << nl
+            << "closest spoke hit: " << closestSpokeHit.hitPoint() << nl
+            << "np: " << surfHit.hitPoint() + np << nl
+            << "ns: " << closestSpokeHit.hitPoint() + na << nl
+            << exit(FatalError);
     }
 
     ns /= mag(ns);
@@ -781,6 +774,8 @@ void Foam::conformalVoronoiMesh::constructFeaturePointLocations()
 
 void Foam::conformalVoronoiMesh::reinsertFeaturePoints()
 {
+    Info<< nl << "    Reinserting stored feature points" << endl;
+
     forAll(featureVertices_, f)
     {
         insertVb(featureVertices_[f]);
@@ -838,31 +833,6 @@ void Foam::conformalVoronoiMesh::insertInitialPoints()
 
     insertPoints(initialPointsMethod_->initialPoints());
 
-    // std::vector<Point> initialPoints = initialPointsMethod_->initialPoints();
-
-    // Info<< "    " << initialPoints.size() << " points to insert..." << endl;
-
-    // label nVert = startOfInternalPoints_;
-
-    // // using the range insert (faster than inserting points one by one)
-    // insert(initialPoints.begin(), initialPoints.end());
-
-    // Info<< "    " << number_of_vertices() - startOfInternalPoints_
-    //     << " points inserted" << endl;
-
-    // for
-    // (
-    //     Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
-    //     vit != finite_vertices_end();
-    //     ++vit
-    // )
-    // {
-    //     if (vit->uninitialised())
-    //     {
-    //         vit->index() = nVert++;
-    //     }
-    // }
-
     writePoints("initialPoints.obj", true);
 }
 
@@ -1791,6 +1761,8 @@ void Foam::conformalVoronoiMesh::conformToSurface()
 {
     Info<< nl << "Conforming to surfaces" << endl;
 
+    timeCheck();
+
     startOfSurfacePointPairs_ = number_of_vertices();
 
     // Initialise containers to store the edge conformation locations
@@ -1887,7 +1859,7 @@ void Foam::conformalVoronoiMesh::conformToSurface()
 
     label iterationNo = 0;
 
-    label maxIterations = 10;
+    label maxIterations = 5;
     Info << "    MAX ITERATIONS HARD CODED TO "<< maxIterations << endl;
 
     // Set totalHits to a positive value to enter the while loop on the first
@@ -1997,12 +1969,6 @@ void Foam::conformalVoronoiMesh::move()
 
     pointField dualVertices(number_of_cells());
 
-    pointField displacementAccumulator(startOfSurfacePointPairs_, vector::zero);
-
-    List<bool> pointToBeRetained(startOfSurfacePointPairs_, true);
-
-    DynamicList<point> newPointsToInsert;
-
     label dualVerti = 0;
 
     // Find the dual point of each tetrahedron and assign it an index.
@@ -2033,8 +1999,6 @@ void Foam::conformalVoronoiMesh::move()
 
     dualVertices.setSize(dualVerti);
 
-    timeCheck();
-
     Info<< nl << "    Calculating target cell alignment and size" << endl;
 
     for
@@ -2054,8 +2018,6 @@ void Foam::conformalVoronoiMesh::move()
         }
     }
 
-    timeCheck();
-
     Info<< nl << "    Looping over all dual faces" << endl;
 
     vectorField cartesianDirections(3);
@@ -2064,6 +2026,18 @@ void Foam::conformalVoronoiMesh::move()
     cartesianDirections[1] = vector(0,1,0);
     cartesianDirections[2] = vector(1,0,0);
 
+    vectorField displacementAccumulator
+    (
+        startOfSurfacePointPairs_,
+        vector::zero
+    );
+
+    PackedBoolList pointToBeRetained(startOfSurfacePointPairs_, true);
+
+    std::vector<Point> pointsToInsert;
+
+    label pointsAdded = 0;
+
     for
     (
         Triangulation::Finite_edges_iterator eit = finite_edges_begin();
@@ -2101,6 +2075,8 @@ void Foam::conformalVoronoiMesh::move()
 
             verticesOnFace.shrink();
 
+            face dualFace(verticesOnFace);
+
             Cell_handle c = eit->first;
             Vertex_handle vA = c->vertex(eit->second);
             Vertex_handle vB = c->vertex(eit->third);
@@ -2169,11 +2145,176 @@ void Foam::conformalVoronoiMesh::move()
                     alignmentDir *= 0.5*targetCellSize;
 
                     vector delta = alignmentDir - 0.5*rAB;
+
+                    scalar faceArea = dualFace.mag(dualVertices);
+
+                    delta *= faceAreaWeightModel_->faceAreaWeight
+                    (
+                        faceArea/targetFaceArea
+                    );
+
+                    if
+                    (
+                        vA->internalPoint()
+                     && vB->internalPoint()
+                     && rABMag
+                          > cvMeshControls().insertionDistCoeff()*targetCellSize
+                     && faceArea
+                          > cvMeshControls().faceAreaRatioCoeff()*targetFaceArea
+                     && alignmentDotProd
+                          > cvMeshControls().cosInsertionAcceptanceAngle()
+                    )
+                    {
+                        // Point insertion
+                        pointsToInsert.push_back
+                        (
+                            toPoint(0.5*(dVA + dVB))
+                        );
+
+                        pointsAdded++;
+                    }
+                    else if
+                    (
+                        (vA->internalPoint() || vB->internalPoint())
+                     && rABMag
+                          < cvMeshControls().removalDistCoeff()*targetCellSize
+                    )
+                    {
+                        // Point removal
+
+                        // Only insert a point at the midpoint of the short edge
+                        // if neither attached point has already been identified
+                        // to be removed.
+                        if
+                        (
+                            pointToBeRetained[vA->index()] == true
+                         && pointToBeRetained[vB->index()] == true
+                        )
+                        {
+                            pointsToInsert.push_back
+                            (
+                                toPoint(0.5*(dVA + dVB))
+                            );
+                        }
+
+                        if (vA->internalPoint())
+                        {
+                            pointToBeRetained[vA->index()] = false;
+                        }
+
+                        if (vB->internalPoint())
+                        {
+                            pointToBeRetained[vB->index()] = false;
+                        }
+                    }
+                    else
+                    {
+                        if (vA->internalPoint())
+                        {
+                            displacementAccumulator[vA->index()] += delta;
+                        }
+                        if (vB->internalPoint())
+                        {
+                            displacementAccumulator[vB->index()] += -delta;
                         }
                     }
                 }
             }
         }
+    }
+
+    // Clip displacements that pierce, or get too close to the surface
+
+    Info<< "ARBITRARY CLIP TO 2*aveDisp" << endl;
+
+    scalar aveDisp =
+        sum(mag(displacementAccumulator))/displacementAccumulator.size();
+
+    forAll(displacementAccumulator, i)
+    {
+        if (mag(displacementAccumulator[i]) > 2*aveDisp)
+        {
+            displacementAccumulator[i] *=
+                2*aveDisp/mag(displacementAccumulator[i]);
+
+            // for
+            // (
+            //     Triangulation::Finite_vertices_iterator vit =
+            //         finite_vertices_begin();
+            //     vit != finite_vertices_end();
+            //     ++vit
+            // )
+            // {
+            //     if (vit->index() == i)
+            //     {
+            //         Info<< topoint(vit->point()) << nl
+            //             << topoint(vit->point()) + displacementAccumulator[i]
+            //             <<endl;
+            //     }
+            // }
+        }
+    }
+
+    // scalarField bins(100);
+
+    // scalar maxDisp = max(mag(displacementAccumulator));
+
+    // forAll(bins, i)
+    // {
+    //     bins[i] = i*1.01*maxDisp/bins.size();
+    // }
+
+    // Histogram<scalarField> dispHisto(bins, mag(displacementAccumulator));
+
+    // Info<< nl << bins << nl << dispHisto.counts() << endl;
+
+    vector totalDisp = sum(displacementAccumulator);
+    scalar totalDist = sum(mag(displacementAccumulator));
+
+    Info<< nl
+        << "    Total displacement = " << totalDisp << nl
+        << "    Total distance = " << totalDist << nl
+        << "    Points added = " << pointsAdded
+        << endl;
+
+    // Relax the calculated displacement
+    displacementAccumulator *= relaxation;
+
+    for
+    (
+        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        ++vit
+    )
+    {
+        if (vit->internalPoint())
+        {
+            if (pointToBeRetained[vit->index()] == true)
+            {
+                pointsToInsert.push_back
+                (
+                    vit->point()
+                  + toCGALVector(displacementAccumulator[vit->index()])
+                );
+            }
+        }
+    }
+
+    // Write the mesh before clearing it
+    if (runTime_.outputTime())
+    {
+        writeMesh(false);
+    }
+
+    // Remove the entire tessellation
+    this->clear();
+
+    reinsertFeaturePoints();
+    startOfInternalPoints_ = number_of_vertices();
+
+    Info<< nl << "Reinserting entire tessellation" << endl;
+    insertPoints(pointsToInsert);
+}
 
 
 // ************************************************************************* //

src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H

@@ -48,6 +48,7 @@ SourceFiles
 #include "cellSizeControlSurfaces.H"
 #include "cvControls.H"
 #include "DynamicList.H"
+#include "PackedBoolList.H"
 #include "Time.H"
 #include "polyMesh.H"
 #include "plane.H"
@@ -60,6 +61,7 @@ SourceFiles
 #include "transform.H"
 #include "volFields.H"
 #include "fvMesh.H"
+#include "Histogram.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -172,7 +174,8 @@ class conformalVoronoiMesh
             const label type
         );
 
-        //- Insert a range of points using the CGAL range insertion method
+        //- Insert a std::vector of CGAL Points using the CGAL range insertion
+        //  method
         inline void insertPoints(const std::vector<Point>& points);
 
         //- Insert a point-pair at a ppDist distance either side of
@@ -386,6 +389,9 @@ public:
         inline pointFromPoint topoint(const Point&) const;
         inline PointFrompoint toPoint(const point&) const;
 
+        typedef K::Vector_3 CGALVector;
+
+        inline CGALVector toCGALVector(const point& pt) const;
 
         // Access
 
@@ -414,8 +420,12 @@ public:
                 const List<point>& points
             ) const;
 
+            //- Write the internal Delaunay vertices of the tessellation as a
+            //  pointField that may be used to restart the meshing process
+            void writeInternalDelaunayVertices() const;
+
             //- Write polyMesh
-            void writeMesh();
+            void writeMesh(bool writeToConstant = true);
 
             //- Write dual points and faces as .obj file
             void writeDual

src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H

@@ -158,7 +158,7 @@ inline void Foam::conformalVoronoiMesh::insertPoints
     const std::vector<Point>& points
 )
 {
-    Info<< "    " << points.size() << " points to insert..." << endl;
+    Info<< nl << "    " << points.size() << " points to insert..." << endl;
 
     label nVert = number_of_vertices();
 
@@ -277,6 +277,13 @@ Foam::conformalVoronoiMesh::toPoint
 #endif
 
 
+inline Foam::conformalVoronoiMesh::CGALVector
+Foam::conformalVoronoiMesh::toCGALVector(const point& pt) const
+{
+    return CGALVector(pt.x(), pt.y(), pt.z());
+}
+
+
 inline const Foam::Time& Foam::conformalVoronoiMesh::time() const
 {
     return runTime_;

src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C

@@ -77,7 +77,48 @@ void Foam::conformalVoronoiMesh::writePoints
 }
 
 
-void Foam::conformalVoronoiMesh::writeMesh()
+void Foam::conformalVoronoiMesh::writeInternalDelaunayVertices() const
+{
+    Info<< nl << "    Writing internal Delaunay vertices to pointField "
+        << "ADD NAME CHOICE ARGUMENT" << endl;
+
+    pointField internalDelaunayVertices(number_of_vertices());
+
+    label vertI = 0;
+
+    for
+    (
+        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        ++vit
+    )
+    {
+        if (vit->internalPoint())
+        {
+            internalDelaunayVertices[vertI++] = topoint(vit->point());
+        }
+    }
+
+    internalDelaunayVertices.setSize(vertI);
+
+    pointIOField internalDVs
+    (
+        IOobject
+        (
+            "internalDelaunayVertices",
+            runTime_.timeName(),
+            runTime_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        internalDelaunayVertices
+    );
+
+    internalDVs.write();
+}
+
+
+void Foam::conformalVoronoiMesh::writeMesh(bool writeToConstant)
 {
     pointField points(0);
     faceList faces(0);
@@ -87,6 +128,8 @@ void Foam::conformalVoronoiMesh::writeMesh()
     labelList patchSizes(0);
     labelList patchStarts(0);
 
+    writeInternalDelaunayVertices();
+
     calcDualMesh
     (
         points,
@@ -101,6 +144,19 @@ void Foam::conformalVoronoiMesh::writeMesh()
     writeDual(points, faces, "dualMesh.obj");
 
     IOobject io
+    (
+        Foam::polyMesh::defaultRegion,
+        runTime_.timeName(),
+        runTime_,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    );
+
+    if (writeToConstant)
+    {
+        Info<< nl << "    Writing polyMesh to constant." << endl;
+
+        io = IOobject
         (
             Foam::polyMesh::defaultRegion,
             runTime_.constant(),
@@ -109,7 +165,12 @@ void Foam::conformalVoronoiMesh::writeMesh()
             IOobject::AUTO_WRITE
         );
 
-    Info<< nl << "Writing polyMesh to constant." << endl;
+    }
+    else
+    {
+        Info<< nl << "    Writing polyMesh to time directory "
+            << runTime_.timeName() << endl;
+    }
 
     polyMesh pMesh
     (
@@ -142,6 +203,8 @@ void Foam::conformalVoronoiMesh::writeMesh()
         << "Failed writing polyMesh."
             << exit(FatalError);
     }
+
+    writeTargetCellSize();
 }
 
 

src/conformalVoronoiMesh/cvControls/cvControls.C

@@ -97,6 +97,7 @@ Foam::cvControls::cvControls
     cosAlignmentAcceptanceAngle_ = cos
     (
         readScalar(motionDict.lookup("alignmentAcceptanceAngle"))
+       *mathematicalConstant::pi/180
     );
 
     // Point removal criteria
@@ -106,7 +107,7 @@ Foam::cvControls::cvControls
         motionDict.subDict("pointInsertionCriteria")
     );
 
-    cellCentreInsertionDistCoeff_ = readScalar
+    insertionDistCoeff_ = readScalar
     (
         insertionDict.lookup("cellCentreDistCoeff")
     );
@@ -119,6 +120,7 @@ Foam::cvControls::cvControls
     cosInsertionAcceptanceAngle_ = cos
     (
         readScalar(insertionDict.lookup("acceptanceAngle"))
+       *mathematicalConstant::pi/180
     );
 
     // Point removal criteria
@@ -128,7 +130,7 @@ Foam::cvControls::cvControls
         motionDict.subDict("pointRemovalCriteria")
     );
 
-    cellCentreRemovalDistCoeff_ = readScalar
+    removalDistCoeff_ = readScalar
     (
         removalDict.lookup("cellCentreDistCoeff")
     );

src/conformalVoronoiMesh/cvControls/cvControls.H

@@ -116,7 +116,7 @@ class cvControls
 
             //- Length between Delaunay vertices above which a new Dv should be
             //  inserted - fraction of the local target cell size
-            scalar cellCentreInsertionDistCoeff_;
+            scalar insertionDistCoeff_;
 
             //- Minimum dual face area corresponding to long Delaunay edge where
             //  a new Dv is to be inserted - fraction of the local target cell
@@ -132,7 +132,7 @@ class cvControls
 
             //- Length between Delaunay vertices below which a Dv should be
             //  removed - fraction of the local target cell size
-            scalar cellCentreRemovalDistCoeff_;
+            scalar removalDistCoeff_;
 
 
     // polyMesh filtering controls
@@ -198,8 +198,8 @@ public:
             //- Return the cosAlignmentAcceptanceAngle
             inline scalar cosAlignmentAcceptanceAngle() const;
 
-            //- Return the cellCentreInsertionDistCoeff
+            //- Return the insertionDistCoeff
-            inline scalar cellCentreInsertionDistCoeff() const;
+            inline scalar insertionDistCoeff() const;
 
             //- Return the faceAreaRatioCoeff
             inline scalar faceAreaRatioCoeff() const;
@@ -207,8 +207,8 @@ public:
             //- Return the cosInsertionAcceptanceAngle
             inline scalar cosInsertionAcceptanceAngle() const;
 
-            //- Return cellCentreRemovalDistCoeff
+            //- Return removalDistCoeff
-            inline scalar cellCentreRemovalDistCoeff() const;
+            inline scalar removalDistCoeff() const;
 
             //- Return the minimumEdgeLengthCoeff
             inline scalar minimumEdgeLengthCoeff() const;

src/conformalVoronoiMesh/cvControls/cvControlsI.H

@@ -97,9 +97,9 @@ inline Foam::scalar Foam::cvControls::cosAlignmentAcceptanceAngle() const
     return cosAlignmentAcceptanceAngle_;
 }
 
-inline Foam::scalar Foam::cvControls::cellCentreInsertionDistCoeff() const
+inline Foam::scalar Foam::cvControls::insertionDistCoeff() const
 {
-    return cellCentreInsertionDistCoeff_;
+    return insertionDistCoeff_;
 }
 
 inline Foam::scalar Foam::cvControls::faceAreaRatioCoeff() const
@@ -112,9 +112,9 @@ inline Foam::scalar Foam::cvControls::cosInsertionAcceptanceAngle() const
     return cosInsertionAcceptanceAngle_;
 }
 
-inline Foam::scalar Foam::cvControls::cellCentreRemovalDistCoeff() const
+inline Foam::scalar Foam::cvControls::removalDistCoeff() const
 {
-    return cellCentreRemovalDistCoeff_;
+    return removalDistCoeff_;
 }
 
 inline Foam::scalar Foam::cvControls::minimumEdgeLengthCoeff() const