Commit bb1a2476 authored by graham's avatar graham
Browse files

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.
parent 3cf4d721
......@@ -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,
......
......@@ -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
if ((pC & fN) > -planarCosAngle_)
{
// Concave or planar face
concave = true;
concave = true;
if (setPtr)
{
setPtr->insert(cellI);
}
if (setPtr)
{
setPtr->insert(cellI);
}
nConcaveCells++;
nConcaveCells++;
break;
break;
}
}
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment