From 1ca0932b9be1abda44278b04d1bc39d46d9cf8d3 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Thu, 7 Jan 2010 17:13:13 +0000
Subject: [PATCH] Stabilising use of detJ for negative or zero values.

Moving removeUnusedPoints to calcDualMesh from
createFacesOwnerNeighbourAndPatches to allow the latter to be used to
build an intermediate mesh.

Adding the use of motionSmoother::checkMesh (as per snappyHexMesh) to
a new function: checkPolyMeshQuality.  Will be the analysis component
of a filtering back-tracking scheme.
---
 .../conformalVoronoiMesh.H                    |   5 +-
 .../conformalVoronoiMeshCalcDualMesh.C        | 117 ++++++++++++++++--
 2 files changed, 108 insertions(+), 14 deletions(-)

diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 94823908b1..5ec0a1bc78 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -513,13 +513,16 @@ private:
             scalar collapseSizeLimitCoeff
         ) const;
 
+        //- Create a polyMesh and check its quality, reports which
+        //  elements damage the mesh quality, allowing backtracking.
+        label checkPolyMeshQuality(const pointField& pts) const;
+
         //- Re-index all of the the Delaunay cells
         void reindexDualVertices(const Map<label>& dualPtIndexMap);
 
         //- Create all of the internal and boundary faces
         void createFacesOwnerNeighbourAndPatches
         (
-            pointField& points,
             faceList& faces,
             labelList& owner,
             labelList& neighbour,
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index f9d4845868..05c8ffbb93 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -25,6 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "conformalVoronoiMesh.H"
+#include "motionSmoother.H"
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
@@ -354,13 +355,16 @@ void Foam::conformalVoronoiMesh::calcDualMesh
 
     collapseFaces(points);
 
+    label nBadQualityFaces = checkPolyMeshQuality(points);
+
+    Info<< "Found " << nBadQualityFaces << " bad quality faces" << endl;
+
     // Final dual face and owner neighbour construction
 
     timeCheck();
 
     createFacesOwnerNeighbourAndPatches
     (
-        points,
         faces,
         owner,
         neighbour,
@@ -369,6 +373,10 @@ void Foam::conformalVoronoiMesh::calcDualMesh
         patchStarts,
         false
     );
+
+    removeUnusedPoints(faces, points);
+
+    timeCheck();
 }
 
 
@@ -730,7 +738,7 @@ Foam::label Foam::conformalVoronoiMesh::smoothSurfaceDualFaces
 {
     label nCollapsedFaces = 0;
 
-    const scalar cosPerpendicularToleranceAngle = cos(degToRad(75));
+    const scalar cosPerpendicularToleranceAngle = cos(degToRad(80));
 
     for
     (
@@ -1288,7 +1296,9 @@ bool Foam::conformalVoronoiMesh::collapseFaceToEdge
 
     scalar aspectRatio = 1;
 
-    scalar detJ = det(J);
+    // Calculating determinant, including stabilisation for zero or
+    // small negative values
+    scalar detJ = max(det(J), SMALL);
 
     if (detJ < 1e-5)
     {
@@ -1321,11 +1331,6 @@ bool Foam::conformalVoronoiMesh::collapseFaceToEdge
 
         // Empirical correlation for high aspect ratio faces
 
-        if (detJ < VSMALL)
-        {
-            Info<< "if (detJ < VSMALL) " << detJ << endl;
-        }
-
         aspectRatio = sqrt(0.35/detJ);
 
         // Info<< "# Longest edge determined collapseAxis" << endl;
@@ -1585,6 +1590,97 @@ bool Foam::conformalVoronoiMesh::collapseFaceToEdge
 }
 
 
+Foam::label Foam::conformalVoronoiMesh::checkPolyMeshQuality
+(
+    const pointField& pts
+) const
+{
+    faceList faces;
+    labelList owner;
+    labelList neighbour;
+    wordList patchNames;
+    labelList patchSizes;
+    labelList patchStarts;
+
+    timeCheck();
+
+    Info<< "    Creating polyMesh to assess quality" << endl;
+
+    createFacesOwnerNeighbourAndPatches
+    (
+        faces,
+        owner,
+        neighbour,
+        patchNames,
+        patchSizes,
+        patchStarts,
+        false
+    );
+
+    IOobject io
+    (
+        "cvMesh_temporary",
+        runTime_.constant(),
+        runTime_,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    );
+
+    polyMesh pMesh
+    (
+        io,
+        xferCopy(pts),
+        xferMove(faces),
+        xferMove(owner),
+        xferMove(neighbour)
+    );
+
+    List<polyPatch*> patches(patchStarts.size());
+
+    forAll (patches, p)
+    {
+        patches[p] = new polyPatch
+        (
+            patchNames[p],
+            patchSizes[p],
+            patchStarts[p],
+            p,
+            pMesh.boundaryMesh()
+        );
+    }
+
+    pMesh.addPatches(patches);
+
+    timeCheck();
+
+    labelHashSet wrongFaces(pMesh.nFaces()/100);
+
+    motionSmoother::checkMesh
+    (
+        false,
+        pMesh,
+        cvMeshControls().cvMeshDict().subDict("meshQualityControls"),
+        wrongFaces
+    );
+
+    forAllConstIter(labelHashSet, wrongFaces, iter)
+    {
+        label faceI = iter.key();
+
+        Info<< faceI << " " << pMesh.faces()[faceI] << endl;
+    }
+
+    return wrongFaces.size();
+
+    // For parallel running:
+    // return returnReduce
+    // (
+    //     wrongFaces.size(),
+    //     sumOp<label>()
+    // );
+}
+
+
 void Foam::conformalVoronoiMesh::reindexDualVertices
 (
     const Map<label>& dualPtIndexMap
@@ -1607,7 +1703,6 @@ void Foam::conformalVoronoiMesh::reindexDualVertices
 
 void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
 (
-    pointField& points,
     faceList& faces,
     labelList& owner,
     labelList& neighbour,
@@ -1779,10 +1874,6 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
         patchOwners,
         false
     );
-
-    removeUnusedPoints(faces, points);
-
-    timeCheck();
 }