From 4db98cd236b60df0415af133be36c5e913d69577 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Mon, 30 Nov 2009 14:38:39 +0000
Subject: [PATCH] Adding a simple surface smoothing function.  Needs
 functionality added to properly collapse small boundary normal faces.

Also adding a function to determine if a Delaunay edge corresponds to
a boundary dual face.
---
 .../conformalVoronoiMesh.C                    | 94 ++++++++++++++++++-
 .../conformalVoronoiMesh.H                    | 16 ++++
 .../conformalVoronoiMeshI.H                   | 25 +++++
 3 files changed, 134 insertions(+), 1 deletion(-)

diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 223599a41c..9d2a73349c 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -2079,6 +2079,25 @@ void Foam::conformalVoronoiMesh::calcDualMesh
 
     } while (nPtsMerged > 0);
 
+    // Smooth the surface of the mesh
+
+    Info<< nl << "    Smoothing surface" << endl;
+
+    label nSmoothedVertices = 0;
+
+    do
+    {
+        Map<label> dualPtIndexMap;
+
+        nSmoothedVertices = smoothSurfaceDualFaces(points, dualPtIndexMap);
+
+        Info<< "        Smoothed " << nPtsMerged << " points (0 HARD CODED)"
+            << endl;
+
+        reindexDualVertices(dualPtIndexMap);
+
+    } while (nSmoothedVertices > 0);
+
     // Assess faces for collapse
 
     Info<< nl << "    Collapsing unnecessary faces" << endl;
@@ -2415,6 +2434,79 @@ Foam::label Foam::conformalVoronoiMesh::mergeCloseDualVertices
     return nPtsMerged;
 }
 
+Foam::label Foam::conformalVoronoiMesh::smoothSurfaceDualFaces
+(
+    pointField& pts,
+    Map<label>& dualPtIndexMap
+)
+{
+    label nSmoothedVertices = 0;
+
+    for
+    (
+        Triangulation::Finite_edges_iterator eit = finite_edges_begin();
+        eit != finite_edges_end();
+        ++eit
+    )
+    {
+        Cell_circulator ccStart = incident_cells(*eit);
+        Cell_circulator cc = ccStart;
+
+        do
+        {
+            if (dualPtIndexMap.found(cc->cellIndex()))
+            {
+                // One of the points of this face has already been
+                // collapsed this sweep, leave for next sweep
+                continue;
+            }
+
+        } while (++cc != ccStart);
+
+        // Cell_handle c = eit->first;
+        // Vertex_handle vA = c->vertex(eit->second);
+        // Vertex_handle vB = c->vertex(eit->third);
+
+        if (isBoundaryDualFace(eit))
+        {
+            face dualFace = buildDualFace(eit);
+
+            if (dualFace.size() < 3)
+            {
+                // This face has been collapsed already
+                continue;
+            }
+
+            forAll(dualFace, fPtI)
+            {
+                label ptI = dualFace[fPtI];
+
+                point& pt = pts[ptI];
+
+                pointIndexHit surfHit;
+                label hitSurface;
+
+                geometryToConformTo_.findSurfaceNearest
+                (
+                    pt,
+                    cvMeshControls().spanSqr(),
+                    surfHit,
+                    hitSurface
+                );
+
+                if (surfHit.hit())
+                {
+                    pt = surfHit.hitPoint();
+
+                    // dualPtIndexMap.insert(ptI, ptI);
+                }
+            }
+        }
+    }
+
+    return nSmoothedVertices;
+}
+
 
 Foam::label Foam::conformalVoronoiMesh::collapseFaces
 (
@@ -2724,7 +2816,7 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
         if
         (
             vA->internalOrBoundaryPoint()
-            || vB->internalOrBoundaryPoint()
+         || vB->internalOrBoundaryPoint()
         )
         {
             face newDualFace = buildDualFace(eit);
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 6ed61a115d..2a3b0369b5 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -456,6 +456,13 @@ private:
             const Triangulation::Finite_edges_iterator& eit
         ) const;
 
+        //- Determines if the dual face constructed by the Delaunay
+        //  edge is a boundary face
+        inline bool isBoundaryDualFace
+        (
+            const Triangulation::Finite_edges_iterator& eit
+        ) const;
+
         //- Merge vertices that are very close together
         label mergeCloseDualVertices
         (
@@ -463,6 +470,15 @@ private:
             Map<label>& dualPtIndexMap
         );
 
+        //- Smooth the dual vertices of the dual faces on the boundary
+        //  so that they conform to the surface and remove any
+        //  small, normal oriented faces
+        label smoothSurfaceDualFaces
+        (
+            pointField& pts,
+            Map<label>& dualPtIndexMap
+        );
+
         //- Collapse a face to a point or an edge, modifying and
         //  mapping the points, returns the true if the face was
         //  collapsed in this operation
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
index 81cd6aae50..7b7ccf1929 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
@@ -337,6 +337,31 @@ inline void Foam::conformalVoronoiMesh::insertVb(const Vb& v, label offset)
 }
 
 
+inline bool Foam::conformalVoronoiMesh::isBoundaryDualFace
+(
+    const Triangulation::Finite_edges_iterator& eit
+) const
+{
+    Cell_handle c = eit->first;
+    Vertex_handle vA = c->vertex(eit->second);
+    Vertex_handle vB = c->vertex(eit->third);
+
+    // A dual face on the boundary will result from one Dv inside and
+    // one outside
+    return
+    (
+        (
+            vA->internalOrBoundaryPoint()
+         || vB->internalOrBoundaryPoint()
+        )
+     && (
+            !vA->internalOrBoundaryPoint()
+         || !vB->internalOrBoundaryPoint()
+        )
+    );
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 #ifdef CGAL_INEXACT