From f85edb01fe02ffb4e88be01ceff958c65a2becb4 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Tue, 18 Feb 2020 12:24:43 +0000 Subject: [PATCH] face, primitiveMesh: Corrected face-centre calculations 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 --- src/OpenFOAM/meshes/meshShapes/face/face.C | 159 ++++++++---------- .../meshes/meshShapes/face/faceTemplates.C | 96 ++++++----- .../primitiveMeshFaceCentresAndAreas.C | 64 ++++--- 3 files changed, 168 insertions(+), 151 deletions(-) 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; } } }