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

face, primitiveMesh: Corrected face-centre calculations

Browse Source 
Face centres are calculated by area-weighting the centres of triangles
that are formed by connecting each edge to a common point. The common
point is taken to be the average of all the face vertices, though in
principle the choice is arbitrary.

The areas used to perform the weighting are now taken to be the
projected areas of the triangles in the direction of the face normal
(i.e., the dot product between the triangle area and the face normal).
Previously the triangle area-magnitudes were used.

The new approach results in a centre that for flat faces is independent
of the choice of common point. It also means that concave faces have
contributions to the weighted sum from reversed triangles correctly
subtracted from the total. The centre of warped faces still changes with
the choice of common point, but that variation is now reduced to be
only in the direction of the face normal.

Preliminary results suggest a positive effect of this change on the
convergence of simulations on meshes with significantly distorted faces
and cells. The simpleFoam motorBike tutorial now converges with
residuals approximately half that previously observed.

Resolves bug report https://bugs.openfoam.org/view.php?id=1993
This commit is contained in:
Will Bainbridge
2020-02-18 12:24:43 +00:00
parent 5ce2130f1a
commit f85edb01fe
3 changed files with 168 additions and 151 deletions
Download Patch File Download Diff File Expand all files Collapse all files

159
src/OpenFOAM/meshes/meshShapes/face/face.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -485,134 +485,123 @@ void Foam::face::flip()
}
Foam::point Foam::face::centre(const pointField& points) const
Foam::point Foam::face::centre(const pointField& ps) const
{
// Calculate the centre by breaking the face into triangles and
// area-weighted averaging their centres
const label nPoints = size();
// If the face is a triangle, do a direct calculation
if (nPoints == 3)
if (size() == 3)
{
return
(1.0/3.0)
*(
points[operator[](0)]
+ points[operator[](1)]
+ points[operator[](2)]
ps[operator[](0)]
+ ps[operator[](1)]
+ ps[operator[](2)]
);
}
// For more complex faces, decompose into triangles ...
point centrePoint = Zero;
for (label pI=0; pI<nPoints; ++pI)
// Compute an estimate of the centre as the average of the points
point pAvg = Zero;
forAll(*this, pi)
{
centrePoint += points[operator[](pI)];
pAvg += ps[operator[](pi)];
}
centrePoint /= nPoints;
pAvg /= size();
scalar sumA = 0;
vector sumAc = Zero;
for (label pI=0; pI<nPoints; ++pI)
// 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(*this, pi)
{
const point& nextPoint = points[operator[]((pI + 1) % nPoints)];
const point& p = ps[operator[](pi)];
const point& pNext = ps[operator[](fcIndex(pi))];
// Calculate 3*triangle centre
const vector ttc
(
points[operator[](pI)]
+ nextPoint
+ centrePoint
);
const vector a = (pNext - p)^(pAvg - p);
// Calculate 2*triangle area
const scalar ta = Foam::mag
(
(points[operator[](pI)] - centrePoint)
^ (nextPoint - centrePoint)
);
sumA += a;
}
const vector sumAHat = normalised(sumA);
sumA += ta;
sumAc += ta*ttc;
// 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;
vector sumAnc = Zero;
forAll(*this, pi)
{
const point& p = ps[operator[](pi)];
const point& pNext = ps[operator[](fcIndex(pi))];
const vector a = (pNext - p)^(pAvg - p);
const vector c = p + pNext + pAvg;
const scalar an = a & sumAHat;
sumAn += an;
sumAnc += an*c;
}
if (sumA > vSmall)
// 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)
{
return sumAc/(3.0*sumA);
return (1.0/3.0)*sumAnc/sumAn;
}
else
{
return centrePoint;
return pAvg;
}
}
Foam::vector Foam::face::area(const pointField& p) const
Foam::vector Foam::face::area(const pointField& ps) const
{
const label nPoints = size();
// Calculate the area by summing the face triangle areas.
// Changed to deal with small concavity by using a central decomposition
//
// If the face is a triangle, do a direct calculation to avoid round-off
// error-related problems
//
if (nPoints == 3)
// If the face is a triangle, do a direct calculation
if (size() == 3)
{
return triPointRef
(
p[operator[](0)],
p[operator[](1)],
p[operator[](2)]
).area();
return
0.5
*(
(ps[operator[](1)] - ps[operator[](0)])
^(ps[operator[](2)] - ps[operator[](0)])
);
}
label pI;
// For more complex faces, decompose into triangles ...
point centrePoint = Zero;
for (pI = 0; pI < nPoints; ++pI)
// Compute an estimate of the centre as the average of the points
point pAvg = Zero;
forAll(*this, pi)
{
centrePoint += p[operator[](pI)];
pAvg += ps[operator[](pi)];
}
centrePoint /= nPoints;
pAvg /= size();
vector a = Zero;
point nextPoint = centrePoint;
for (pI = 0; pI < nPoints; ++pI)
// 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(*this, pi)
{
if (pI < nPoints - 1)
{
nextPoint = p[operator[](pI + 1)];
}
else
{
nextPoint = p[operator[](0)];
}
const point& p = ps[operator[](pi)];
const point& pNext = ps[operator[](fcIndex(pi))];
// Note: for best accuracy, centre point always comes last
//
a += triPointRef
(
p[operator[](pI)],
nextPoint,
centrePoint
).area();
const vector a = (pNext - p)^(pAvg - p);
sumA += a;
}
return a;
return 0.5*sumA;
}
Foam::vector Foam::face::normal(const pointField& points) const
{
const vector a = area(points);
const scalar maga = Foam::mag(a);
return maga > 0 ? a/maga : Zero;
return normalised(area(points));
}

96
src/OpenFOAM/meshes/meshShapes/face/faceTemplates.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -49,13 +49,10 @@ Foam::label Foam::face::triangles
template<class Type>
Type Foam::face::average
(
const pointField& meshPoints,
const pointField& ps,
const Field<Type>& fld
) const
{
// Calculate the average by breaking the face into triangles and
// area-weighted averaging their averages
// If the face is a triangle, do a direct calculation
if (size() == 3)
{
@ -68,51 +65,66 @@ Type Foam::face::average
);
}
label nPoints = size();
// For more complex faces, decompose into triangles ...
point centrePoint = Zero;
Type cf = Zero;
for (label pI=0; pI<nPoints; pI++)
// Compute an estimate of the centre and the field average as the average
// of the point values
point pAvg = Zero;
Type fldAvg = Zero;
forAll(*this, pi)
{
centrePoint += meshPoints[operator[](pI)];
cf += fld[operator[](pI)];
pAvg += ps[operator[](pi)];
fldAvg += fld[operator[](pi)];
}
pAvg /= size();
fldAvg /= size();
// 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(*this, pi)
{
const point& p = ps[operator[](pi)];
const point& pNext = ps[operator[](fcIndex(pi))];
const vector a = (pNext - p)^(pAvg - p);
sumA += a;
}
const vector sumAHat = normalised(sumA);
// Compute the area-weighted sum of the average field values on each
// triangle
scalar sumAn = 0;
Type sumAnf = Zero;
forAll(*this, pi)
{
const point& p = ps[operator[](pi)];
const point& pNext = ps[operator[](fcIndex(pi))];
const vector a = (pNext - p)^(pAvg - p);
const Type f =
fld[operator[](pi)]
+ fld[operator[](fcIndex(pi))]
+ fldAvg;
const scalar an = a & sumAHat;
sumAn += an;
sumAnf += an*f;
}
centrePoint /= nPoints;
cf /= nPoints;
scalar sumA = 0;
Type sumAf = Zero;
for (label pI=0; pI<nPoints; pI++)
// Complete calculating the average. If the face is too small for the sums
// to be reliably divided then just set the average to the initial
// estimate.
if (sumAn > vSmall)
{
// Calculate 3*triangle centre field value
Type ttcf =
(
fld[operator[](pI)]
+ fld[operator[]((pI + 1) % nPoints)]
+ cf
);
// Calculate 2*triangle area
scalar ta = Foam::mag
(
(meshPoints[operator[](pI)] - centrePoint)
^ (meshPoints[operator[]((pI + 1) % nPoints)] - centrePoint)
);
sumA += ta;
sumAf += ta*ttcf;
}
if (sumA > vSmall)
{
return sumAf/(3*sumA);
return (1.0/3.0)*sumAnf/sumAn;
}
else
{
return cf;
return fldAvg;
}
}

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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -91,45 +91,61 @@ void Foam::primitiveMesh::makeFaceCentresAndAreas
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
{
vector sumN = Zero;
scalar sumA = 0.0;
vector sumAc = Zero;
point fCentre = p[f[0]];
// Compute an estimate of the centre as the average of the points
point pAvg = p[f[0]];
for (label pi = 1; pi < nPoints; pi++)
{
fCentre += p[f[pi]];
pAvg += p[f[pi]];
}
pAvg /= nPoints;
fCentre /= nPoints;
for (label pi = 0; pi < nPoints; pi++)
// 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 point& nextPoint = p[f[(pi + 1) % nPoints]];
const vector a =
(p[f[f.fcIndex(i)]] - p[f[i]])^(pAvg - p[f[i]]);
vector c = p[f[pi]] + nextPoint + fCentre;
vector n = (nextPoint - p[f[pi]])^(fCentre - p[f[pi]]);
scalar a = mag(n);
sumN += n;
sumA += a;
sumAc += a*c;
}
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;
}
// This is to deal with zero-area faces. Mark very small faces
// to be detected in e.g., processorPolyPatch.
if (sumA < rootVSmall)
// 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] = fCentre;
fAreas[facei] = Zero;
fCtrs[facei] = (1.0/3.0)*sumAnc/sumAn;
}
else
{
fCtrs[facei] = (1.0/3.0)*sumAc/sumA;
fAreas[facei] = 0.5*sumN;
fCtrs[facei] = pAvg;
}
fAreas[facei] = 0.5*sumA;
}
}
}