ENH: Add concave cell checking

applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C

@@ -291,6 +291,22 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         }
     }
 
+    if (allGeometry)
+    {
+        cellSet cells(mesh, "concaveCells", mesh.nCells()/100);
+        if (mesh.checkConcaveCells(true, &cells))
+        {
+            noFailedChecks++;
+
+            label nCells = returnReduce(cells.size(), sumOp<label>());
+
+            Info<< "  <<Writing " << nCells
+                << " concave cells to set " << cells.name() << endl;
+            cells.instance() = mesh.pointsInstance();
+            cells.write();
+        }
+    }
+
 
     return noFailedChecks;
 }

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H

@@ -8,10 +8,10 @@
 License
     This file is part of OpenFOAM.
 
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
 
     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
@@ -19,8 +19,7 @@ License
     for more details.
 
     You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
     Foam::primitiveMesh
@@ -292,6 +291,9 @@ class primitiveMesh
             //- Skewness warning threshold
             static scalar skewThreshold_;
 
+            //- Threshold where faces are considered coplanar
+            static scalar planarCosAngle_;
+
 
 protected:
 
@@ -647,6 +649,13 @@ public:
                     labelHashSet* setPtr = NULL
                 ) const;
 
+                //- Check for concave cells
+                bool checkConcaveCells
+                (
+                    const bool report = false,
+                    labelHashSet* setPtr = NULL
+                ) const;
+
 
             //- Check mesh topology for correctness.
             //  Returns false for no error.

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C

@@ -8,10 +8,10 @@
 License
     This file is part of OpenFOAM.
 
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
 
     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
@@ -19,8 +19,7 @@ License
     for more details.
 
     You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 \*---------------------------------------------------------------------------*/
 
@@ -38,6 +37,7 @@ Foam::scalar Foam::primitiveMesh::closedThreshold_  = 1.0e-6;
 Foam::scalar Foam::primitiveMesh::aspectThreshold_  = 1000;
 Foam::scalar Foam::primitiveMesh::nonOrthThreshold_ = 70;    // deg
 Foam::scalar Foam::primitiveMesh::skewThreshold_    = 4;
+Foam::scalar Foam::primitiveMesh::planarCosAngle_   = 1.0e-6;
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -2094,6 +2094,129 @@ bool Foam::primitiveMesh::checkCellDeterminant
 }
 
 
+bool Foam::primitiveMesh::checkConcaveCells
+(
+    const bool report,
+    labelHashSet* setPtr
+) const
+{
+    if (debug)
+    {
+        Info<< "bool primitiveMesh::checkConcaveCells(const bool"
+            << ", labelHashSet*) const: "
+            << "checking for concave cells" << endl;
+    }
+
+    const cellList& c = cells();
+    const labelList& fOwner = faceOwner();
+    const vectorField& fAreas = faceAreas();
+    const pointField& fCentres = faceCentres();
+    const pointField& pts = points();
+
+    label nConcaveCells = 0;
+
+    forAll (c, cellI)
+    {
+        const cell& cFaces = c[cellI];
+
+        const labelList& cPoints = cellPoints()[cellI];
+
+        bool concave = false;
+
+        forAll(cFaces, i)
+        {
+            if (concave)
+            {
+                break;
+            }
+
+            label fI = cFaces[i];
+
+            const point& fC = fCentres[fI];
+
+            const face& f = faces()[fI];
+
+            vector fN = fAreas[fI];
+
+            fN /= max(mag(fN), VSMALL);
+
+            // Flip normal if required so that it is always pointing out of
+            // the cell
+            if (fOwner[fI] != cellI)
+            {
+                fN *= -1;
+            }
+
+            // Is any vertex of the cell on the wrong side of the
+            // plane of this face?
+
+            forAll(cPoints, cPtI)
+            {
+                label ptI = cPoints[cPtI];
+
+                // Skip points that are on this face
+                if (findIndex(f, ptI) > -1)
+                {
+                    continue;
+                }
+
+                const point& pt = pts[ptI];
+
+                // If the cell is concave, the point will be on the
+                // positive normal side of the plane of f, defined by
+                // its centre and normal, and the angle between (pt -
+                // fC) and fN will be less than 90 degrees, so the dot
+                // product will be positive.
+
+                vector pC = (pt - fC);
+
+                pC /= max(mag(pC), VSMALL);
+
+                if ((pC & fN) > -planarCosAngle_)
+                {
+                    // Concave or planar face
+
+                    concave = true;
+
+                    if (setPtr)
+                    {
+                        setPtr->insert(cellI);
+                    }
+
+                    nConcaveCells++;
+
+                    break;
+                }
+            }
+        }
+    }
+
+    reduce(nConcaveCells, sumOp<label>());
+
+    if (nConcaveCells > 0)
+    {
+        if (debug || report)
+        {
+            Info<< " ***Concave cells found, number of cells: "
+                << nConcaveCells << endl;
+        }
+
+        return true;
+    }
+    else
+    {
+        if (debug || report)
+        {
+            Info<< "    Concave cell check OK." << endl;
+        }
+
+        return false;
+    }
+
+    return false;
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::primitiveMesh::checkTopology(const bool report) const