From 1f61715738cf805bbbf9b5c126060aba5c802c22 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Thu, 16 Jun 2011 14:57:54 +0100
Subject: [PATCH] ENH: Distributing and storing feature points, surface
 conformation and sizes.

---
 .../conformalVoronoiMesh.C                    | 178 +++++++++++++-----
 .../conformalVoronoiMesh.H                    |  40 +++-
 .../conformalVoronoiMeshConformToSurface.C    |  51 ++++-
 .../conformalVoronoiMeshI.H                   |  60 +++---
 .../conformalVoronoiMesh/indexedVertex.H      |  67 +++----
 5 files changed, 266 insertions(+), 130 deletions(-)

diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 5157f3c580..208dba1083 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -283,7 +283,7 @@ void Foam::conformalVoronoiMesh::insertPoints
     //     ++pit
     // )
     // {
-    //     insertVb(*pit);
+    //     insertPoint(topoint(*pit), Vb::ptInternalPoint);
     // }
 
     label nInserted(number_of_vertices() - preInsertionSize);
@@ -318,6 +318,10 @@ void Foam::conformalVoronoiMesh::insertPoints
     bool distribute
 )
 {
+    // The pts, indices and types lists must be intact and up-to-date at the
+    // end of this function as they may also be used by other functions
+    // subsequently.
+
     if (Pstream::parRun() && distribute)
     {
         // The link between vertices that form the boundary via pairs cannot be
@@ -363,7 +367,6 @@ void Foam::conformalVoronoiMesh::insertPoints
         if (type > Vb::ptFarPoint)
         {
             // This is a member of a point pair, don't use the type directly
-
             type += number_of_vertices();
         }
 
@@ -750,6 +753,7 @@ void Foam::conformalVoronoiMesh::insertFeaturePoints()
 
     createMixedFeaturePoints(pts, indices, types);
 
+    // Insert the created points, distributing to the appropriate processor
     insertPoints(pts, indices, types, true);
 
     if(cvMeshControls().objOutput())
@@ -770,27 +774,15 @@ void Foam::conformalVoronoiMesh::insertFeaturePoints()
             << " feature vertices" << endl;
     }
 
-    Info<< "SORT OUT FEATURE POINT DISTRIBUTION AND STORAGE" << endl;
+    featureVertices_.clear();
 
-    // featureVertices_.setSize(number_of_vertices());
-
-    // label featPtI = 0;
-
-    // for
-    // (
-    //     Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-    //     vit != finite_vertices_end();
-    //     vit++
-    // )
-    // {
-    //     featureVertices_[featPtI] = Vb(vit->point());
-
-    //     featureVertices_[featPtI].index() = vit->index();
-
-    //     featureVertices_[featPtI].type() = vit->type();
-
-    //     featPtI++;
-    // }
+    forAll(pts, pI)
+    {
+        featureVertices_.push_back
+        (
+            Vb(toPoint(pts[pI]), indices[pI], types[pI])
+        );
+    }
 
     constructFeaturePointLocations();
 }
@@ -1066,14 +1058,72 @@ void Foam::conformalVoronoiMesh::constructFeaturePointLocations()
 }
 
 
-void Foam::conformalVoronoiMesh::reinsertFeaturePoints()
+void Foam::conformalVoronoiMesh::reinsertFeaturePoints(bool distribute)
 {
     Info<< nl << "Reinserting stored feature points" << endl;
 
-    forAll(featureVertices_, f)
+    label preReinsertionSize(number_of_vertices());
+
+    if (distribute)
     {
-        insertVb(featureVertices_[f]);
+        DynamicList<Foam::point> pointsToInsert;
+        DynamicList<label> indices;
+        DynamicList<label> types;
+
+        for
+        (
+            std::list<Vb>::iterator vit=featureVertices_.begin();
+            vit != featureVertices_.end();
+            ++vit
+        )
+        {
+            pointsToInsert.append(topoint(vit->point()));
+            indices.append(vit->index());
+            types.append(vit->type());
+        }
+
+        // Insert distributed points
+        insertPoints(pointsToInsert, indices, types, true);
+
+        // Save points in new distribution
+        featureVertices_.clear();
+
+        forAll(pointsToInsert, pI)
+        {
+            featureVertices_.push_back
+            (
+                Vb(toPoint(pointsToInsert[pI]), indices[pI], types[pI])
+            );
+        }
     }
+    else
+    {
+        for
+        (
+            std::list<Vb>::iterator vit=featureVertices_.begin();
+            vit != featureVertices_.end();
+            ++vit
+        )
+        {
+            // Assuming that all of the reinsertions are pair points, and that
+            // the index and type are relative, i.e. index 0 and type relative
+            // to it.
+            insertPoint
+            (
+                vit->point(),
+                vit->index() + number_of_vertices(),
+                vit->type() + number_of_vertices()
+            );
+        }
+    }
+
+    Info<< "    Reinserted "
+        << returnReduce
+        (
+            label(number_of_vertices()) - preReinsertionSize,
+            sumOp<label>()
+        )
+        << " vertices" << endl;
 }
 
 
@@ -1195,7 +1245,8 @@ bool Foam::conformalVoronoiMesh::distributeBackground()
             vit++
         )
         {
-            if (vit->real())
+            // Only store real vertices that are not feature vertices
+            if (vit->real() && vit->index() >= startOfInternalPoints_)
             {
                 nRealVertices++;
             }
@@ -1259,7 +1310,8 @@ bool Foam::conformalVoronoiMesh::distributeBackground()
             vit++
         )
         {
-            if (vit->real())
+            // Only store real vertices that are not feature vertices
+            if (vit->real() && vit->index() >= startOfInternalPoints_)
             {
                 Foam::point v = topoint(vit->point());
 
@@ -1308,8 +1360,8 @@ bool Foam::conformalVoronoiMesh::distributeBackground()
         // Remove the entire tessellation
         reset();
 
-        Info<< "NEED TO MAP FEATURE POINTS" << endl;
-        // reinsertFeaturePoints();
+        // Reinsert feature points, distributing them as necessary.
+        reinsertFeaturePoints(true);
 
         startOfInternalPoints_ = number_of_vertices();
 
@@ -1326,6 +1378,8 @@ bool Foam::conformalVoronoiMesh::distributeBackground()
             forAll(cellVertices[cI], cVPI)
             {
                 pointsToInsert.append(cellVertices[cI][cVPI]);
+
+                // All insertions relative to index of zero
                 indices.append(0);
 
                 label type = cellVertexTypes[cI][cVPI];
@@ -1369,6 +1423,27 @@ bool Foam::conformalVoronoiMesh::distributeBackground()
 }
 
 
+void Foam::conformalVoronoiMesh::storeSizesAndAlignments()
+{
+    std::list<Point> storePts;
+
+    for
+    (
+        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        vit++
+    )
+    {
+        if (vit->internalPoint())
+        {
+            storePts.push_back(vit->point());
+        }
+    }
+
+    storeSizesAndAlignments(storePts);
+}
+
+
 void Foam::conformalVoronoiMesh::storeSizesAndAlignments
 (
     const std::list<Point>& storePts
@@ -1427,6 +1502,8 @@ void Foam::conformalVoronoiMesh::updateSizesAndAlignments
     )
     {
         storeSizesAndAlignments(storePts);
+
+        timeCheck("Updated sizes and alignments");
     }
 }
 
@@ -1794,6 +1871,8 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
 
     insertBoundingPoints();
 
+    insertFeaturePoints();
+
     insertInitialPoints();
 
     // Improve the guess that the backgroundMeshDecomposition makes with the
@@ -1801,15 +1880,6 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
     // better balance the surface conformation load.
     distributeBackground();
 
-    insertFeaturePoints();
-
-    Info<< "NEED TO MAP FEATURE POINTS AFTER DISTRIBUTE" << endl;
-    Info<< "NEED TO ENSURE THAT FEATURE POINTS ARE INSERTED ON THE "
-        << "CORRECT PROCESSOR" << endl;
-
-    Info<< "NEED TO CHANGE storeSizesAndAlignments AFTER DISTRIBUTE" << endl;
-    // storeSizesAndAlignments(initPts);
-
     buildSurfaceConformation(rmCoarse);
 
     // The introduction of the surface conformation may have distorted the
@@ -1822,6 +1892,21 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
         buildParallelInterface("rebuild");
     }
 
+    // Do not store the surface conformation until after it has been
+    // (potentially) redistributed.
+    storeSurfaceConformation();
+
+    // Use storeSizesAndAlignments with no feed points because all background
+    // points may have been distributed.  It is a requirement that none of the
+    // preceding functions requires look up of sizes or alignments from the
+    // Delaunay vertices, i.e. setVertexSizeAndAlignment cannot be called
+    // before this point.
+    storeSizesAndAlignments();
+
+    // Report any Delaunay vertices that do not think that they are in the
+    // domain the processor they are on.
+    reportProcessorOccupancy();
+
     if(cvMeshControls().objOutput())
     {
         writePoints("allInitialPoints.obj", false);
@@ -2174,13 +2259,6 @@ void Foam::conformalVoronoiMesh::move()
         }
     }
 
-    // Write the mesh before clearing it.  Beware that writeMesh destroys the
-    // indexing of the tessellation.  Do not filter the dual faces.
-    if (runTime_.outputTime())
-    {
-        writeMesh(runTime_.timeName(), false);
-    }
-
     // Remove the entire tessellation
     reset();
 
@@ -2207,14 +2285,20 @@ void Foam::conformalVoronoiMesh::move()
 
     updateSizesAndAlignments(pointsToInsert);
 
+    buildParallelInterface("move_" + runTime_.timeName());
+
+    // Write the intermediate mesh, do not filter the dual faces.
+    if (runTime_.outputTime())
+    {
+        writeMesh(runTime_.timeName(), false);
+    }
+
     Info<< nl
         << "Total displacement = " << totalDisp << nl
         << "Total distance = " << totalDist << nl
         << "Points added = " << pointsAdded << nl
         << "Points removed = " << pointsRemoved
         << endl;
-
-    timeCheck("Updated sizes and alignments (if required), end of move");
 }
 
 
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index ffc87006d6..8a83f44a41 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -138,7 +138,7 @@ private:
 
         //- Store the feature constraining points to be reinserted after a
         //  triangulation clear
-        List<Vb> featureVertices_;
+        std::list<Vb> featureVertices_;
 
         //- Storing the locations of all of the features to be conformed to.
         //  Single pointField required by the featurePointTree.
@@ -240,17 +240,32 @@ private:
         //- Return the required alignment directions at the given location
         tensor requiredAlignment(const Foam::point& pt) const;
 
-        //- Insert point and return its index
+        //- Insert Foam::point and return its auto-generated index
         inline label insertPoint
         (
-            const Foam::point& pt,
+            const Foam::point& p,
             const label type
         );
 
-        //- Insert point and return its index
+        //- Insert Point and return its auto-generated index
+        inline label insertPoint
+        (
+            const Point& P,
+            const label type
+        );
+
+        //- Insert Foam::point with specified index and type
         inline void insertPoint
         (
-            const Foam::point& pt,
+            const Foam::point& p,
+            const label index,
+            const label type
+        );
+
+        //- Insert Point with specified index and type
+        inline void insertPoint
+        (
+            const Point& P,
             const label index,
             const label type
         );
@@ -295,10 +310,6 @@ private:
             DynamicList<label>& types
         );
 
-        //- Insert a Vb (a typedef of CGAL::indexedVertex<K>) with the
-        //  possibility of an offset for the index and the type
-        inline void insertVb(const Vb& v, label offset = 0);
-
         //- Insert pairs of points on the surface with the given normals, at the
         //  specified spacing
         void insertSurfacePointPairs
@@ -420,7 +431,7 @@ private:
         void constructFeaturePointLocations();
 
         //- Reinsert stored feature point defining points
-        void reinsertFeaturePoints();
+        void reinsertFeaturePoints(bool distribute = false);
 
         //- Demand driven construction of octree for feature points
         const indexedOctree<treeDataPoint>& featurePointTree() const;
@@ -445,6 +456,12 @@ private:
         //  referred vertices, so the parallel interface may need rebuilt.
         bool distributeBackground();
 
+        //- Store data for sizeAndAlignmentLocations_, storedSizes_ and
+        //  storedAlignments_ and initialise the sizeAndAlignmentTreePtr_,
+        //  determining the appropriate sizeAndAlignmentLocations_
+        //  automatically
+        void storeSizesAndAlignments();
+
         //- Store data for sizeAndAlignmentLocations_, storedSizes_ and
         //  storedAlignments_ and initialise the sizeAndAlignmentTreePtr_
         void storeSizesAndAlignments(const std::list<Point>& storePts);
@@ -587,6 +604,9 @@ private:
             label& hitSurfaceLargest
         ) const;
 
+        //- Write out vertex-processor occupancy information for debugging
+        void reportProcessorOccupancy();
+
         //- Write out debugging information about the surface conformation
         //  quality
         void reportSurfaceConformationQuality();
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index 1b88564fa6..e129296adc 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -41,7 +41,11 @@ void Foam::conformalVoronoiMesh::conformToSurface()
     {
         // Rebuild, insert and store new surface conformation
         buildSurfaceConformation(reconfMode);
+
+        storeSurfaceConformation();
     }
+
+    // reportSurfaceConformationQuality();
 }
 
 
@@ -411,10 +415,6 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
                 << "), stopping iterations" << endl;
         }
     }
-
-    // reportSurfaceConformationQuality();
-
-    storeSurfaceConformation();
 }
 
 
@@ -620,6 +620,8 @@ void Foam::conformalVoronoiMesh::buildParallelInterface
         receivedVertices,
         outputName
     );
+
+    timeCheck("After buildParallelInterface");
 }
 
 
@@ -1346,6 +1348,27 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
 }
 
 
+void Foam::conformalVoronoiMesh::reportProcessorOccupancy()
+{
+    for
+    (
+        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        vit++
+    )
+    {
+        if (vit->real())
+        {
+            if (!positionOnThisProc(topoint(vit->point())))
+            {
+                Pout<< topoint(vit->point()) << " is not on this processor "
+                    << endl;
+            }
+        }
+    }
+}
+
+
 void Foam::conformalVoronoiMesh::reportSurfaceConformationQuality()
 {
     Info<< nl << "Check surface conformation quality" << endl;
@@ -1605,6 +1628,14 @@ void Foam::conformalVoronoiMesh::buildEdgeLocationTree
 
 void Foam::conformalVoronoiMesh::buildSizeAndAlignmentTree() const
 {
+    if (sizeAndAlignmentLocations_.empty())
+    {
+        FatalErrorIn("buildSizeAndAlignmentTree()")
+            << "sizeAndAlignmentLocations empty, must be populated before "
+            << "sizeAndAlignmentTree can be built."
+            << exit(FatalError);
+    }
+
     treeBoundBox overallBb
     (
         geometryToConformTo_.globalBounds().extend(rndGen_, 1e-4)
@@ -1785,6 +1816,9 @@ void Foam::conformalVoronoiMesh::reinsertSurfaceConformation()
 
     label preReinsertionSize(number_of_vertices());
 
+    // It is assumed that the stored surface conformation is on the correct
+    // processor and does not need distributed
+
     for
     (
         std::list<Vb>::iterator vit=surfaceConformationVertices_.begin();
@@ -1792,7 +1826,14 @@ void Foam::conformalVoronoiMesh::reinsertSurfaceConformation()
         ++vit
     )
     {
-        insertVb(*vit, number_of_vertices());
+        // Assuming that all of the reinsertions are pair points, and that the
+        // index and type are relative, i.e. index 0 and type relative to it.
+        insertPoint
+        (
+            vit->point(),
+            vit->index() + number_of_vertices(),
+            vit->type() + number_of_vertices()
+        );
     }
 
     Info<< "    Reinserted "
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
index f620c7cd08..140c1d7070 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
@@ -215,14 +215,24 @@ inline Foam::label Foam::conformalVoronoiMesh::insertPoint
     const Foam::point& p,
     const label type
 )
+{
+    return insertPoint(toPoint(p), type);
+}
+
+
+inline Foam::label Foam::conformalVoronoiMesh::insertPoint
+(
+    const Point& P,
+    const label type
+)
 {
     uint nVert = number_of_vertices();
 
-    Vertex_handle vh = insert(toPoint(p));
+    Vertex_handle vh = insert(P);
 
     if (nVert == number_of_vertices())
     {
-        Pout << "Failed to insert point " << p << endl;
+        Pout << "Failed to insert point " << topoint(P) << endl;
     }
     else
     {
@@ -240,14 +250,25 @@ inline void Foam::conformalVoronoiMesh::insertPoint
     const label index,
     const label type
 )
+{
+    insertPoint(toPoint(p), index, type);
+}
+
+
+inline void Foam::conformalVoronoiMesh::insertPoint
+(
+    const Point& P,
+    const label index,
+    const label type
+)
 {
     uint nVert = number_of_vertices();
 
-    Vertex_handle vh = insert(toPoint(p));
+    Vertex_handle vh = insert(P);
 
     if (nVert == number_of_vertices())
     {
-        Pout << "Failed to insert point " << p << endl;
+        Pout << "Failed to insert point " << topoint(P) << endl;
     }
     else
     {
@@ -311,37 +332,6 @@ inline void Foam::conformalVoronoiMesh::createPointPair
 }
 
 
-inline void Foam::conformalVoronoiMesh::insertVb(const Vb& v, label offset)
-{
-    const Point& Pt(v.point());
-
-    uint nVert = number_of_vertices();
-
-    Vertex_handle vh = insert(Pt);
-
-    if (nVert == number_of_vertices())
-    {
-        FatalErrorIn("Foam::conformalVoronoiMesh::insertVb(const Vb& v")
-            << "Failed to reinsert Vb at " << topoint(Pt)
-            << nl
-            << exit(FatalError);
-    }
-    else
-    {
-        vh->index() = v.index() + offset;
-
-        if (v.pairPoint())
-        {
-            vh->type() = v.type() + offset;
-        }
-        else
-        {
-            vh->type() = v.type();
-        }
-    }
-}
-
-
 inline bool Foam::conformalVoronoiMesh::isBoundaryDualFace
 (
     const Delaunay::Finite_edges_iterator& eit
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
index f423d17cf6..492cf10889 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
@@ -342,6 +342,7 @@ public:
         return internalPoint(type) || ppMaster(index, type);
     }
 
+
     //- Is point near the boundary or part of the boundary definition
     inline bool nearOrOnBoundary() const
     {
@@ -349,43 +350,43 @@ public:
     }
 
 
-    //- Do the two given vertices consitute a boundary point-pair
-    inline friend bool pointPair
-    (
-        const indexedVertex& v0,
-        const indexedVertex& v1
-    )
-    {
-        return v0.index_ == v1.type_ || v1.index_ == v0.type_;
-    }
+    // //- Do the two given vertices consitute a boundary point-pair
+    // inline friend bool pointPair
+    // (
+    //     const indexedVertex& v0,
+    //     const indexedVertex& v1
+    // )
+    // {
+    //     return v0.index_ == v1.type_ || v1.index_ == v0.type_;
+    // }
 
 
-    //- Do the three given vertices consitute a boundary triangle
-    inline friend bool boundaryTriangle
-    (
-        const indexedVertex& v0,
-        const indexedVertex& v1,
-        const indexedVertex& v2
-    )
-    {
-        return (v0.pairPoint() && pointPair(v1, v2))
-            || (v1.pairPoint() && pointPair(v2, v0))
-            || (v2.pairPoint() && pointPair(v0, v1));
-    }
+    // //- Do the three given vertices consitute a boundary triangle
+    // inline friend bool boundaryTriangle
+    // (
+    //     const indexedVertex& v0,
+    //     const indexedVertex& v1,
+    //     const indexedVertex& v2
+    // )
+    // {
+    //     return (v0.pairPoint() && pointPair(v1, v2))
+    //         || (v1.pairPoint() && pointPair(v2, v0))
+    //         || (v2.pairPoint() && pointPair(v0, v1));
+    // }
 
 
-    //- Do the three given vertices consitute an outside triangle
-    inline friend bool outsideTriangle
-    (
-        const indexedVertex& v0,
-        const indexedVertex& v1,
-        const indexedVertex& v2
-    )
-    {
-        return (v0.farPoint() || v0.ppSlave())
-            || (v1.farPoint() || v1.ppSlave())
-            || (v2.farPoint() || v2.ppSlave());
-    }
+    // //- Do the three given vertices consitute an outside triangle
+    // inline friend bool outsideTriangle
+    // (
+    //     const indexedVertex& v0,
+    //     const indexedVertex& v1,
+    //     const indexedVertex& v2
+    // )
+    // {
+    //     return (v0.farPoint() || v0.ppSlave())
+    //         || (v1.farPoint() || v1.ppSlave())
+    //         || (v2.farPoint() || v2.ppSlave());
+    // }
 
 
     // inline void operator=(const Delaunay::Finite_vertices_iterator vit)