ENH: add primitiveMeshTools support for face lists

- allows reuse by finiteArea, for example.
- simplify edge looping with face thisLabel/nextLabel method

ENH: additional storage checks for mesh weights (faMesh + fvMesh)

- allow finite-area field decomposition without edge weights.

STYLE: use tmp New in various places. Simpler updateGeom check

src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,7 +31,7 @@ License
 #include "pyramidPointFaceRef.H"
 #include "primitiveMeshTools.H"
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
 Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceOrthogonality
 (
@@ -44,8 +44,8 @@ Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceOrthogonality
     const labelList& nei = mesh.faceNeighbour();
     const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 
-    tmp<scalarField> tortho(new scalarField(mesh.nFaces(), 1.0));
+    auto tortho = tmp<scalarField>::New(mesh.nFaces(), scalar(1));
-    scalarField& ortho = tortho.ref();
+    auto& ortho = tortho.ref();
 
     // Internal faces
     forAll(nei, facei)
@@ -64,9 +64,8 @@ Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceOrthogonality
     pointField neighbourCc;
     syncTools::swapBoundaryCellPositions(mesh, cc, neighbourCc);
 
-    forAll(pbm, patchi)
+    for (const polyPatch& pp : pbm)
     {
-        const polyPatch& pp = pbm[patchi];
         if (pp.coupled())
         {
             forAll(pp, i)
@@ -99,16 +98,17 @@ Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceSkewness
 {
     const labelList& own = mesh.faceOwner();
     const labelList& nei = mesh.faceNeighbour();
+    const faceList& faces = mesh.faces();
     const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 
-    tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
+    auto tskew = tmp<scalarField>::New(mesh.nFaces());
-    scalarField& skew = tskew.ref();
+    auto& skew = tskew.ref();
 
     forAll(nei, facei)
     {
         skew[facei] = primitiveMeshTools::faceSkewness
         (
-            mesh,
+            faces,
             p,
             fCtrs,
             fAreas,
@@ -126,9 +126,8 @@ Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceSkewness
     pointField neighbourCc;
     syncTools::swapBoundaryCellPositions(mesh, cellCtrs, neighbourCc);
 
-    forAll(pbm, patchi)
+    for (const polyPatch& pp : pbm)
     {
-        const polyPatch& pp = pbm[patchi];
         if (pp.coupled())
         {
             forAll(pp, i)
@@ -138,7 +137,7 @@ Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceSkewness
 
                 skew[facei] = primitiveMeshTools::faceSkewness
                 (
-                    mesh,
+                    faces,
                     p,
                     fCtrs,
                     fAreas,
@@ -157,7 +156,7 @@ Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceSkewness
 
                 skew[facei] = primitiveMeshTools::boundaryFaceSkewness
                 (
-                    mesh,
+                    faces,
                     p,
                     fCtrs,
                     fAreas,
@@ -185,8 +184,8 @@ Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceWeights
     const labelList& nei = mesh.faceNeighbour();
     const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 
-    tmp<scalarField> tweight(new scalarField(mesh.nFaces(), 1.0));
+    auto tweight = tmp<scalarField>::New(mesh.nFaces(), scalar(1));
-    scalarField& weight = tweight.ref();
+    auto& weight = tweight.ref();
 
     // Internal faces
     forAll(nei, facei)
@@ -206,9 +205,8 @@ Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceWeights
     pointField neiCc;
     syncTools::swapBoundaryCellPositions(mesh, cellCtrs, neiCc);
 
-    forAll(pbm, patchi)
+    for (const polyPatch& pp : pbm)
     {
-        const polyPatch& pp = pbm[patchi];
         if (pp.coupled())
         {
             forAll(pp, i)
@@ -241,8 +239,8 @@ Foam::tmp<Foam::scalarField> Foam::polyMeshTools::volRatio
     const labelList& nei = mesh.faceNeighbour();
     const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 
-    tmp<scalarField> tratio(new scalarField(mesh.nFaces(), 1.0));
+    auto tratio = tmp<scalarField>::New(mesh.nFaces(), scalar(1));
-    scalarField& ratio = tratio.ref();
+    auto& ratio = tratio.ref();
 
     // Internal faces
     forAll(nei, facei)
@@ -259,9 +257,8 @@ Foam::tmp<Foam::scalarField> Foam::polyMeshTools::volRatio
     scalarField neiVol;
     syncTools::swapBoundaryCellList(mesh, vol, neiVol);
 
-    forAll(pbm, patchi)
+    for (const polyPatch& pp : pbm)
     {
-        const polyPatch& pp = pbm[patchi];
         if (pp.coupled())
         {
             forAll(pp, i)
@@ -305,10 +302,8 @@ Foam::polyMesh::readUpdateState Foam::polyMeshTools::combine
     {
         return state1;
     }
-    else
+
-    {
+    return state0;
-        return state0;
-    }
 }
 
 

src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H

@@ -24,7 +24,7 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Namespace
+Class
     Foam::polyMeshTools
 
 Description
@@ -36,8 +36,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef polyMeshTools_H
+#ifndef Foam_polyMeshTools_H
-#define polyMeshTools_H
+#define Foam_polyMeshTools_H
 
 #include "polyMesh.H"
 #include "primitiveMeshTools.H"
@@ -48,18 +48,17 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                        Namespace polyMeshTools Declaration
+                        Class polyMeshTools Declaration
 \*---------------------------------------------------------------------------*/
 
 class polyMeshTools
 :
     public primitiveMeshTools
 {
-
 public:
 
-    //- Generate orthogonality field. (1 for fully orthogonal, < 1 for
+    //- Generate orthogonality field.
-    //  non-orthogonal)
+    //- (1 for fully orthogonal, < 1 for non-orthogonal)
     static tmp<scalarField> faceOrthogonality
     (
         const polyMesh& mesh,
@@ -93,8 +92,8 @@ public:
         const scalarField& vol
     );
 
-    //- Combine readUpdateState. topo change trumps geom-only
+    //- Combine readUpdateState.
-    //  change etc.
+    //- topo change trumps geom-only change etc.
     static polyMesh::readUpdateState combine
     (
         const polyMesh::readUpdateState& state0,

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C

@@ -156,8 +156,8 @@ Foam::PatchTools::pointNormals
     // 1. Start off with local normals (note:without calculating pointNormals
     //    to avoid them being stored)
 
-    tmp<pointField> textrudeN(new pointField(p.nPoints(), Zero));
+    auto textrudeN = tmp<pointField>::New(p.nPoints(), Zero);
-    pointField& extrudeN = textrudeN.ref();
+    auto& extrudeN = textrudeN.ref();
     {
         const faceList& localFaces = p.localFaces();
         const vectorField& faceNormals = p.faceNormals();

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C

@@ -333,14 +333,14 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMesh::movePoints
     }
 
     // Create swept volumes
-    const faceList& f = faces();
+    const faceList& fcs = faces();
 
-    tmp<scalarField> tsweptVols(new scalarField(f.size()));
+    auto tsweptVols = tmp<scalarField>::New(fcs.size());
-    scalarField& sweptVols = tsweptVols.ref();
+    auto& sweptVols = tsweptVols.ref();
 
-    forAll(f, facei)
+    forAll(fcs, facei)
     {
-        sweptVols[facei] = f[facei].sweptVol(oldPoints, newPoints);
+        sweptVols[facei] = fcs[facei].sweptVol(oldPoints, newPoints);
     }
 
     // Force recalculation of all geometric data with new points
@@ -364,20 +364,15 @@ const Foam::cellShapeList& Foam::primitiveMesh::cellShapes() const
 
 void Foam::primitiveMesh::updateGeom()
 {
-    if (!faceCentresPtr_)
+    if (!faceCentresPtr_ || !faceAreasPtr_)
     {
+        // These are always calculated in tandem, but only once
         calcFaceCentresAndAreas();
     }
-    if (!faceAreasPtr_)
+
-    {
+    if (!cellCentresPtr_ || !cellVolumesPtr_)
-        calcFaceCentresAndAreas();
-    }
-    if (!cellCentresPtr_)
-    {
-        calcCellCentresAndVols();
-    }
-    if (!cellVolumesPtr_)
     {
+        // These are always calculated in tandem, but only once
         calcCellCentresAndVols();
     }
 }

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H

@@ -50,8 +50,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef primitiveMesh_H
+#ifndef Foam_primitiveMesh_H
-#define primitiveMesh_H
+#define Foam_primitiveMesh_H
 
 #include "DynamicList.H"
 #include "edgeList.H"
@@ -188,7 +188,7 @@ class primitiveMesh
         void operator=(const primitiveMesh&) = delete;
 
 
-        // Topological calculations
+        // Topological Calculations
 
             //- Calculate cell shapes
             void calcCellShapes() const;
@@ -809,8 +809,8 @@ public:
             inline bool hasPointPoints() const noexcept;
             inline bool hasCellPoints() const noexcept;
             inline bool hasCellCentres() const noexcept;
-            inline bool hasFaceCentres() const noexcept;
             inline bool hasCellVolumes() const noexcept;
+            inline bool hasFaceCentres() const noexcept;
             inline bool hasFaceAreas() const noexcept;
 
             // On-the-fly addressing calculation. These functions return either

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C

@@ -39,16 +39,15 @@ void Foam::primitiveMesh::calcCellCentresAndVols() const
     if (debug)
     {
         Pout<< "primitiveMesh::calcCellCentresAndVols() : "
-            << "Calculating cell centres and cell volumes"
+            << "Calculating cell centres and volumes"
             << endl;
     }
 
-    // It is an error to attempt to recalculate cellCentres
+    // These are always calculated in tandem, but only once.
-    // if the pointer is already set
     if (cellCentresPtr_ || cellVolumesPtr_)
     {
         FatalErrorInFunction
-            << "Cell centres or cell volumes already calculated"
+            << "Cell centres or volumes already calculated"
             << abort(FatalError);
     }
 
@@ -73,7 +72,7 @@ void Foam::primitiveMesh::calcCellCentresAndVols() const
     if (debug)
     {
         Pout<< "primitiveMesh::calcCellCentresAndVols() : "
-            << "Finished calculating cell centres and cell volumes"
+            << "Finished calculating cell centres and volumes"
             << endl;
     }
 }

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,21 +32,23 @@ License
 #include "pyramidPointFaceRef.H"
 #include "PrecisionAdaptor.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
 void Foam::primitiveMeshTools::makeFaceCentresAndAreas
 (
-    const primitiveMesh& mesh,
+    const UList<face>& fcs,
     const pointField& p,
     vectorField& fCtrs,
     vectorField& fAreas
 )
 {
-    const faceList& fs = mesh.faces();
+    // Safety first - ensure properly sized
+    fCtrs.resize_nocopy(fcs.size());
+    fAreas.resize_nocopy(fcs.size());
 
-    forAll(fs, facei)
+    forAll(fcs, facei)
     {
-        const labelList& f = fs[facei];
+        const face& f = fcs[facei];
         const label nPoints = f.size();
 
         // If the face is a triangle, do a direct calculation for efficiency
@@ -59,26 +61,25 @@ void Foam::primitiveMeshTools::makeFaceCentresAndAreas
         else
         {
             solveVector sumN = Zero;
-            solveScalar sumA = 0.0;
+            solveScalar sumA = Zero;
             solveVector sumAc = Zero;
 
             solveVector fCentre = p[f[0]];
-            for (label pi = 1; pi < nPoints; pi++)
+            for (label pi = 1; pi < nPoints; ++pi)
             {
                 fCentre += solveVector(p[f[pi]]);
             }
-
             fCentre /= nPoints;
 
-            for (label pi = 0; pi < nPoints; pi++)
+            for (label pi = 0; pi < nPoints; ++pi)
             {
-                const label nextPi(pi == nPoints-1 ? 0 : pi+1);
+                const solveVector thisPoint(p[f.thisLabel(pi)]);
-                const solveVector nextPoint(p[f[nextPi]]);
+                const solveVector nextPoint(p[f.nextLabel(pi)]);
-                const solveVector thisPoint(p[f[pi]]);
 
                 solveVector c = thisPoint + nextPoint + fCentre;
                 solveVector n = (nextPoint - thisPoint)^(fCentre - thisPoint);
                 solveScalar a = mag(n);
+
                 sumN += n;
                 sumA += a;
                 sumAc += a*c;
@@ -101,6 +102,18 @@ void Foam::primitiveMeshTools::makeFaceCentresAndAreas
 }
 
 
+void Foam::primitiveMeshTools::makeFaceCentresAndAreas
+(
+    const primitiveMesh& mesh,
+    const pointField& p,
+    vectorField& fCtrs,
+    vectorField& fAreas
+)
+{
+    makeFaceCentresAndAreas(mesh.faces(), p, fCtrs, fAreas);
+}
+
+
 void Foam::primitiveMeshTools::makeCellCentresAndVols
 (
     const primitiveMesh& mesh,
@@ -117,7 +130,7 @@ void Foam::primitiveMeshTools::makeCellCentresAndVols
 
     // Clear the fields for accumulation
     cellCtrs = Zero;
-    cellVols = 0.0;
+    cellVols = Zero;
 
     const labelList& own = mesh.faceOwner();
     const labelList& nei = mesh.faceNeighbour();
@@ -201,7 +214,7 @@ void Foam::primitiveMeshTools::makeCellCentresAndVols
 
 Foam::scalar Foam::primitiveMeshTools::faceSkewness
 (
-    const primitiveMesh& mesh,
+    const UList<face>& fcs,
     const pointField& p,
     const vectorField& fCtrs,
     const vectorField& fAreas,
@@ -211,6 +224,8 @@ Foam::scalar Foam::primitiveMeshTools::faceSkewness
     const point& neiCc
 )
 {
+    const face& f = fcs[facei];
+
     vector Cpf = fCtrs[facei] - ownCc;
     vector d = neiCc - ownCc;
 
@@ -224,7 +239,6 @@ Foam::scalar Foam::primitiveMeshTools::faceSkewness
     // from the face centre to the edge of the face in the direction
     // of the skewness
     scalar fd = 0.2*mag(d) + ROOTVSMALL;
-    const face& f = mesh.faces()[facei];
     forAll(f, pi)
     {
         fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[facei])));
@@ -235,6 +249,60 @@ Foam::scalar Foam::primitiveMeshTools::faceSkewness
 }
 
 
+Foam::scalar Foam::primitiveMeshTools::boundaryFaceSkewness
+(
+    const UList<face>& fcs,
+    const pointField& p,
+    const vectorField& fCtrs,
+    const vectorField& fAreas,
+
+    const label facei,
+    const point& ownCc
+)
+{
+    const face& f = fcs[facei];
+
+    vector Cpf = fCtrs[facei] - ownCc;
+
+    vector normal = normalised(fAreas[facei]);
+    vector d = normal*(normal & Cpf);
+
+    // Skewness vector
+    vector sv =
+        Cpf
+      - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL))*d;
+    vector svHat = sv/(mag(sv) + ROOTVSMALL);
+
+    // Normalisation distance calculated as the approximate distance
+    // from the face centre to the edge of the face in the direction
+    // of the skewness
+    scalar fd = 0.4*mag(d) + ROOTVSMALL;
+    forAll(f, pi)
+    {
+        fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[facei])));
+    }
+
+    // Normalised skewness
+    return mag(sv)/fd;
+}
+
+
+Foam::scalar Foam::primitiveMeshTools::faceSkewness
+(
+    const primitiveMesh& mesh,
+    const pointField& p,
+    const vectorField& fCtrs,
+    const vectorField& fAreas,
+
+    const label facei,
+    const point& ownCc,
+    const point& neiCc
+)
+{
+    return faceSkewness(mesh.faces(), p, fCtrs, fAreas, facei, ownCc, neiCc);
+}
+
+
 Foam::scalar Foam::primitiveMeshTools::boundaryFaceSkewness
 (
     const primitiveMesh& mesh,
@@ -246,30 +314,7 @@ Foam::scalar Foam::primitiveMeshTools::boundaryFaceSkewness
     const point& ownCc
 )
 {
-    vector Cpf = fCtrs[facei] - ownCc;
+    return boundaryFaceSkewness(mesh.faces(), p, fCtrs, fAreas, facei, ownCc);
-
-    vector normal = normalised(fAreas[facei]);
-    vector d = normal*(normal & Cpf);
-
-
-    // Skewness vector
-    vector sv =
-        Cpf
-      - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL))*d;
-    vector svHat = sv/(mag(sv) + ROOTVSMALL);
-
-    // Normalisation distance calculated as the approximate distance
-    // from the face centre to the edge of the face in the direction
-    // of the skewness
-    scalar fd = 0.4*mag(d) + ROOTVSMALL;
-    const face& f = mesh.faces()[facei];
-    forAll(f, pi)
-    {
-        fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[facei])));
-    }
-
-    // Normalised skewness
-    return mag(sv)/fd;
 }
 
 
@@ -298,8 +343,8 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceOrthogonality
     const labelList& own = mesh.faceOwner();
     const labelList& nei = mesh.faceNeighbour();
 
-    tmp<scalarField> tortho(new scalarField(mesh.nInternalFaces()));
+    auto tortho = tmp<scalarField>::New(mesh.nInternalFaces());
-    scalarField& ortho = tortho.ref();
+    auto& ortho = tortho.ref();
 
     // Internal faces
     forAll(nei, facei)
@@ -327,15 +372,17 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceSkewness
 {
     const labelList& own = mesh.faceOwner();
     const labelList& nei = mesh.faceNeighbour();
+    const faceList& faces = mesh.faces();
 
-    tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
+    auto tskew = tmp<scalarField>::New(mesh.nFaces());
-    scalarField& skew = tskew.ref();
+    auto& skew = tskew.ref();
 
-    forAll(nei, facei)
+    // Internal faces
+    for (label facei = 0; facei < mesh.nInternalFaces(); ++facei)
     {
         skew[facei] = faceSkewness
         (
-            mesh,
+            faces,
             p,
             fCtrs,
             fAreas,
@@ -350,11 +397,11 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceSkewness
     // Boundary faces: consider them to have only skewness error.
     // (i.e. treat as if mirror cell on other side)
 
-    for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
+    for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei)
     {
         skew[facei] = boundaryFaceSkewness
         (
-            mesh,
+            faces,
             p,
             fCtrs,
             fAreas,
@@ -513,15 +560,14 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceConcavity
     vectorField faceNormals(faceAreas);
     faceNormals /= mag(faceNormals) + ROOTVSMALL;
 
-    tmp<scalarField> tfaceAngles(new scalarField(mesh.nFaces()));
+    auto tfaceAngles = tmp<scalarField>::New(mesh.nFaces());
-    scalarField& faceAngles = tfaceAngles.ref();
+    auto&& faceAngles = tfaceAngles.ref();
-
 
     forAll(fcs, facei)
     {
         const face& f = fcs[facei];
 
-        // Get edge from f[0] to f[size-1];
+        // Normalized vector from f[size-1] to f[0];
         vector ePrev(p[f.first()] - p[f.last()]);
         scalar magEPrev = mag(ePrev);
         ePrev /= magEPrev + ROOTVSMALL;
@@ -530,11 +576,8 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceConcavity
 
         forAll(f, fp0)
         {
-            // Get vertex after fp
-            label fp1 = f.fcIndex(fp0);
-
             // Normalized vector between two consecutive points
-            vector e10(p[f[fp1]] - p[f[fp0]]);
+            vector e10(p[f.nextLabel(fp0)] - p[f.thisLabel(fp0)]);
             scalar magE10 = mag(e10);
             e10 /= magE10 + ROOTVSMALL;
 
@@ -584,8 +627,8 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceFlatness
     // primitiveMeshFaceCentresAndAreas.C)
     scalarField magAreas(mag(faceAreas));
 
-    tmp<scalarField> tfaceFlatness(new scalarField(mesh.nFaces(), 1.0));
+    auto tfaceFlatness = tmp<scalarField>::New(mesh.nFaces(), scalar(1));
-    scalarField& faceFlatness = tfaceFlatness.ref();
+    auto& faceFlatness = tfaceFlatness.ref();
 
     forAll(fcs, facei)
     {
@@ -641,8 +684,8 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::cellDeterminant
         }
     }
 
-    tmp<scalarField> tcellDeterminant(new scalarField(mesh.nCells()));
+    auto tcellDeterminant = tmp<scalarField>::New(mesh.nCells());
-    scalarField& cellDeterminant = tcellDeterminant.ref();
+    auto& cellDeterminant = tcellDeterminant.ref();
 
     const cellList& c = mesh.cells();
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,7 +24,7 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Namespace
+Class
     Foam::primitiveMeshTools
 
 Description
@@ -35,8 +35,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef primitiveMeshTools_H
+#ifndef Foam_primitiveMeshTools_H
-#define primitiveMeshTools_H
+#define Foam_primitiveMeshTools_H
 
 #include "bitSet.H"
 #include "primitiveFields.H"
@@ -47,24 +47,36 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declarations
+// Forward Declarations
+class face;
 class primitiveMesh;
 
 /*---------------------------------------------------------------------------*\
-                   Namespace primitiveMeshTools Declaration
+                     Class primitiveMeshTools Declaration
 \*---------------------------------------------------------------------------*/
 
 class primitiveMeshTools
 {
 public:
 
-    //- Calculate face centres and areas
+    //- Calculate face centres and areas for specified faces.
+    //  Adjusts the lengths of centres and area normals if required.
     static void makeFaceCentresAndAreas
     (
-        const primitiveMesh& mesh,
+        const UList<face>& faces,   //!< The faces to consider
-        const pointField& p,
+        const pointField& p,        //!< Support points for faces
-        vectorField& fCtrs,
+        vectorField& fCtrs,         //!< [out] The face centres
-        vectorField& fAreas
+        vectorField& fAreas         //!< [out] The face area normals
+    );
+
+    //- Calculate face centres and areas for all mesh faces.
+    //  Adjusts the lengths of centres and area normals if required.
+    static void makeFaceCentresAndAreas
+    (
+        const primitiveMesh& mesh,  //!< The mesh faces
+        const pointField& p,        //!< Support points for faces
+        vectorField& fCtrs,         //!< [out] The face centres
+        vectorField& fAreas         //!< [out] The face area normals
     );
 
     //- Calculate cell centres and volumes from face properties
@@ -117,7 +129,7 @@ public:
     );
 
     //- Generate face concavity field. Returns per face the (sin of the)
-    //  most concave angle between two consecutive edges
+    //- most concave angle between two consecutive edges
     static tmp<scalarField> faceConcavity
     (
         const scalar maxSin,
@@ -127,8 +139,8 @@ public:
     );
 
     //- Generate face flatness field. Compares the individual triangles'
-    //  normals against the face average normal. Between 0 (fully warped)
+    //- normals against the face average normal. Between 0 (fully warped)
-    //  and 1 (fully flat)
+    //- and 1 (fully flat)
     static tmp<scalarField> faceFlatness
     (
         const primitiveMesh& mesh,
@@ -138,7 +150,7 @@ public:
     );
 
     //- Generate edge alignment field. Is per face the minimum aligned edge
-    //  (does not use edge addressing)
+    //- (does not use edge addressing)
     static tmp<scalarField> edgeAlignment
     (
         const primitiveMesh& mesh,
@@ -161,8 +173,8 @@ public:
         //- Skewness of single face
         static scalar faceSkewness
         (
-            const primitiveMesh& mesh,
+            const UList<face>& faces,   //!< The faces to consider
-            const pointField& p,
+            const pointField& p,        //!< Support points for faces
             const vectorField& fCtrs,
             const vectorField& fAreas,
 
@@ -174,8 +186,33 @@ public:
         //- Skewness of single boundary face
         static scalar boundaryFaceSkewness
         (
-            const primitiveMesh& mesh,
+            const UList<face>& faces,   //!< The faces to consider
-            const pointField& p,
+            const pointField& p,        //!< Support points for faces
+            const vectorField& fCtrs,
+            const vectorField& fAreas,
+
+            const label facei,
+            const point& ownCc
+        );
+
+        //- Skewness of single face
+        static scalar faceSkewness
+        (
+            const primitiveMesh& mesh,  //!< The mesh faces
+            const pointField& p,        //!< Support points for faces
+            const vectorField& fCtrs,
+            const vectorField& fAreas,
+
+            const label facei,
+            const point& ownCc,
+            const point& neiCc
+        );
+
+        //- Skewness of single boundary face
+        static scalar boundaryFaceSkewness
+        (
+            const primitiveMesh& mesh,  //!< The mesh faces
+            const pointField& p,        //!< Support points for faces
             const vectorField& fCtrs,
             const vectorField& fAreas,
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshClear.C

@@ -106,21 +106,20 @@ void Foam::primitiveMesh::printAllocated() const
         Pout<< "    Cell-centres" << endl;
     }
 
-    if (faceCentresPtr_)
-    {
-        Pout<< "    Face-centres" << endl;
-    }
-
     if (cellVolumesPtr_)
     {
         Pout<< "    Cell-volumes" << endl;
     }
 
+    if (faceCentresPtr_)
+    {
+        Pout<< "    Face-centres" << endl;
+    }
+
     if (faceAreasPtr_)
     {
         Pout<< "    Face-areas" << endl;
     }
-
 }
 
 
@@ -134,8 +133,8 @@ void Foam::primitiveMesh::clearGeom()
     }
 
     deleteDemandDrivenData(cellCentresPtr_);
-    deleteDemandDrivenData(faceCentresPtr_);
     deleteDemandDrivenData(cellVolumesPtr_);
+    deleteDemandDrivenData(faceCentresPtr_);
     deleteDemandDrivenData(faceAreasPtr_);
 }
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C

@@ -42,16 +42,15 @@ void Foam::primitiveMesh::calcFaceCentresAndAreas() const
     if (debug)
     {
         Pout<< "primitiveMesh::calcFaceCentresAndAreas() : "
-            << "Calculating face centres and face areas"
+            << "Calculating face centres and areas"
             << endl;
     }
 
-    // It is an error to attempt to recalculate faceCentres
+    // These are always calculated in tandem, but only once.
-    // if the pointer is already set
     if (faceCentresPtr_ || faceAreasPtr_)
     {
         FatalErrorInFunction
-            << "Face centres or face areas already calculated"
+            << "Face centres and areas already calculated"
             << abort(FatalError);
     }
 
@@ -66,7 +65,7 @@ void Foam::primitiveMesh::calcFaceCentresAndAreas() const
     if (debug)
     {
         Pout<< "primitiveMesh::calcFaceCentresAndAreas() : "
-            << "Finished calculating face centres and face areas"
+            << "Finished calculating face centres and areas"
             << endl;
     }
 }

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshI.H

@@ -192,18 +192,18 @@ inline bool Foam::primitiveMesh::hasCellCentres() const noexcept
 }
 
 
-inline bool Foam::primitiveMesh::hasFaceCentres() const noexcept
-{
-    return bool(faceCentresPtr_);
-}
-
-
 inline bool Foam::primitiveMesh::hasCellVolumes() const noexcept
 {
     return bool(cellVolumesPtr_);
 }
 
 
+inline bool Foam::primitiveMesh::hasFaceCentres() const noexcept
+{
+    return bool(faceCentresPtr_);
+}
+
+
 inline bool Foam::primitiveMesh::hasFaceAreas() const noexcept
 {
     return bool(faceAreasPtr_);

src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2019 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,12 +50,12 @@ void Foam::polyMeshGeometry::updateFaceCentresAndAreas
     const labelList& changedFaces
 )
 {
-    const faceList& fs = mesh_.faces();
+    const faceList& fcs = mesh_.faces();
 
     for (const label facei : changedFaces)
     {
-        const labelList& f = fs[facei];
+        const face& f = fcs[facei];
-        label nPoints = f.size();
+        const label nPoints = f.size();
 
         // If the face is a triangle, do a direct calculation for efficiency
         // and to avoid round-off error-related problems
@@ -66,33 +66,41 @@ void Foam::polyMeshGeometry::updateFaceCentresAndAreas
         }
         else
         {
-            vector sumN = Zero;
+            solveVector sumN = Zero;
-            scalar sumA = Zero;
+            solveScalar sumA = Zero;
-            vector sumAc = Zero;
+            solveVector sumAc = Zero;
 
-            point fCentre = p[f[0]];
+            solveVector fCentre = p[f[0]];
             for (label pi = 1; pi < nPoints; ++pi)
             {
-                fCentre += p[f[pi]];
+                fCentre += solveVector(p[f[pi]]);
             }
-
             fCentre /= nPoints;
 
             for (label pi = 0; pi < nPoints; ++pi)
             {
-                const point& nextPoint = p[f[(pi + 1) % nPoints]];
+                const solveVector thisPoint(p[f.thisLabel(pi)]);
+                const solveVector nextPoint(p[f.nextLabel(pi)]);
 
-                vector c = p[f[pi]] + nextPoint + fCentre;
+                solveVector c = thisPoint + nextPoint + fCentre;
-                vector n = (nextPoint - p[f[pi]])^(fCentre - p[f[pi]]);
+                solveVector n = (nextPoint - thisPoint)^(fCentre - thisPoint);
-                scalar a = mag(n);
+                solveScalar a = mag(n);
 
                 sumN += n;
                 sumA += a;
                 sumAc += a*c;
             }
 
-            faceCentres_[facei] = (1.0/3.0)*sumAc/(sumA + VSMALL);
+            if (sumA < ROOTVSMALL)
-            faceAreas_[facei] = 0.5*sumN;
+            {
+                faceCentres_[facei] = fCentre;
+                faceAreas_[facei] = Zero;
+            }
+            else
+            {
+                faceCentres_[facei] = (1.0/3.0)*sumAc/sumA;
+                faceAreas_[facei] = 0.5*sumN;
+            }
         }
     }
 }
@@ -946,6 +954,7 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
 
     const labelList& own = mesh.faceOwner();
     const labelList& nei = mesh.faceNeighbour();
+    const faceList& faces = mesh.faces();
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
     // Calculate coupled cell centre
@@ -962,7 +971,7 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
         {
             scalar skewness = primitiveMeshTools::faceSkewness
             (
-                mesh,
+                faces,
                 points,
                 faceCentres,
                 faceAreas,
@@ -996,7 +1005,7 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
         {
             scalar skewness = primitiveMeshTools::faceSkewness
             (
-                mesh,
+                faces,
                 points,
                 faceCentres,
                 faceAreas,
@@ -1030,7 +1039,7 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
         {
             scalar skewness = primitiveMeshTools::boundaryFaceSkewness
             (
-                mesh,
+                faces,
                 points,
                 faceCentres,
                 faceAreas,
@@ -1072,7 +1081,7 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
 
         scalar skewness = primitiveMeshTools::faceSkewness
         (
-            mesh,
+            faces,
             points,
             faceCentres,
             faceAreas,
@@ -1449,18 +1458,15 @@ bool Foam::polyMeshGeometry::checkFaceAngles
 
         const vector faceNormal = normalised(faceAreas[facei]);
 
-        // Get edge from f[0] to f[size-1];
+        // Normalized vector from f[size-1] to f[0];
         vector ePrev(p[f.first()] - p[f.last()]);
         scalar magEPrev = mag(ePrev);
         ePrev /= magEPrev + VSMALL;
 
         forAll(f, fp0)
         {
-            // Get vertex after fp
-            label fp1 = f.fcIndex(fp0);
-
             // Normalized vector between two consecutive points
-            vector e10(p[f[fp1]] - p[f[fp0]]);
+            vector e10(p[f.nextLabel(fp0)] - p[f.thisLabel(fp0)]);
             scalar magE10 = mag(e10);
             e10 /= magE10 + VSMALL;
 

src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -61,18 +61,15 @@ bool Foam::combineFaces::convexFace
     // Get outwards pointing normal of f, only the sign matters.
     const vector areaNorm = f.areaNormal(points);
 
-    // Get edge from f[0] to f[size-1];
+    // Normalized vector from f[size-1] to f[0];
     vector ePrev(points[f.first()] - points[f.last()]);
     scalar magEPrev = mag(ePrev);
     ePrev /= magEPrev + VSMALL;
 
     forAll(f, fp0)
     {
-        // Get vertex after fp
-        label fp1 = f.fcIndex(fp0);
-
         // Normalized vector between two consecutive points
-        vector e10(points[f[fp1]] - points[f[fp0]]);
+        vector e10(points[f.nextLabel(fp0)] - points[f.thisLabel(fp0)]);
         scalar magE10 = mag(e10);
         e10 /= magE10 + VSMALL;
 
@@ -138,12 +135,10 @@ void Foam::combineFaces::regioniseFaces
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    forAll(cEdges, i)
+    for (const label edgei : cEdges)
     {
-        const label edgeI = cEdges[i];
-
         label f0, f1;
-        meshTools::getEdgeFaces(mesh_, celli, edgeI, f0, f1);
+        meshTools::getEdgeFaces(mesh_, celli, edgei, f0, f1);
 
         const vector& a0 = mesh_.faceAreas()[f0];
         const vector& a1 = mesh_.faceAreas()[f1];

src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C

@@ -59,10 +59,10 @@ Foam::edgeInterpolation::edgeInterpolation(const faMesh& fam)
     lPN_(nullptr),
     weightingFactors_(nullptr),
     differenceFactors_(nullptr),
-    orthogonal_(false),
     correctionVectors_(nullptr),
-    skew_(true),
+    skewCorrectionVectors_(nullptr),
-    skewCorrectionVectors_(nullptr)
+    orthogonal_(false),
+    skew_(true)
 {}
 
 

src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.H

@@ -38,8 +38,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef edgeInterpolation_H
+#ifndef Foam_edgeInterpolation_H
-#define edgeInterpolation_H
+#define Foam_edgeInterpolation_H
 
 #include "tmp.H"
 #include "scalar.H"
@@ -52,6 +52,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class polyMesh;
 
 /*---------------------------------------------------------------------------*\
@@ -60,7 +61,7 @@ class polyMesh;
 
 class edgeInterpolation
 {
-    // Private data
+    // Private Data
 
         // Reference to faMesh
         const faMesh& faMesh_;
@@ -77,20 +78,20 @@ class edgeInterpolation
             //- Face-gradient difference factors
             mutable edgeScalarField* differenceFactors_;
 
-            //- Is mesh orthogonal
-            mutable bool orthogonal_;
-
             //- Non-orthogonality correction vectors
             mutable edgeVectorField* correctionVectors_;
 
-            //- Is mesh skew
-            mutable bool skew_;
-
             //- Skew correction vectors
             mutable edgeVectorField* skewCorrectionVectors_;
 
+            //- Is mesh orthogonal
+            mutable bool orthogonal_;
 
-    // Private member functions
+            //- Is mesh skew
+            mutable bool skew_;
+
+
+    // Private Member Functions
 
         //- Construct geodesic distance between P and N
         void makeLPN() const;
@@ -107,15 +108,12 @@ class edgeInterpolation
         //- Construct skewness correction vectors
         void makeSkewCorrectionVectors() const;
 
-//         //- Construct Least-squares gradient vectors
-//         void makeLeastSquareVectors() const;
-
 
 protected:
 
-    // Protected member functions
+    // Protected Member Functions
 
-        // Storage management
+        // Storage Management
 
             //- Clear all geometry and addressing
             void clearOut();
@@ -137,10 +135,10 @@ public:
     ~edgeInterpolation();
 
 
-    // Member functions
+    // Member Functions
 
         //- Return mesh reference
-        const faMesh& mesh() const
+        const faMesh& mesh() const noexcept
         {
             return faMesh_;
         }
@@ -154,21 +152,29 @@ public:
         //- Return reference to difference factors array
         const edgeScalarField& deltaCoeffs() const;
 
-        //- Return whether mesh is orthogonal or not
-        bool orthogonal() const;
-
         //- Return reference to non-orthogonality correction vectors array
         const edgeVectorField& correctionVectors() const;
 
-        //- Return whether mesh is skew or not
-        bool skew() const;
-
         //- Return reference to skew vectors array
         const edgeVectorField& skewCorrectionVectors() const;
 
+        //- Return whether mesh is orthogonal or not
+        bool orthogonal() const;
+
+        //- Return whether mesh is skew or not
+        bool skew() const;
+
+
+    // Mesh Motion
 
         //- Do what is necessary if the mesh has moved
         bool movePoints() const;
+
+
+    // Storage Management
+
+        //- True if weights exist
+        bool hasWeights() const noexcept { return bool(weightingFactors_); }
 };
 
 

src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.H

@@ -35,8 +35,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef surfaceInterpolation_H
+#ifndef Foam_surfaceInterpolation_H
-#define surfaceInterpolation_H
+#define Foam_surfaceInterpolation_H
 
 #include "tmp.H"
 #include "scalar.H"
@@ -49,6 +49,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class fvMesh;
 class fvGeometryScheme;
 
@@ -58,7 +59,7 @@ class fvGeometryScheme;
 
 class surfaceInterpolation
 {
-    // Private data
+    // Private Data
 
         // Reference to fvMesh
         const fvMesh& mesh_;
@@ -108,7 +109,7 @@ public:
     virtual ~surfaceInterpolation();
 
 
-    // Member functions
+    // Member Functions
 
         //- Return reference to geometry calculation scheme
         virtual const fvGeometryScheme& geometry() const;
@@ -134,6 +135,12 @@ public:
 
         //- Update all geometric data
         virtual void updateGeom();
+
+
+    // Storage Management
+
+        //- Has weights
+        bool hasWeights() const noexcept { return bool(weights_); }
 };
 
 

src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C

@@ -763,18 +763,15 @@ bool Foam::primitiveMeshGeometry::checkFaceAngles
 
         const vector faceNormal = normalised(faceAreas[facei]);
 
-        // Get edge from f[0] to f[size-1];
+        // Normalized vector from f[size-1] to f[0];
         vector ePrev(p[f.first()] - p[f.last()]);
         scalar magEPrev = mag(ePrev);
         ePrev /= magEPrev + VSMALL;
 
         forAll(f, fp0)
         {
-            // Get vertex after fp
-            label fp1 = f.fcIndex(fp0);
-
             // Normalized vector between two consecutive points
-            vector e10(p[f[fp1]] - p[f[fp0]]);
+            vector e10(p[f.nextLabel(fp0)] - p[f.thisLabel(fp0)]);
             scalar magE10 = mag(e10);
             e10 /= magE10 + VSMALL;
 

src/parallel/decompose/faDecompose/faFieldDecomposer.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -126,7 +126,11 @@ processorAreaPatchFieldDecomposer
         mesh.edgeOwner(),
         mesh.edgeNeighbour(),
         addressingSlice,
-        mesh.weights().internalField()
+        (
+            mesh.hasWeights()
+          ? mesh.weights().primitiveField()
+          : scalarField::null()
+        )
     )
 {}
 
@@ -310,7 +314,9 @@ void Foam::faFieldDecomposer::reset(const faMesh& completeMesh)
     processorEdgePatchFieldDecomposerPtrs_.resize(nMappers);
 
     // Create weightings now - needed for proper parallel synchronization
-    (void)completeMesh.weights();
+    //// (void)completeMesh.weights();
+    // Disabled the above (2022-04-04)
+    // Use weights if they already exist, otherwise simply ignore
 
     // faPatches don't have their own start() - so these are invariant
     const labelList completePatchStarts