From cf0c71789ae9f7e53f6a9ab89cd6c685092ffaac Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Wed, 13 Jan 2010 16:12:05 +0000
Subject: [PATCH] Adding filter limit scalars to the Delaunay cells.  After
 initial face check, identifying faces that have bad quality, identifying the
 points that support the faces, and 'limiting' them.

Wrong faces are also those that are attached to a cell with 1-3 faces.

At each iteration of filtering, the Delaunay cells are all reindexed.
Now a separate function.

Currently 'limiting' means: do not collapse any face that has any
vertex that has been limited at all and do not move a vertex to the
boundary during the surface smooth if it has been limited.  This is
the most agressive form of limiting and leaves faces in the mesh that
might be safely removed with an less aggressive filter.

Limiting of all points connected to a cell which has a limited
point could be used (comment out code).

Removed some redundant code from indexedCell.
---
 .../conformalVoronoiMesh.C                    |  40 +++
 .../conformalVoronoiMesh.H                    |  12 +
 .../conformalVoronoiMeshCalcDualMesh.C        | 271 ++++++++++++++----
 .../conformalVoronoiMesh/indexedCell.H        |  60 ++--
 .../conformalVoronoiMesh/indexedVertex.H      |  28 +-
 src/mesh/conformalVoronoiMesh/log             |  20 --
 6 files changed, 321 insertions(+), 110 deletions(-)
 delete mode 100644 src/mesh/conformalVoronoiMesh/log

diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index fd1e5c4d42..e146e488a0 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -1030,6 +1030,46 @@ Foam::face Foam::conformalVoronoiMesh::buildDualFace
 }
 
 
+Foam::scalar Foam::conformalVoronoiMesh::minFilterLimit
+(
+    const Triangulation::Finite_edges_iterator& eit
+) const
+{
+    Cell_circulator ccStart = incident_cells(*eit);
+    Cell_circulator cc = ccStart;
+
+    scalar minFilterLimit = GREAT;
+
+    do
+    {
+        if (cc->cellIndex() < 0)
+        {
+            Cell_handle c = eit->first;
+            Vertex_handle vA = c->vertex(eit->second);
+            Vertex_handle vB = c->vertex(eit->third);
+
+            FatalErrorIn("Foam::conformalVoronoiMesh::buildDualFace")
+                << "Dual face uses circumcenter defined by a "
+                << "Delaunay tetrahedron with no internal "
+                << "or boundary points.  Defining Delaunay edge ends: "
+                << topoint(vA->point()) << " "
+                << topoint(vB->point()) << nl
+                << exit(FatalError);
+        }
+
+        if (cc->filterLimit() < minFilterLimit)
+        {
+            minFilterLimit = cc->filterLimit();
+        }
+
+        cc++;
+
+    } while (cc != ccStart);
+
+    return minFilterLimit;
+}
+
+
 bool Foam::conformalVoronoiMesh::ownerAndNeighbour
 (
     Vertex_handle vA,
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 9e9c45fee7..095ddc5a18 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -357,6 +357,14 @@ private:
             const Triangulation::Finite_edges_iterator& eit
         ) const;
 
+        //- Finds the minimum filterLimit of the dual vertices
+        //  (Delaunay cells) that form the dual face produced by the
+        //  supplied edge
+        scalar minFilterLimit
+        (
+            const Triangulation::Finite_edges_iterator& eit
+        ) const;
+
         //- Determines the owner and neighbour labels for dual cells
         //  corresponding to the dual face formed by the supplied
         //  Delaunay vertices.  If the dual face is a boundary face
@@ -552,6 +560,10 @@ private:
         //  elements damage the mesh quality, allowing backtracking.
         label checkPolyMeshQuality(const pointField& pts) const;
 
+        //- Index all of the the Delaunay cells and calculate their
+        //- dual points
+        void indexDualVertices(pointField& pts);
+
         //- Re-index all of the the Delaunay cells
         void reindexDualVertices(const Map<label>& dualPtIndexMap);
 
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index b0f9adbf78..489401a77d 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -76,53 +76,28 @@ void Foam::conformalVoronoiMesh::calcDualMesh
         }
     }
 
-    // Indexing Delaunay cells, which are the dual vertices
+    indexDualVertices(points);
 
-    label dualVertI = 0;
-
-    points.setSize(number_of_cells());
-
-    for
-    (
-        Triangulation::Finite_cells_iterator cit = finite_cells_begin();
-        cit != finite_cells_end();
-        ++cit
-    )
     {
-        if
-        (
-            cit->vertex(0)->internalOrBoundaryPoint()
-         || cit->vertex(1)->internalOrBoundaryPoint()
-         || cit->vertex(2)->internalOrBoundaryPoint()
-         || cit->vertex(3)->internalOrBoundaryPoint()
-        )
-        {
-            cit->cellIndex() = dualVertI;
-            points[dualVertI] = topoint(dual(cit));
-            dualVertI++;
-        }
-        else
-        {
-            cit->cellIndex() = -1;
-        }
+        // No-risk face filtering to get rid of zero area faces and
+        // establish if the mesh can be produced at all to the
+        // specified criteria
+
+        Info<< nl << "    Merging close points" << endl;
+
+        // There is no guarantee that a merge of close points is
+        // no-risk, but it seems to work using 1e-4 as the mergeClosenessCoeff
+        mergeCloseDualVertices(points);
     }
 
-    points.setSize(dualVertI);
-
-    // No-risk face filtering to get rid of zero area faces and
-    // establish if the mesh can be produced at all to the
-    // specified criteria
-
-    Info<< nl << "    Merging close points" << endl;
-
-    mergeCloseDualVertices(points);
-
     label nInitialBadQualityFaces = checkPolyMeshQuality(points);
 
     Info<< "Initial check before face collapse, found "
         << nInitialBadQualityFaces << " bad quality faces"
         << endl;
 
+    HashSet<labelPair, labelPair::Hash<> > deferredCollapseFaces;
+
     if (nInitialBadQualityFaces > 0)
     {
         Info<< "A mesh could not be produced to satisfy the specified quality "
@@ -132,24 +107,41 @@ void Foam::conformalVoronoiMesh::calcDualMesh
             << "be applied in areas where problems are occurring."
             << endl;
     }
-
-    HashSet<labelPair, labelPair::Hash<> > deferredCollapseFaces;
-
+    else
     {
-        // Risky and undo-able face filtering to reduce the face count
-        // as much as possible staying within the specified criteria
+        label nBadQualityFaces = 0;
 
-        Info<< nl << "    Smoothing surface" << endl;
+        do
+        {
+            // Reindexing the Delaunay cells and regenerating the
+            // points resets the mesh to the starting condition.
 
-        smoothSurface(points);
+            indexDualVertices(points);
 
-        Info<< nl << "    Collapsing unnecessary faces" << endl;
+            {
+                Info<< nl << "    Merging close points" << endl;
 
-        collapseFaces(points, deferredCollapseFaces);
+                mergeCloseDualVertices(points);
+            }
 
-        label nBadQualityFaces = checkPolyMeshQuality(points);
+            {
+                // Risky and undo-able face filtering to reduce the face count
+                // as much as possible, staying within the specified criteria
 
-        Info<< "Found " << nBadQualityFaces << " bad quality faces" << endl;
+                Info<< nl << "    Smoothing surface" << endl;
+
+                smoothSurface(points);
+
+                Info<< nl << "    Collapsing unnecessary faces" << endl;
+
+                collapseFaces(points, deferredCollapseFaces);
+            }
+
+            nBadQualityFaces = checkPolyMeshQuality(points);
+
+            Info<< "Found " << nBadQualityFaces << " bad quality faces" << endl;
+
+        } while (nBadQualityFaces > 0);
     }
 
     // Final dual face and owner neighbour construction
@@ -485,6 +477,12 @@ void Foam::conformalVoronoiMesh::smoothSurface(pointField& pts)
     {
         label ptI = cit->cellIndex();
 
+        if (cit->filterLimit() < 1.0)
+        {
+            // This vertex has been limited, skip
+            continue;
+        }
+
         // Only cells with indices > -1 are valid
         if (ptI > -1)
         {
@@ -580,6 +578,14 @@ Foam::label Foam::conformalVoronoiMesh::smoothSurfaceDualFaces
                 continue;
             }
 
+            scalar minFL = minFilterLimit(eit);
+
+            if (minFL < 1.0)
+            {
+                // A vertex on this face has been limited, skip
+                continue;
+            }
+
             pointIndexHit surfHit;
             label hitSurface;
 
@@ -749,6 +755,14 @@ Foam::label Foam::conformalVoronoiMesh::collapseFaces
                 continue;
             }
 
+            scalar minFL = minFilterLimit(eit);
+
+            if (minFL < 1.0)
+            {
+                // A vertex on this face has been limited, skip
+                continue;
+            }
+
             scalar targetFaceSize = averageAnyCellSize(vA, vB);
 
             faceCollapseMode mode = collapseFace
@@ -1256,31 +1270,138 @@ Foam::label Foam::conformalVoronoiMesh::checkPolyMeshQuality
 
     motionSmoother::checkMesh
     (
-        false,
+        false,  // report
         pMesh,
         cvMeshControls().cvMeshDict().subDict("meshQualityControls"),
         wrongFaces
     );
 
+    label nInvalidPolyhedra = 0;
+
     const cellList& cells = pMesh.cells();
 
     forAll(cells, cI)
     {
-        if (cells[cI].size() < 4)
+        if (cells[cI].size() < 4 && cells[cI].size() > 0)
         {
-            Info<< "cell " << cI
-                << " has " << cells[cI].size() << " faces."
-                << endl;
+            // Info<< "cell " << cI << " " << cells[cI]
+            //     << " has " << cells[cI].size() << " faces."
+            //     << endl;
+
+            nInvalidPolyhedra++;
+
+            forAll(cells[cI], cFI)
+            {
+                wrongFaces.insert(cells[cI][cFI]);
+            }
         }
     }
 
+    Info<< "Cells with fewer than 4 faces                              : "
+        << nInvalidPolyhedra << endl;
+
+    PackedBoolList ptToBeLimited(pts.size(), false);
+
+    forAllConstIter(labelHashSet, wrongFaces, iter)
+    {
+        const face f = pMesh.faces()[iter.key()];
+
+        forAll(f, fPtI)
+        {
+            ptToBeLimited[f[fPtI]] = true;
+        }
+    }
+
+    // // Limit connected cells
+
+    // labelHashSet limitCells(pMesh.nCells()/100);
+
+    // const labelListList& ptCells = pMesh.pointCells();
+
     // forAllConstIter(labelHashSet, wrongFaces, iter)
     // {
-    //     label faceI = iter.key();
+    //     const face f = pMesh.faces()[iter.key()];
 
-    //     Info<< faceI << " " << pMesh.faces()[faceI] << endl;
+    //     forAll(f, fPtI)
+    //     {
+    //         label ptI = f[fPtI];
+
+    //         const labelList& pC = ptCells[ptI];
+
+    //         forAll(pC, pCI)
+    //         {
+    //             limitCells.insert(pC[pCI]);
+    //         }
+    //     }
     // }
 
+    // const labelListList& cellPts = pMesh.cellPoints();
+
+    // forAllConstIter(labelHashSet, limitCells, iter)
+    // {
+    //     label cellI = iter.key();
+
+    //     const labelList& cP = cellPts[cellI];
+
+    //     forAll(cP, cPI)
+    //     {
+    //         ptToBeLimited[cP[cPI]] = true;
+    //     }
+    // }
+
+
+    // Limit Delaunay cell filter values
+
+    for
+    (
+        Triangulation::Finite_cells_iterator cit = finite_cells_begin();
+        cit != finite_cells_end();
+        ++cit
+    )
+    {
+        label cI = cit->cellIndex();
+
+        if (cI >= 0)
+        {
+            if (ptToBeLimited[cI] == true)
+            {
+                cit->filterLimit() *= 0.75;
+            }
+        }
+    }
+
+    // Determine the minimum minFilterLimit
+
+    scalar minMinFilterLimit = GREAT;
+
+    for
+    (
+        Triangulation::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->internalOrBoundaryPoint()
+         || vB->internalOrBoundaryPoint()
+        )
+        {
+            scalar minFL = minFilterLimit(eit);
+
+            if (minFL < minMinFilterLimit)
+            {
+                minMinFilterLimit = minFL;
+            }
+        }
+    }
+
+    Info<< "minMinFilterLimit " << minMinFilterLimit << endl;
+
     return wrongFaces.size();
 
     // For parallel running:
@@ -1292,6 +1413,46 @@ Foam::label Foam::conformalVoronoiMesh::checkPolyMeshQuality
 }
 
 
+void Foam::conformalVoronoiMesh::indexDualVertices
+(
+    pointField& pts
+)
+{
+    // Indexing Delaunay cells, which are the dual vertices
+
+    label dualVertI = 0;
+
+    pts.setSize(number_of_cells());
+
+    for
+    (
+        Triangulation::Finite_cells_iterator cit = finite_cells_begin();
+        cit != finite_cells_end();
+        ++cit
+    )
+    {
+        if
+        (
+            cit->vertex(0)->internalOrBoundaryPoint()
+         || cit->vertex(1)->internalOrBoundaryPoint()
+         || cit->vertex(2)->internalOrBoundaryPoint()
+         || cit->vertex(3)->internalOrBoundaryPoint()
+        )
+        {
+            cit->cellIndex() = dualVertI;
+            pts[dualVertI] = topoint(dual(cit));
+            dualVertI++;
+        }
+        else
+        {
+            cit->cellIndex() = -1;
+        }
+    }
+
+    pts.setSize(dualVertI);
+}
+
+
 void Foam::conformalVoronoiMesh::reindexDualVertices
 (
     const Map<label>& dualPtIndexMap
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
index 54775ed4a7..91ac79c57a 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
@@ -36,6 +36,7 @@ Description
 
 #include <CGAL/Triangulation_3.h>
 #include "indexedVertex.H"
+#include "scalar.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,7 +44,7 @@ namespace CGAL
 {
 
 /*---------------------------------------------------------------------------*\
-                       Class indexedCell Declaration
+                          Class indexedCell Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Gt, class Cb=CGAL::Triangulation_cell_base_3<Gt> >
@@ -54,20 +55,15 @@ class indexedCell
     // Private data
 
         //- The index for this tetrahedral cell
-        //  -1: tetrahedral and changed and associated data needs updating
-        // >=0: index of cells, set by external update algorithm
         int index_;
 
+        //- The applicable fraction of the face filtering mechanism
+        //  controls
+        Foam::scalar filterLimit_;
+
 
 public:
 
-    enum cellTypes
-    {
-        ctUnchanged           = 0,
-        ctChanged             = -1,
-        ctSaveChanged         = -2
-    };
-
     typedef typename Cb::Vertex_handle                   Vertex_handle;
     typedef typename Cb::Cell_handle                     Cell_handle;
 
@@ -82,18 +78,22 @@ public:
     indexedCell()
     :
         Cb(),
-        index_(ctChanged)
+        index_(-1),
+        filterLimit_(1.0)
     {}
 
+
     indexedCell
     (
         Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3
     )
     :
         Cb(v0, v1, v2, v3),
-        index_(ctChanged)
+        index_(-1),
+        filterLimit_(1.0)
     {}
 
+
     indexedCell
     (
         Vertex_handle v0,
@@ -107,39 +107,33 @@ public:
     )
     :
         Cb(v0, v1, v2, v3, n0, n1, n2, n3),
-        index_(ctChanged)
+        index_(-1),
+        filterLimit_(1.0)
     {}
 
-    void set_vertex(int i, Vertex_handle v)
-    {
-        index_ = ctChanged;
-        Cb::set_vertex(i, v);
-    }
-
-    void set_vertices()
-    {
-        index_ = ctChanged;
-        Cb::set_vertices();
-    }
-
-    void set_vertices
-    (
-        Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3
-    )
-    {
-        index_ = ctChanged;
-        Cb::set_vertices(v0, v1, v2, v3);
-    }
 
     int& cellIndex()
     {
         return index_;
     }
 
+
     int cellIndex() const
     {
         return index_;
     }
+
+
+    inline Foam::scalar& filterLimit()
+    {
+        return filterLimit_;
+    }
+
+
+    inline Foam::scalar filterLimit() const
+    {
+        return filterLimit_;
+    }
 };
 
 
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
index 93fcb12346..a9116054cd 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
@@ -35,7 +35,6 @@ Description
 #define indexedVertex_H
 
 #include <CGAL/Triangulation_3.h>
-#include "List.H"
 #include "tensor.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -44,7 +43,7 @@ namespace CGAL
 {
 
 /*---------------------------------------------------------------------------*\
-                       Class indexedVertex Declaration
+                         Class indexedVertex Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Gt, class Vb=CGAL::Triangulation_vertex_base_3<Gt> >
@@ -91,6 +90,7 @@ public:
         typedef indexedVertex<Gt,Vb2>                           Other;
     };
 
+
     indexedVertex()
     :
         Vb(),
@@ -100,6 +100,7 @@ public:
         targetCellSize_(0.0)
     {}
 
+
     indexedVertex(const Point& p)
     :
         Vb(p),
@@ -109,6 +110,7 @@ public:
         targetCellSize_(0.0)
     {}
 
+
     indexedVertex(const Point& p, Cell_handle f)
     :
         Vb(f, p),
@@ -118,6 +120,7 @@ public:
         targetCellSize_(0.0)
     {}
 
+
     indexedVertex(Cell_handle f)
     :
         Vb(f),
@@ -133,70 +136,83 @@ public:
         return index_;
     }
 
+
     int index() const
     {
         return index_;
     }
 
+
     int& type()
     {
         return type_;
     }
 
+
     int type() const
     {
         return type_;
     }
 
+
     inline Foam::tensor& alignment()
     {
         return alignment_;
     }
 
+
     inline const Foam::tensor& alignment() const
     {
         return alignment_;
     }
 
+
     inline Foam::scalar& targetCellSize()
     {
         return targetCellSize_;
     }
 
+
     inline Foam::scalar targetCellSize() const
     {
         return targetCellSize_;
     }
 
+
     inline bool uninitialised() const
     {
         return type_ == ptInternalPoint && index_ == ptInternalPoint;
     }
 
+
     //- Is point a far-point
     inline bool farPoint() const
     {
         return type_ == ptFarPoint;
     }
 
+
     //- Is point internal, i.e. not on boundary
     inline bool internalPoint() const
     {
         return type_ <= ptInternalPoint;
     }
 
+
     //- Set the point to be internal
     inline void setInternal()
     {
         type_ = ptInternalPoint;
     }
 
+
     //- Is point internal and near the boundary
     inline bool nearBoundary() const
     {
         return type_ == ptNearBoundaryPoint;
     }
 
+
     //- Set the point to be near the boundary
     inline void setNearBoundary()
     {
@@ -209,12 +225,14 @@ public:
         return type_ == ptMirrorPoint;
     }
 
+
     //- Either master or slave of pointPair.
     inline bool pairPoint() const
     {
         return type_ >= 0;
     }
 
+
     //- Master of a pointPair is the lowest numbered one.
     inline bool ppMaster() const
     {
@@ -228,6 +246,7 @@ public:
         }
     }
 
+
     //- Slave of a pointPair is the highest numbered one.
     inline bool ppSlave() const
     {
@@ -241,18 +260,21 @@ public:
         }
     }
 
+
     //- Either original internal point or master of pointPair.
     inline bool internalOrBoundaryPoint() const
     {
         return internalPoint() || ppMaster();
     }
 
+
     //- Is point near the boundary or part of the boundary definition
     inline bool nearOrOnBoundary() const
     {
         return pairPoint() || mirrorPoint() || nearBoundary();
     }
 
+
     //- Do the two given vertices consitute a boundary point-pair
     inline friend bool pointPair
     (
@@ -263,6 +285,7 @@ public:
         return v0.index_ == v1.type_ || v1.index_ == v0.type_;
     }
 
+
     //- Do the three given vertices consitute a boundary triangle
     inline friend bool boundaryTriangle
     (
@@ -276,6 +299,7 @@ public:
             || (v2.pairPoint() && pointPair(v0, v1));
     }
 
+
     //- Do the three given vertices consitute an outside triangle
     inline friend bool outsideTriangle
     (
diff --git a/src/mesh/conformalVoronoiMesh/log b/src/mesh/conformalVoronoiMesh/log
deleted file mode 100644
index 0e31b724c8..0000000000
--- a/src/mesh/conformalVoronoiMesh/log
+++ /dev/null
@@ -1,20 +0,0 @@
-/*---------------------------------------------------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-Build  : dev-79ee164b7a49
-Exec   : cvMesh
-Date   : Jan 12 2010
-Time   : 19:49:15
-Host   : rockall
-PID    : 14798
-Case   : /home/graham/OpenFOAM/OpenFOAM-dev/src/mesh/conformalVoronoiMesh
-nProcs : 1
-SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-Create time
-