From cf5ae6d629843f86b18d900b4a53fc7c62ad5ccd Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 31 Mar 2010 15:53:07 +0100
Subject: [PATCH] ENH: Use concavity check based on face planes.

A cell is concave if there is a face whose face-centre is outside the plane of any other face on the cell.
Another option is to check for any vertices being         outside      ,,    ,,
but this is not consistent with current tracking.
---
 .../primitiveMeshCheck/primitiveMeshCheck.C   | 64 ++++++++-----------
 1 file changed, 28 insertions(+), 36 deletions(-)

diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
index 30b051d1e9..f28e9affd0 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
@@ -2111,7 +2111,6 @@ bool Foam::primitiveMesh::checkConcaveCells
     const labelList& fOwner = faceOwner();
     const vectorField& fAreas = faceAreas();
     const pointField& fCentres = faceCentres();
-    const pointField& pts = points();
 
     label nConcaveCells = 0;
 
@@ -2119,8 +2118,6 @@ bool Foam::primitiveMesh::checkConcaveCells
     {
         const cell& cFaces = c[cellI];
 
-        const labelList& cPoints = cellPoints()[cellI];
-
         bool concave = false;
 
         forAll(cFaces, i)
@@ -2134,8 +2131,6 @@ bool Foam::primitiveMesh::checkConcaveCells
 
             const point& fC = fCentres[fI];
 
-            const face& f = faces()[fI];
-
             vector fN = fAreas[fI];
 
             fN /= max(mag(fN), VSMALL);
@@ -2147,45 +2142,42 @@ bool Foam::primitiveMesh::checkConcaveCells
                 fN *= -1;
             }
 
-            // Is any vertex of the cell on the wrong side of the
-            // plane of this face?
+            // Is the centre of any other face of the cell on the
+            // wrong side of the plane of this face?
 
-            forAll(cPoints, cPtI)
+            forAll(cFaces, j)
             {
-                label ptI = cPoints[cPtI];
-
-                // Skip points that are on this face
-                if (findIndex(f, ptI) > -1)
+                if (j != i)
                 {
-                    continue;
-                }
+                    label fJ = cFaces[j];
 
-                const point& pt = pts[ptI];
+                    const point& pt = fCentres[fJ];
 
-                // 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.
+                    // 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);
+                    vector pC = (pt - fC);
 
-                pC /= max(mag(pC), VSMALL);
+                    pC /= max(mag(pC), VSMALL);
 
-                if ((pC & fN) > -planarCosAngle_)
-                {
-                    // Concave or planar face
-
-                    concave = true;
-
-                    if (setPtr)
+                    if ((pC & fN) > -planarCosAngle_)
                     {
-                        setPtr->insert(cellI);
+                        // Concave or planar face
+
+                        concave = true;
+
+                        if (setPtr)
+                        {
+                            setPtr->insert(cellI);
+                        }
+
+                        nConcaveCells++;
+
+                        break;
                     }
-
-                    nConcaveCells++;
-
-                    break;
                 }
             }
         }
@@ -2197,8 +2189,8 @@ bool Foam::primitiveMesh::checkConcaveCells
     {
         if (debug || report)
         {
-            Info<< " ***Concave cells found, number of cells: "
-                << nConcaveCells << endl;
+            Info<< " ***Concave cells (using face planes) found,"
+                << " number of cells: " << nConcaveCells << endl;
         }
 
         return true;