zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

primitiveMesh: Utilise geometry calculations defined in the face class

Browse Source
This commit is contained in:
Will Bainbridge
2023-09-05 10:24:33 +01:00
parent 53d92ef476
commit c73cbf6160
1 changed files with 4 additions and 66 deletions
Download Patch File Download Diff File Expand all files Collapse all files

70
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C
View File

@ -85,73 +85,11 @@ void Foam::primitiveMesh::makeFaceCentresAndAreas
forAll(fs, facei)
{
const labelList& f = fs[facei];
label nPoints = f.size();
// If the face is a triangle, do a direct calculation for efficiency
// and to avoid round-off error-related problems
if (nPoints == 3)
{
fCtrs[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
fAreas[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
}
// For more complex faces, decompose into triangles
else
{
// Compute an estimate of the centre as the average of the points
point pAvg = p[f[0]];
for (label pi = 1; pi < nPoints; pi++)
{
pAvg += p[f[pi]];
}
pAvg /= nPoints;
// Compute the face area normal and unit normal by summing up the
// normals of the triangles formed by connecting each edge to the
// point average.
vector sumA = Zero;
forAll(f, i)
{
const vector a =
(p[f[f.fcIndex(i)]] - p[f[i]])^(pAvg - p[f[i]]);
sumA += a;
}
const vector sumAHat = normalised(sumA);
// Compute the area-weighted sum of the triangle centres. Note use
// the triangle area projected in the direction of the face normal
// as the weight, *not* the triangle area magnitude. Only the
// former makes the calculation independent of the initial estimate.
scalar sumAn = 0.0;
vector sumAnc = Zero;
forAll(f, i)
{
const vector a =
(p[f[f.fcIndex(i)]] - p[f[i]])^(pAvg - p[f[i]]);
const vector c = p[f[i]] + p[f[f.fcIndex(i)]] + pAvg;
const scalar an = a & sumAHat;
sumAn += an;
sumAnc += an*c;
}
// Complete calculating centres and areas. If the face is too small
// for the sums to be reliably divided then just set the centre to
// the initial estimate.
if (sumAn > vSmall)
{
fCtrs[facei] = (1.0/3.0)*sumAnc/sumAn;
}
else
{
fCtrs[facei] = pAvg;
}
fAreas[facei] = 0.5*sumA;
}
const Tuple2<vector, point> areaAndCentre =
face::areaAndCentre(UIndirectList<point>(p, fs[facei]));
fCtrs[facei] = areaAndCentre.second();
fAreas[facei] = areaAndCentre.first();
magfAreas[facei] = max(mag(fAreas[facei]), rootVSmall);
}
}