From 449b7d77315f4539ec6ab3b1962dacdaa3bce88b Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 22 Nov 2012 09:11:39 +0000
Subject: [PATCH] ENH: polyMeshGeometry: added flatness check

---
 .../motionSmoother/motionSmootherCheck.C      |  51 ++++++-
 .../polyMeshGeometry/polyMeshGeometry.C       | 137 +++++++++++++++++-
 .../polyMeshGeometry/polyMeshGeometry.H       |  25 +++-
 3 files changed, 206 insertions(+), 7 deletions(-)

diff --git a/src/dynamicMesh/motionSmoother/motionSmootherCheck.C b/src/dynamicMesh/motionSmoother/motionSmootherCheck.C
index da20ea741e..77ce89ca71 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherCheck.C
+++ b/src/dynamicMesh/motionSmoother/motionSmootherCheck.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,6 +104,8 @@ bool Foam::motionSmoother::checkMesh
     (
         readScalar(dict.lookup("minTriangleTwist", true))
     );
+    scalar minFaceFlatness = -1.0;
+    dict.readIfPresent("minFaceFlatness", minFaceFlatness, true);
     const scalar minDet
     (
         readScalar(dict.lookup("minDeterminant", true))
@@ -356,6 +358,30 @@ bool Foam::motionSmoother::checkMesh
         nWrongFaces = nNewWrongFaces;
     }
 
+    if (minFaceFlatness > -SMALL)
+    {
+        polyMeshGeometry::checkFaceFlatness
+        (
+            report,
+            minFaceFlatness,
+            mesh,
+            mesh.faceAreas(),
+            mesh.faceCentres(),
+            mesh.points(),
+            checkFaces,
+            &wrongFaces
+        );
+
+        label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
+
+        Info<< "    faces with flatness < "
+            << setw(5) << minFaceFlatness
+            << "                      : "
+            << nNewWrongFaces-nWrongFaces << endl;
+
+        nWrongFaces = nNewWrongFaces;
+    }
+
     if (minDet > -1)
     {
         polyMeshGeometry::checkCellDeterminant
@@ -479,6 +505,8 @@ bool Foam::motionSmoother::checkMesh
     (
         readScalar(dict.lookup("minTriangleTwist", true))
     );
+    scalar minFaceFlatness = -1.0;
+    dict.readIfPresent("minFaceFlatness", minFaceFlatness, true);
     const scalar minDet
     (
         readScalar(dict.lookup("minDeterminant", true))
@@ -697,6 +725,27 @@ bool Foam::motionSmoother::checkMesh
         nWrongFaces = nNewWrongFaces;
     }
 
+    if (minFaceFlatness > -1)
+    {
+        meshGeom.checkFaceFlatness
+        (
+            report,
+            minFaceFlatness,
+            meshGeom.mesh().points(),
+            checkFaces,
+            &wrongFaces
+        );
+
+        label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
+
+        Info<< "    faces with flatness < "
+            << setw(5) << minFaceFlatness
+            << "                      : "
+            << nNewWrongFaces-nWrongFaces << endl;
+
+        nWrongFaces = nNewWrongFaces;
+    }
+
     if (minDet > -1)
     {
         meshGeom.checkCellDeterminant
diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
index 1670de1c92..13786a1c1c 100644
--- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
+++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
@@ -341,8 +341,8 @@ bool Foam::polyMeshGeometry::checkFaceTet
             {
                 Pout<< "bool polyMeshGeometry::checkFaceTets("
                     << "const bool, const scalar, const pointField&"
-                    << ", const pointField&, const labelList&,"
-                    << " labelHashSet*): "
+                    << ", const pointField&"
+                    << ", const labelList&, labelHashSet*) : "
                     << "face " << faceI
                     << " has a triangle that points the wrong way."
                      << endl
@@ -373,9 +373,6 @@ Foam::polyMeshGeometry::polyMeshGeometry(const polyMesh& mesh)
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 //- Take over properties from mesh
@@ -2014,6 +2011,113 @@ bool Foam::polyMeshGeometry::checkTriangleTwist
 }
 
 
+bool Foam::polyMeshGeometry::checkFaceFlatness
+(
+    const bool report,
+    const scalar minFlatness,
+    const polyMesh& mesh,
+    const vectorField& faceAreas,
+    const vectorField& faceCentres,
+    const pointField& p,
+    const labelList& checkFaces,
+    labelHashSet* setPtr
+)
+{
+    if (minFlatness < -SMALL || minFlatness > 1+SMALL)
+    {
+        FatalErrorIn
+        (
+            "polyMeshGeometry::checkFaceFlatness"
+            "(const bool, const scalar, const polyMesh&, const pointField&"
+            ", const pointField&, const labelList&, labelHashSet*)"
+        )   << "minFlatness should be [0..1] but is now " << minFlatness
+            << abort(FatalError);
+    }
+
+    const faceList& fcs = mesh.faces();
+
+    label nWarped = 0;
+
+    forAll(checkFaces, i)
+    {
+        label faceI = checkFaces[i];
+
+        const face& f = fcs[faceI];
+
+        if (f.size() > 3)
+        {
+            const point& fc = faceCentres[faceI];
+
+            // Sum triangle areas
+            scalar sumArea = 0.0;
+
+            forAll(f, fp)
+            {
+                sumArea += triPointRef
+                (
+                    p[f[fp]],
+                    p[f.nextLabel(fp)],
+                    fc
+                ).mag();
+            }
+
+            if (sumArea/mag(faceAreas[faceI]) < minFlatness)
+            {
+                nWarped++;
+
+                if (setPtr)
+                {
+                    setPtr->insert(faceI);
+                }
+            }
+        }
+    }
+
+    reduce(nWarped, sumOp<label>());
+
+    if (report)
+    {
+        if (nWarped> 0)
+        {
+            Info<< "There are " << nWarped
+                << " faces with area of invidual triangles"
+                << " compared to overall area less than "
+                << minFlatness << nl << endl;
+        }
+        else
+        {
+            Info<< "All faces are flat in that the area of invidual triangles"
+                << " compared to overall area is less than "
+                << minFlatness << nl << endl;
+        }
+    }
+
+    if (nWarped > 0)
+    {
+        if (report)
+        {
+            WarningIn
+            (
+                "polyMeshGeometry::checkFaceFlatness"
+                "(const bool, const scalar, const polyMesh&"
+                ", const pointField&, const pointField&, const labelList&"
+                ", labelHashSet*)"
+            )   << nWarped  << " non-flat faces "
+                << "(area of invidual triangles"
+                << " compared to overall area"
+                << " < " << minFlatness << ") found.\n"
+                << endl;
+        }
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 bool Foam::polyMeshGeometry::checkFaceArea
 (
     const bool report,
@@ -2398,6 +2502,29 @@ bool Foam::polyMeshGeometry::checkTriangleTwist
 }
 
 
+bool Foam::polyMeshGeometry::checkFaceFlatness
+(
+    const bool report,
+    const scalar minFlatness,
+    const pointField& p,
+    const labelList& checkFaces,
+    labelHashSet* setPtr
+) const
+{
+    return checkFaceFlatness
+    (
+        report,
+        minFlatness,
+        mesh_,
+        faceAreas_,
+        faceCentres_,
+        p,
+        checkFaces,
+        setPtr
+    );
+}
+
+
 bool Foam::polyMeshGeometry::checkFaceArea
 (
     const bool report,
diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.H b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.H
index 149fca822a..2ef8604cac 100644
--- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.H
+++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -303,6 +303,19 @@ public:
                 labelHashSet* setPtr
             );
 
+            //- Area of faces v.s. sum of triangle areas
+            static bool checkFaceFlatness
+            (
+                const bool report,
+                const scalar minFlatness,
+                const polyMesh&,
+                const vectorField& faceAreas,
+                const vectorField& faceCentres,
+                const pointField& p,
+                const labelList& checkFaces,
+                labelHashSet* setPtr
+            );
+
             //- Small faces
             static bool checkFaceArea
             (
@@ -314,6 +327,7 @@ public:
                 labelHashSet* setPtr
             );
 
+            //- Area of internal faces v.s. boundary faces
             static bool checkCellDeterminant
             (
                 const bool report,
@@ -413,6 +427,15 @@ public:
                 labelHashSet* setPtr
             ) const;
 
+            bool checkFaceFlatness
+            (
+                const bool report,
+                const scalar minFlatness,
+                const pointField& p,
+                const labelList& checkFaces,
+                labelHashSet* setPtr
+            ) const;
+
             bool checkFaceArea
             (
                 const bool report,