diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.C b/src/OpenFOAM/meshes/meshShapes/face/face.C index 3f1ba47dee..5162fc82ac 100644 --- a/src/OpenFOAM/meshes/meshShapes/face/face.C +++ b/src/OpenFOAM/meshes/meshShapes/face/face.C @@ -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 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)); } diff --git a/src/OpenFOAM/meshes/meshShapes/face/faceTemplates.C b/src/OpenFOAM/meshes/meshShapes/face/faceTemplates.C index 06f0f9c19d..c5f72f212b 100644 --- a/src/OpenFOAM/meshes/meshShapes/face/faceTemplates.C +++ b/src/OpenFOAM/meshes/meshShapes/face/faceTemplates.C @@ -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 Type Foam::face::average ( - const pointField& meshPoints, + const pointField& ps, const Field& 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 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; } } diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C index 7f17a973ce..c83dde8fee 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C @@ -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; } } }