From bb1a2476fddd7670d81c7bdb0e48170cd3f23fea Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Wed, 31 Mar 2010 18:10:05 +0100
Subject: [PATCH] ENH: Change concave cell checking to be based on face centres
 and face planes, rather than a vertex-by-vertex check.  Assesses concavity as
 relevant to the tracking.

---
 .../meshes/primitiveMesh/primitiveMesh.H      |  2 +-
 .../primitiveMeshCheck/primitiveMeshCheck.C   | 60 ++++++++-----------
 2 files changed, 27 insertions(+), 35 deletions(-)

diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
index b1b54c9aa5..718df480e0 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
@@ -649,7 +649,7 @@ public:
                     labelHashSet* setPtr = NULL
                 ) const;
 
-                //- Check for concave cells
+                //- Check for concave cells by the planes of faces
                 bool checkConcaveCells
                 (
                     const bool report = false,
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
index 30b051d1e9..7dcfceaf10 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;
                 }
             }
         }