Skip to content
Snippets Groups Projects
Commit cf5ae6d6 authored by mattijs's avatar mattijs
Browse files

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.
parent d947aa9c
No related branches found
No related tags found
No related merge requests found
...@@ -2111,7 +2111,6 @@ bool Foam::primitiveMesh::checkConcaveCells ...@@ -2111,7 +2111,6 @@ bool Foam::primitiveMesh::checkConcaveCells
const labelList& fOwner = faceOwner(); const labelList& fOwner = faceOwner();
const vectorField& fAreas = faceAreas(); const vectorField& fAreas = faceAreas();
const pointField& fCentres = faceCentres(); const pointField& fCentres = faceCentres();
const pointField& pts = points();
label nConcaveCells = 0; label nConcaveCells = 0;
...@@ -2119,8 +2118,6 @@ bool Foam::primitiveMesh::checkConcaveCells ...@@ -2119,8 +2118,6 @@ bool Foam::primitiveMesh::checkConcaveCells
{ {
const cell& cFaces = c[cellI]; const cell& cFaces = c[cellI];
const labelList& cPoints = cellPoints()[cellI];
bool concave = false; bool concave = false;
forAll(cFaces, i) forAll(cFaces, i)
...@@ -2134,8 +2131,6 @@ bool Foam::primitiveMesh::checkConcaveCells ...@@ -2134,8 +2131,6 @@ bool Foam::primitiveMesh::checkConcaveCells
const point& fC = fCentres[fI]; const point& fC = fCentres[fI];
const face& f = faces()[fI];
vector fN = fAreas[fI]; vector fN = fAreas[fI];
fN /= max(mag(fN), VSMALL); fN /= max(mag(fN), VSMALL);
...@@ -2147,45 +2142,42 @@ bool Foam::primitiveMesh::checkConcaveCells ...@@ -2147,45 +2142,42 @@ bool Foam::primitiveMesh::checkConcaveCells
fN *= -1; fN *= -1;
} }
// Is any vertex of the cell on the wrong side of the // Is the centre of any other face of the cell on the
// plane of this face? // wrong side of the plane of this face?
forAll(cPoints, cPtI) forAll(cFaces, j)
{ {
label ptI = cPoints[cPtI]; if (j != i)
// Skip points that are on this face
if (findIndex(f, ptI) > -1)
{ {
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 // If the cell is concave, the point will be on the
// positive normal side of the plane of f, defined by // positive normal side of the plane of f, defined by
// its centre and normal, and the angle between (pt - // its centre and normal, and the angle between (pt -
// fC) and fN will be less than 90 degrees, so the dot // fC) and fN will be less than 90 degrees, so the dot
// product will be positive. // 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_) if ((pC & fN) > -planarCosAngle_)
{ {
// Concave or planar face // Concave or planar face
concave = true; concave = true;
if (setPtr) if (setPtr)
{ {
setPtr->insert(cellI); setPtr->insert(cellI);
} }
nConcaveCells++; nConcaveCells++;
break; break;
}
} }
} }
} }
...@@ -2197,8 +2189,8 @@ bool Foam::primitiveMesh::checkConcaveCells ...@@ -2197,8 +2189,8 @@ bool Foam::primitiveMesh::checkConcaveCells
{ {
if (debug || report) if (debug || report)
{ {
Info<< " ***Concave cells found, number of cells: " Info<< " ***Concave cells (using face planes) found,"
<< nConcaveCells << endl; << " number of cells: " << nConcaveCells << endl;
} }
return true; return true;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment