From 683ff304cae957029b0388787c0d2fa1ca7730db Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Fri, 16 Nov 2018 09:18:09 +0100
Subject: [PATCH] ENH: update looping and use of bitSet for motionSmootherAlgo
 (#1075)

---
 .../motionSmoother/motionSmoother.H           |  19 +-
 .../motionSmoother/motionSmootherAlgo.C       | 104 ++----
 .../motionSmoother/motionSmootherAlgo.H       |   9 +-
 .../motionSmootherAlgoTemplates.C             |  62 ++--
 .../motionSmoother/motionSmootherData.H       |  10 +-
 .../polyMeshGeometry/polyMeshGeometry.C       | 327 ++++++------------
 6 files changed, 182 insertions(+), 349 deletions(-)

diff --git a/src/dynamicMesh/motionSmoother/motionSmoother.H b/src/dynamicMesh/motionSmoother/motionSmoother.H
index 580d6f18e67..2b456460049 100644
--- a/src/dynamicMesh/motionSmoother/motionSmoother.H
+++ b/src/dynamicMesh/motionSmoother/motionSmoother.H
@@ -93,8 +93,6 @@ class motionSmoother
     public motionSmootherAlgo
 
 {
-    // Private class
-
 public:
 
     ClassName("motionSmoother");
@@ -105,10 +103,10 @@ public:
         //  Reads displacement field (only boundary conditions used)
         motionSmoother
         (
-            polyMesh&,
-            pointMesh&,
-            indirectPrimitivePatch& pp,         // 'outside' points
-            const labelList& adaptPatchIDs,     // patches forming 'outside'
+            polyMesh& mesh,
+            pointMesh& pMesh,
+            indirectPrimitivePatch& pp,         //!< 'outside' points
+            const labelList& adaptPatchIDs,     //!< patches forming 'outside'
             const dictionary& paramDict
         );
 
@@ -116,13 +114,12 @@ public:
         //  displacementfield (only boundary conditions used)
         motionSmoother
         (
-            polyMesh&,
-            indirectPrimitivePatch& pp,         // 'outside' points
-            const labelList& adaptPatchIDs,     // patches forming 'outside'
-            const pointVectorField&,
+            polyMesh& mesh,
+            indirectPrimitivePatch& pp,         //!< 'outside' points
+            const labelList& adaptPatchIDs,     //!< patches forming 'outside'
+            const pointVectorField& displacement,
             const dictionary& paramDict
         );
-
 };
 
 
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
index e565582e86c..ee8a1137e4d 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
@@ -90,22 +90,6 @@ void Foam::motionSmootherAlgo::checkFld(const pointScalarField& fld)
 }
 
 
-Foam::labelHashSet Foam::motionSmootherAlgo::getPoints
-(
-    const labelUList& faceLabels
-) const
-{
-    labelHashSet usedPoints(mesh_.nPoints()/100);
-
-    for (const label faceId : faceLabels)
-    {
-        usedPoints.insert(mesh_.faces()[faceId]);
-    }
-
-    return usedPoints;
-}
-
-
 Foam::labelHashSet Foam::motionSmootherAlgo::getPoints
 (
     const labelHashSet& faceLabels
@@ -129,12 +113,12 @@ Foam::tmp<Foam::scalarField> Foam::motionSmootherAlgo::calcEdgeWeights
 {
     const edgeList& edges = mesh_.edges();
 
-    tmp<scalarField> twght(new scalarField(edges.size()));
-    scalarField& wght = twght.ref();
+    auto twght = tmp<scalarField>::New(edges.size());
+    auto& wght = twght.ref();
 
-    forAll(edges, edgeI)
+    forAll(edges, edgei)
     {
-        wght[edgeI] = 1.0/(edges[edgeI].mag(points)+SMALL);
+        wght[edgei] = 1.0/(edges[edgei].mag(points)+SMALL);
     }
     return twght;
 }
@@ -157,9 +141,8 @@ void Foam::motionSmootherAlgo::minSmooth
     );
     const pointScalarField& avgFld = tavgFld();
 
-    forAll(meshPoints, i)
+    for (const label pointi : meshPoints)
     {
-        label pointi = meshPoints[i];
         if (isAffectedPoint.test(pointi))
         {
             newFld[pointi] = min
@@ -193,7 +176,7 @@ void Foam::motionSmootherAlgo::minSmooth
 
     forAll(fld, pointi)
     {
-        if (isAffectedPoint.test(pointi) && isInternalPoint(pointi))
+        if (isAffectedPoint.test(pointi) && isInternalPoint_.test(pointi))
         {
             newFld[pointi] = min
             (
@@ -203,7 +186,7 @@ void Foam::motionSmootherAlgo::minSmooth
         }
     }
 
-   // Single and multi-patch constraints
+    // Single and multi-patch constraints
     pointConstraints::New(pMesh()).constrain(newFld, false);
 
 }
@@ -219,7 +202,7 @@ void Foam::motionSmootherAlgo::scaleField
 {
     for (const label pointi : pointLabels)
     {
-        if (isInternalPoint(pointi))
+        if (isInternalPoint_.test(pointi))
         {
             fld[pointi] *= scale;
         }
@@ -259,7 +242,7 @@ void Foam::motionSmootherAlgo::subtractField
 {
     for (const label pointi : pointLabels)
     {
-        if (isInternalPoint(pointi))
+        if (isInternalPoint_.test(pointi))
         {
             fld[pointi] = max(0.0, fld[pointi]-f);
         }
@@ -279,10 +262,8 @@ void Foam::motionSmootherAlgo::subtractField
     pointScalarField& fld
 ) const
 {
-    forAll(meshPoints, i)
+    for (const label pointi : meshPoints)
     {
-        label pointi = meshPoints[i];
-
         if (pointLabels.found(pointi))
         {
             fld[pointi] = max(0.0, fld[pointi]-f);
@@ -291,12 +272,6 @@ void Foam::motionSmootherAlgo::subtractField
 }
 
 
-bool Foam::motionSmootherAlgo::isInternalPoint(const label pointi) const
-{
-    return isInternalPoint_.test(pointi);
-}
-
-
 void Foam::motionSmootherAlgo::getAffectedFacesAndPoints
 (
     const label nPointIter,
@@ -306,8 +281,8 @@ void Foam::motionSmootherAlgo::getAffectedFacesAndPoints
     bitSet& isAffectedPoint
 ) const
 {
+    isAffectedPoint.reset();
     isAffectedPoint.resize(mesh_.nPoints());
-    isAffectedPoint = false;
 
     faceSet nbrFaces(mesh_, "checkFaces", wrongFaces);
 
@@ -318,17 +293,13 @@ void Foam::motionSmootherAlgo::getAffectedFacesAndPoints
 
     for (label i = 0; i < nPointIter; i++)
     {
-        pointSet nbrPoints(mesh_, "grownPoints", getPoints(nbrFaces.toc()));
+        pointSet nbrPoints(mesh_, "grownPoints", getPoints(nbrFaces));
 
         for (const label pointi : nbrPoints)
         {
-            const labelList& pCells = mesh_.pointCells(pointi);
-
-            forAll(pCells, pCelli)
+            for (const label celli : mesh_.pointCells(pointi))
             {
-                const cell& cFaces = mesh_.cells()[pCells[pCelli]];
-
-                nbrFaces.insert(cFaces);
+                nbrFaces.set(mesh_.cells()[celli]); // set multiple
             }
         }
         nbrFaces.sync(mesh_);
@@ -337,8 +308,7 @@ void Foam::motionSmootherAlgo::getAffectedFacesAndPoints
         {
             for (const label facei : nbrFaces)
             {
-                const face& f = mesh_.faces()[facei];
-                isAffectedPoint.set(f);
+                isAffectedPoint.set(mesh_.faces()[facei]);  // set multiple
             }
         }
     }
@@ -436,10 +406,8 @@ void Foam::motionSmootherAlgo::setDisplacementPatchFields
 
     // Adapt the fixedValue bc's (i.e. copy internal point data to
     // boundaryField for all affected patches)
-    forAll(patchIDs, i)
+    for (const label patchi : patchIDs)
     {
-        label patchi = patchIDs[i];
-
         displacementBf[patchi] ==
             displacementBf[patchi].patchInternalField();
     }
@@ -455,7 +423,7 @@ void Foam::motionSmootherAlgo::setDisplacementPatchFields
 
     forAll(patchSchedule, patchEvalI)
     {
-        label patchi = patchSchedule[patchEvalI].patch;
+        const label patchi = patchSchedule[patchEvalI].patch;
 
         if (!adaptPatchSet.found(patchi))
         {
@@ -478,12 +446,9 @@ void Foam::motionSmootherAlgo::setDisplacementPatchFields
     // Adapt the fixedValue bc's (i.e. copy internal point data to
     // boundaryField for all affected patches) to take the changes caused
     // by multi-corner constraints into account.
-    forAll(patchIDs, i)
+    for (const label patchi : patchIDs)
     {
-        label patchi = patchIDs[i];
-
-        displacementBf[patchi] ==
-            displacementBf[patchi].patchInternalField();
+        displacementBf[patchi] == displacementBf[patchi].patchInternalField();
     }
 }
 
@@ -523,14 +488,14 @@ void Foam::motionSmootherAlgo::setDisplacement
             mesh,
             isPatchPoint,
             maxEqOp<unsigned int>(),
-            0
+            0u
         );
-        forAll(cppMeshPoints, i)
+
+        for (const label pointi : cppMeshPoints)
         {
-            label pointI = cppMeshPoints[i];
-            if (isPatchPoint[pointI])
+            if (isPatchPoint.test(pointi))
             {
-                displacement[pointI] = Zero;
+                displacement[pointi] = Zero;
             }
         }
     }
@@ -611,7 +576,7 @@ void Foam::motionSmootherAlgo::correctBoundaryConditions
     // 1. evaluate on adaptPatches
     forAll(patchSchedule, patchEvalI)
     {
-        label patchi = patchSchedule[patchEvalI].patch;
+        const label patchi = patchSchedule[patchEvalI].patch;
 
         if (adaptPatchSet.found(patchi))
         {
@@ -632,7 +597,7 @@ void Foam::motionSmootherAlgo::correctBoundaryConditions
     // 2. evaluate on non-AdaptPatches
     forAll(patchSchedule, patchEvalI)
     {
-        label patchi = patchSchedule[patchEvalI].patch;
+        const label patchi = patchSchedule[patchEvalI].patch;
 
         if (!adaptPatchSet.found(patchi))
         {
@@ -686,9 +651,10 @@ void Foam::motionSmootherAlgo::modifyMotionPoints(pointField& newPoints) const
         const labelList& neIndices = twoDCorrector.normalEdgeIndices();
         const vector& pn = twoDCorrector.planeNormal();
 
+        for (const label edgei : neIndices)
         forAll(neIndices, i)
         {
-            const edge& e = edges[neIndices[i]];
+            const edge& e = edges[edgei];
 
             point& pStart = newPoints[e.start()];
 
@@ -964,7 +930,6 @@ bool Foam::motionSmootherAlgo::scaleMesh
         usedPoints.sync(mesh_);
 
 
-
         // Grow a few layers to determine
         // - points to be smoothed
         // - faces to be checked in next iteration
@@ -999,7 +964,7 @@ bool Foam::motionSmootherAlgo::scaleMesh
 
         scalarField eWeights(calcEdgeWeights(oldPoints_));
 
-        for (label i = 0; i < nSmoothScale; i++)
+        for (label i = 0; i < nSmoothScale; ++i)
         {
             if (adaptPatchIDs_.size())
             {
@@ -1051,10 +1016,8 @@ void Foam::motionSmootherAlgo::updateMesh()
     const pointBoundaryMesh& patches = pMesh_.boundary();
 
     // Check whether displacement has fixed value b.c. on adaptPatchID
-    forAll(adaptPatchIDs_, i)
+    for (const label patchi : adaptPatchIDs_)
     {
-        label patchi = adaptPatchIDs_[i];
-
         if
         (
            !isA<fixedValuePointPatchVectorField>
@@ -1078,12 +1041,7 @@ void Foam::motionSmootherAlgo::updateMesh()
     // Determine internal points. Note that for twoD there are no internal
     // points so we use the points of adaptPatchIDs instead
 
-    const labelList& meshPoints = pp_.meshPoints();
-
-    forAll(meshPoints, i)
-    {
-        isInternalPoint_.unset(meshPoints[i]);
-    }
+    isInternalPoint_.unset(pp_.meshPoints());  // unset multiple
 
     // Calculate master edge addressing
     isMasterEdge_ = syncTools::getMasterEdges(mesh_);
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H
index b6b85974fe5..742962e5bb0 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -104,7 +104,6 @@ class motionSmootherAlgo
         //  zero displacement.
         class maxMagEqOp
         {
-
         public:
 
             void operator()(vector& x, const vector& y) const
@@ -208,9 +207,6 @@ class motionSmootherAlgo
 
         static void checkFld(const pointScalarField&);
 
-        //- Get points used by given faces
-        labelHashSet getPoints(const labelUList& faceLabels) const;
-
         //- Get points used by given faces
         labelHashSet getPoints(const labelHashSet& faceLabels) const;
 
@@ -272,9 +268,6 @@ class motionSmootherAlgo
             pointScalarField&
         ) const;
 
-        //- Helper function. Is point internal?
-        bool isInternalPoint(const label pointi) const;
-
         //- Given a set of faces that cause smoothing and a number of
         //  iterations determine the maximum set of points who are affected
         //  and the accordingly affected faces.
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C
index 178574ee7bb..a228c61bb79 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -44,15 +44,15 @@ void Foam::motionSmootherAlgo::checkConstraints
 
     const polyBoundaryMesh& bm = mesh.boundaryMesh();
 
-    // first count the total number of patch-patch points
+    // First count the total number of patch-patch points
 
     label nPatchPatchPoints = 0;
 
-    forAll(bm, patchi)
+    for (const polyPatch& pp : bm)
     {
-        if (!isA<emptyPolyPatch>(bm[patchi]))
+        if (!isA<emptyPolyPatch>(pp))
         {
-            nPatchPatchPoints += bm[patchi].boundaryPoints().size();
+            nPatchPatchPoints += pp.boundaryPoints().size();
         }
     }
 
@@ -110,9 +110,9 @@ void Foam::motionSmootherAlgo::checkConstraints
             const labelList& bp = bm[patchi].boundaryPoints();
             const labelList& meshPoints = bm[patchi].meshPoints();
 
-            forAll(bp, pointi)
+            for (const label pointi : bp)
             {
-                label ppp = meshPoints[bp[pointi]];
+                const label ppp = meshPoints[pointi];
 
                 const Type& savedVal = boundaryPointValues[nPatchPatchPoints++];
 
@@ -141,24 +141,21 @@ Foam::motionSmootherAlgo::avg
     const scalarField& edgeWeight
 ) const
 {
-    tmp<GeometricField<Type, pointPatchField, pointMesh>> tres
+    auto tres = tmp<GeometricField<Type, pointPatchField, pointMesh>>::New
     (
-        new GeometricField<Type, pointPatchField, pointMesh>
+        IOobject
         (
-            IOobject
-            (
-                "avg("+fld.name()+')',
-                fld.time().timeName(),
-                fld.db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            fld.mesh(),
-            dimensioned<Type>(fld.dimensions(), Zero)
-        )
+            "avg("+fld.name()+')',
+            fld.time().timeName(),
+            fld.db(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        fld.mesh(),
+        dimensioned<Type>(fld.dimensions(), Zero)
     );
-    GeometricField<Type, pointPatchField, pointMesh>& res = tres.ref();
+    auto& res = tres.ref();
 
     const polyMesh& mesh = fld.mesh()();
 
@@ -169,23 +166,20 @@ Foam::motionSmootherAlgo::avg
     // Note: on coupled edges use only one edge (through isMasterEdge)
     // This is done so coupled edges do not get counted double.
 
-    scalarField sumWeight(mesh.nPoints(), 0.0);
+    scalarField sumWeight(mesh.nPoints(), Zero);
 
     const edgeList& edges = mesh.edges();
 
-    forAll(edges, edgeI)
+    for (const label edgei : isMasterEdge_)
     {
-        if (isMasterEdge_.get(edgeI) == 1)
-        {
-            const edge& e = edges[edgeI];
-            const scalar w = edgeWeight[edgeI];
+        const edge& e = edges[edgei];
+        const scalar w = edgeWeight[edgei];
 
-            res[e[0]] += w*fld[e[1]];
-            sumWeight[e[0]] += w;
+        res[e[0]] += w*fld[e[1]];
+        sumWeight[e[0]] += w;
 
-            res[e[1]] += w*fld[e[0]];
-            sumWeight[e[1]] += w;
-        }
+        res[e[1]] += w*fld[e[0]];
+        sumWeight[e[1]] += w;
     }
 
 
@@ -244,7 +238,7 @@ void Foam::motionSmootherAlgo::smooth
 
     forAll(fld, pointi)
     {
-        if (isInternalPoint(pointi))
+        if (isInternalPoint_.test(pointi))
         {
             newFld[pointi] = 0.5*fld[pointi] + 0.5*avgFld[pointi];
         }
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherData.H b/src/dynamicMesh/motionSmoother/motionSmootherData.H
index a4f2fa8a592..e3ea1661228 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherData.H
+++ b/src/dynamicMesh/motionSmoother/motionSmootherData.H
@@ -66,16 +66,10 @@ public:
     // Constructors
 
         //- Construct read
-        motionSmootherData
-        (
-            const pointMesh&
-        );
+        explicit motionSmootherData(const pointMesh& pMesh);
 
         //- Construct from pointDisplacement
-        motionSmootherData
-        (
-            const pointVectorField&
-        );
+        explicit motionSmootherData(const pointVectorField& displacement);
 
 
     // Member Functions
diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
index 46500f13999..c00eb53339e 100644
--- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
+++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
@@ -49,10 +49,8 @@ void Foam::polyMeshGeometry::updateFaceCentresAndAreas
 {
     const faceList& fs = mesh_.faces();
 
-    forAll(changedFaces, i)
+    for (const label facei : changedFaces)
     {
-        label facei = changedFaces[i];
-
         const labelList& f = fs[facei];
         label nPoints = f.size();
 
@@ -66,18 +64,18 @@ void Foam::polyMeshGeometry::updateFaceCentresAndAreas
         else
         {
             vector sumN = Zero;
-            scalar sumA = 0.0;
+            scalar sumA = Zero;
             vector sumAc = Zero;
 
             point fCentre = p[f[0]];
-            for (label pi = 1; pi < nPoints; pi++)
+            for (label pi = 1; pi < nPoints; ++pi)
             {
                 fCentre += p[f[pi]];
             }
 
             fCentre /= nPoints;
 
-            for (label pi = 0; pi < nPoints; pi++)
+            for (label pi = 0; pi < nPoints; ++pi)
             {
                 const point& nextPoint = p[f[(pi + 1) % nPoints]];
 
@@ -108,22 +106,21 @@ void Foam::polyMeshGeometry::updateCellCentresAndVols
 
     // Clear the fields for accumulation
     UIndirectList<vector>(cellCentres_, changedCells) = Zero;
-    UIndirectList<scalar>(cellVolumes_, changedCells) = 0.0;
+    UIndirectList<scalar>(cellVolumes_, changedCells) = Zero;
 
 
     // Re-calculate the changed cell centres and volumes
-    forAll(changedCells, changedCellI)
+    for (const label celli : changedCells)
     {
-        const label cellI(changedCells[changedCellI]);
-        const labelList& cFaces(cells[cellI]);
+        const labelList& cFaces = cells[celli];
 
         // Estimate the cell centre and bounding box using the face centres
         vector cEst(Zero);
         boundBox bb(boundBox::invertedBox);
 
-        forAll(cFaces, cFaceI)
+        for (const label facei : cFaces)
         {
-            const point& fc = faceCentres_[cFaces[cFaceI]];
+            const point& fc = faceCentres_[facei];
             cEst += fc;
             bb.add(fc);
         }
@@ -131,52 +128,50 @@ void Foam::polyMeshGeometry::updateCellCentresAndVols
 
 
         // Sum up the face-pyramid contributions
-        forAll(cFaces, cFaceI)
+        for (const label facei : cFaces)
         {
-            const label faceI(cFaces[cFaceI]);
-
             // Calculate 3* the face-pyramid volume
-            scalar pyr3Vol = faceAreas_[faceI] & (faceCentres_[faceI] - cEst);
+            scalar pyr3Vol = faceAreas_[facei] & (faceCentres_[facei] - cEst);
 
-            if (own[faceI] != cellI)
+            if (own[facei] != celli)
             {
                 pyr3Vol = -pyr3Vol;
             }
 
             // Accumulate face-pyramid volume
-            cellVolumes_[cellI] += pyr3Vol;
+            cellVolumes_[celli] += pyr3Vol;
 
             // Calculate the face-pyramid centre
-            const vector pCtr = (3.0/4.0)*faceCentres_[faceI] + (1.0/4.0)*cEst;
+            const vector pCtr = (3.0/4.0)*faceCentres_[facei] + (1.0/4.0)*cEst;
 
             // Accumulate volume-weighted face-pyramid centre
-            cellCentres_[cellI] += pyr3Vol*pCtr;
+            cellCentres_[celli] += pyr3Vol*pCtr;
         }
 
         // Average the accumulated quantities
 
-        if (mag(cellVolumes_[cellI]) > VSMALL)
+        if (mag(cellVolumes_[celli]) > VSMALL)
         {
-            point cc = cellCentres_[cellI] / cellVolumes_[cellI];
+            point cc = cellCentres_[celli] / cellVolumes_[celli];
 
             // Do additional check for collapsed cells since some volumes
             // (e.g. 1e-33) do not trigger above but do return completely
             // wrong cell centre
             if (bb.contains(cc))
             {
-                cellCentres_[cellI] = cc;
+                cellCentres_[celli] = cc;
             }
             else
             {
-                cellCentres_[cellI] = cEst;
+                cellCentres_[celli] = cEst;
             }
         }
         else
         {
-            cellCentres_[cellI] = cEst;
+            cellCentres_[celli] = cEst;
         }
 
-        cellVolumes_[cellI] *= (1.0/3.0);
+        cellVolumes_[celli] *= (1.0/3.0);
     }
 }
 
@@ -192,10 +187,8 @@ Foam::labelList Foam::polyMeshGeometry::affectedCells
 
     labelHashSet affectedCells(2*changedFaces.size());
 
-    forAll(changedFaces, i)
+    for (const label facei : changedFaces)
     {
-        label facei = changedFaces[i];
-
         affectedCells.insert(own[facei]);
 
         if (mesh.isInternalFace(facei))
@@ -245,7 +238,7 @@ Foam::scalar Foam::polyMeshGeometry::checkNonOrtho
                     << " deg." << endl;
             }
 
-            severeNonOrth++;
+            ++severeNonOrth;
         }
         else
         {
@@ -262,7 +255,7 @@ Foam::scalar Foam::polyMeshGeometry::checkNonOrtho
                     << " deg." << endl;
             }
 
-            errorNonOrth++;
+            ++errorNonOrth;
         }
 
         if (setPtr)
@@ -309,13 +302,11 @@ bool Foam::polyMeshGeometry::checkFaceTet
                     << ", const pointField&"
                     << ", const labelList&, labelHashSet*) : "
                     << "face " << facei
-                    << " has a triangle that points the wrong way."
-                     << endl
+                    << " has a triangle that points the wrong way." << nl
                     << "Tet quality: " << tetQual
                     << " Face " << facei
                     << endl;
             }
-
             if (setPtr)
             {
                 setPtr->insert(facei);
@@ -387,7 +378,7 @@ bool Foam::polyMeshGeometry::checkFaceDotProduct
     // Calculate coupled cell centre
     pointField neiCc(mesh.nBoundaryFaces());
 
-    for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
+    for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei)
     {
         neiCc[facei-mesh.nInternalFaces()] = cellCentres[own[facei]];
     }
@@ -403,10 +394,8 @@ bool Foam::polyMeshGeometry::checkFaceDotProduct
 
     label errorNonOrth = 0;
 
-    forAll(checkFaces, i)
+    for (const label facei : checkFaces)
     {
-        label facei = checkFaces[i];
-
         const point& ownCc = cellCentres[own[facei]];
 
         if (mesh.isInternalFace(facei))
@@ -431,11 +420,11 @@ bool Foam::polyMeshGeometry::checkFaceDotProduct
             }
 
             sumDDotS += dDotS;
-            nDDotS++;
+            ++nDDotS;
         }
         else
         {
-            label patchi = patches.whichPatch(facei);
+            const label patchi = patches.whichPatch(facei);
 
             if (patches[patchi].coupled())
             {
@@ -459,15 +448,15 @@ bool Foam::polyMeshGeometry::checkFaceDotProduct
                 }
 
                 sumDDotS += dDotS;
-                nDDotS++;
+                ++nDDotS;
             }
         }
     }
 
-    forAll(baffles, i)
+    for (const labelPair& baffle : baffles)
     {
-        label face0 = baffles[i].first();
-        label face1 = baffles[i].second();
+        const label face0 = baffle.first();
+        const label face1 = baffle.second();
 
         const point& ownCc = cellCentres[own[face0]];
 
@@ -491,7 +480,7 @@ bool Foam::polyMeshGeometry::checkFaceDotProduct
         }
 
         sumDDotS += dDotS;
-        nDDotS++;
+        ++nDDotS;
     }
 
     reduce(minDDotS, minOp<scalar>());
@@ -564,10 +553,8 @@ bool Foam::polyMeshGeometry::checkFacePyramids
 
     label nErrorPyrs = 0;
 
-    forAll(checkFaces, i)
+    for (const label facei : checkFaces)
     {
-        label facei = checkFaces[i];
-
         // Create the owner pyramid - it will have negative volume
         scalar pyrVol = pyramidPointFaceRef
         (
@@ -577,6 +564,8 @@ bool Foam::polyMeshGeometry::checkFacePyramids
 
         if (pyrVol > -minPyrVol)
         {
+            ++nErrorPyrs;
+
             if (report)
             {
                 Pout<< "bool polyMeshGeometry::checkFacePyramids("
@@ -590,14 +579,10 @@ bool Foam::polyMeshGeometry::checkFacePyramids
                     << mesh.cells()[own[facei]].labels(f)
                     << endl;
             }
-
-
             if (setPtr)
             {
                 setPtr->insert(facei);
             }
-
-            nErrorPyrs++;
         }
 
         if (mesh.isInternalFace(facei))
@@ -608,6 +593,8 @@ bool Foam::polyMeshGeometry::checkFacePyramids
 
             if (pyrVol < minPyrVol)
             {
+                ++nErrorPyrs;
+
                 if (report)
                 {
                     Pout<< "bool polyMeshGeometry::checkFacePyramids("
@@ -621,25 +608,22 @@ bool Foam::polyMeshGeometry::checkFacePyramids
                         << mesh.cells()[nei[facei]].labels(f)
                         << endl;
                 }
-
                 if (setPtr)
                 {
                     setPtr->insert(facei);
                 }
-
-                nErrorPyrs++;
             }
         }
     }
 
-    forAll(baffles, i)
+    for (const labelPair& baffle : baffles)
     {
-        label face0 = baffles[i].first();
-        label face1 = baffles[i].second();
+        const label face0 = baffle.first();
+        const label face1 = baffle.second();
 
         const point& ownCc = cellCentres[own[face0]];
 
-       // Create the owner pyramid - it will have negative volume
+        // Create the owner pyramid - it will have negative volume
         scalar pyrVolOwn = pyramidPointFaceRef
         (
             f[face0],
@@ -648,6 +632,8 @@ bool Foam::polyMeshGeometry::checkFacePyramids
 
         if (pyrVolOwn > -minPyrVol)
         {
+            ++nErrorPyrs;
+
             if (report)
             {
                 Pout<< "bool polyMeshGeometry::checkFacePyramids("
@@ -661,14 +647,10 @@ bool Foam::polyMeshGeometry::checkFacePyramids
                     << mesh.cells()[own[face0]].labels(f)
                     << endl;
             }
-
-
             if (setPtr)
             {
                 setPtr->insert(face0);
             }
-
-            nErrorPyrs++;
         }
 
         // Create the neighbour pyramid - it will have positive volume
@@ -677,6 +659,8 @@ bool Foam::polyMeshGeometry::checkFacePyramids
 
         if (pyrVolNbr < minPyrVol)
         {
+            ++nErrorPyrs;
+
             if (report)
             {
                 Pout<< "bool polyMeshGeometry::checkFacePyramids("
@@ -690,13 +674,10 @@ bool Foam::polyMeshGeometry::checkFacePyramids
                     << mesh.cells()[own[face1]].labels(f)
                     << endl;
             }
-
             if (setPtr)
             {
                 setPtr->insert(face0);
             }
-
-            nErrorPyrs++;
         }
     }
 
@@ -747,7 +728,7 @@ bool Foam::polyMeshGeometry::checkFaceTets
     // Calculate coupled cell centre
     pointField neiCc(mesh.nBoundaryFaces());
 
-    for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
+    for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei)
     {
         neiCc[facei - mesh.nInternalFaces()] = cellCentres[own[facei]];
     }
@@ -756,10 +737,8 @@ bool Foam::polyMeshGeometry::checkFaceTets
 
     label nErrorTets = 0;
 
-    forAll(checkFaces, i)
+    for (const label facei : checkFaces)
     {
-        label facei = checkFaces[i];
-
         // Create the owner pyramid - note: exchange cell and face centre
         // to get positive volume.
         bool tetError = checkFaceTet
@@ -776,7 +755,7 @@ bool Foam::polyMeshGeometry::checkFaceTets
 
         if (tetError)
         {
-            nErrorTets++;
+            ++nErrorTets;
         }
 
         if (mesh.isInternalFace(facei))
@@ -796,7 +775,7 @@ bool Foam::polyMeshGeometry::checkFaceTets
 
             if (tetError)
             {
-                nErrorTets++;
+                ++nErrorTets;
             }
 
             if
@@ -810,12 +789,11 @@ bool Foam::polyMeshGeometry::checkFaceTets
                 ) == -1
             )
             {
+                ++nErrorTets;
                 if (setPtr)
                 {
                     setPtr->insert(facei);
                 }
-
-                nErrorTets++;
             }
         }
         else
@@ -836,12 +814,11 @@ bool Foam::polyMeshGeometry::checkFaceTets
                     ) == -1
                 )
                 {
+                    ++nErrorTets;
                     if (setPtr)
                     {
                         setPtr->insert(facei);
                     }
-
-                    nErrorTets++;
                 }
             }
             else
@@ -857,21 +834,20 @@ bool Foam::polyMeshGeometry::checkFaceTets
                     ) == -1
                 )
                 {
+                    ++nErrorTets;
                     if (setPtr)
                     {
                         setPtr->insert(facei);
                     }
-
-                    nErrorTets++;
                 }
             }
         }
     }
 
-    forAll(baffles, i)
+    for (const labelPair& baffle : baffles)
     {
-        label face0 = baffles[i].first();
-        label face1 = baffles[i].second();
+        const label face0 = baffle.first();
+        const label face1 = baffle.second();
 
         bool tetError = checkFaceTet
         (
@@ -887,7 +863,7 @@ bool Foam::polyMeshGeometry::checkFaceTets
 
         if (tetError)
         {
-            nErrorTets++;
+            ++nErrorTets;
         }
 
         // Create the neighbour tets - they will have positive volume
@@ -905,7 +881,7 @@ bool Foam::polyMeshGeometry::checkFaceTets
 
         if (tetError)
         {
-            nErrorTets++;
+            ++nErrorTets;
         }
 
         if
@@ -920,12 +896,11 @@ bool Foam::polyMeshGeometry::checkFaceTets
             ) == -1
         )
         {
+            ++nErrorTets;
             if (setPtr)
             {
                 setPtr->insert(face0);
             }
-
-            nErrorTets++;
         }
     }
 
@@ -984,10 +959,8 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
 
     label nWarnSkew = 0;
 
-    forAll(checkFaces, i)
+    for (const label facei : checkFaces)
     {
-        label facei = checkFaces[i];
-
         if (mesh.isInternalFace(facei))
         {
             scalar skewness = primitiveMeshTools::faceSkewness
@@ -1007,18 +980,17 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
             // mesh.
             if (skewness > internalSkew)
             {
+                ++nWarnSkew;
+
                 if (report)
                 {
                     Pout<< "Severe skewness for face " << facei
                         << " skewness = " << skewness << endl;
                 }
-
                 if (setPtr)
                 {
                     setPtr->insert(facei);
                 }
-
-                nWarnSkew++;
             }
 
             maxSkew = max(maxSkew, skewness);
@@ -1042,18 +1014,17 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
             // mesh.
             if (skewness > internalSkew)
             {
+                ++nWarnSkew;
+
                 if (report)
                 {
                     Pout<< "Severe skewness for coupled face " << facei
                         << " skewness = " << skewness << endl;
                 }
-
                 if (setPtr)
                 {
                     setPtr->insert(facei);
                 }
-
-                nWarnSkew++;
             }
 
             maxSkew = max(maxSkew, skewness);
@@ -1077,28 +1048,27 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
             // mesh.
             if (skewness > boundarySkew)
             {
+                ++nWarnSkew;
+
                 if (report)
                 {
                     Pout<< "Severe skewness for boundary face " << facei
                         << " skewness = " << skewness << endl;
                 }
-
                 if (setPtr)
                 {
                     setPtr->insert(facei);
                 }
-
-                nWarnSkew++;
             }
 
             maxSkew = max(maxSkew, skewness);
         }
     }
 
-    forAll(baffles, i)
+    for (const labelPair& baffle : baffles)
     {
-        label face0 = baffles[i].first();
-        label face1 = baffles[i].second();
+        const label face0 = baffle.first();
+        const label face1 = baffle.second();
 
         const point& ownCc = cellCentres[own[face0]];
         const point& neiCc = cellCentres[own[face1]];
@@ -1120,18 +1090,17 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
         // mesh.
         if (skewness > internalSkew)
         {
+            ++nWarnSkew;
+
             if (report)
             {
                 Pout<< "Severe skewness for face " << face0
                     << " skewness = " << skewness << endl;
             }
-
             if (setPtr)
             {
                 setPtr->insert(face0);
             }
-
-            nWarnSkew++;
         }
 
         maxSkew = max(maxSkew, skewness);
@@ -1189,7 +1158,7 @@ bool Foam::polyMeshGeometry::checkFaceWeights
     // Calculate coupled cell centre
     pointField neiCc(mesh.nBoundaryFaces());
 
-    for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
+    for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei)
     {
         neiCc[facei-mesh.nInternalFaces()] = cellCentres[own[facei]];
     }
@@ -1200,10 +1169,8 @@ bool Foam::polyMeshGeometry::checkFaceWeights
 
     label nWarnWeight = 0;
 
-    forAll(checkFaces, i)
+    for (const label facei : checkFaces)
     {
-        label facei = checkFaces[i];
-
         const point& fc = faceCentres[facei];
         const vector& fa = faceAreas[facei];
 
@@ -1216,18 +1183,17 @@ bool Foam::polyMeshGeometry::checkFaceWeights
 
             if (weight < warnWeight)
             {
+                ++nWarnWeight;
+
                 if (report)
                 {
                     Pout<< "Small weighting factor for face " << facei
                         << " weight = " << weight << endl;
                 }
-
                 if (setPtr)
                 {
                     setPtr->insert(facei);
                 }
-
-                nWarnWeight++;
             }
 
             minWeight = min(minWeight, weight);
@@ -1243,18 +1209,17 @@ bool Foam::polyMeshGeometry::checkFaceWeights
 
                 if (weight < warnWeight)
                 {
+                    ++nWarnWeight;
+
                     if (report)
                     {
                         Pout<< "Small weighting factor for face " << facei
                             << " weight = " << weight << endl;
                     }
-
                     if (setPtr)
                     {
                         setPtr->insert(facei);
                     }
-
-                    nWarnWeight++;
                 }
 
                 minWeight = min(minWeight, weight);
@@ -1262,10 +1227,10 @@ bool Foam::polyMeshGeometry::checkFaceWeights
         }
     }
 
-    forAll(baffles, i)
+    for (const labelPair& baffle : baffles)
     {
-        label face0 = baffles[i].first();
-        label face1 = baffles[i].second();
+        const label face0 = baffle.first();
+        const label face1 = baffle.second();
 
         const point& ownCc = cellCentres[own[face0]];
         const point& fc = faceCentres[face0];
@@ -1277,18 +1242,17 @@ bool Foam::polyMeshGeometry::checkFaceWeights
 
         if (weight < warnWeight)
         {
+            ++nWarnWeight;
+
             if (report)
             {
                 Pout<< "Small weighting factor for face " << face0
                     << " weight = " << weight << endl;
             }
-
             if (setPtr)
             {
                 setPtr->insert(face0);
             }
-
-            nWarnWeight++;
         }
 
         minWeight = min(minWeight, weight);
@@ -1342,7 +1306,7 @@ bool Foam::polyMeshGeometry::checkVolRatio
     // Calculate coupled cell vol
     scalarField neiVols(mesh.nBoundaryFaces());
 
-    for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
+    for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei)
     {
         neiVols[facei-mesh.nInternalFaces()] = cellVolumes[own[facei]];
     }
@@ -1353,10 +1317,8 @@ bool Foam::polyMeshGeometry::checkVolRatio
 
     label nWarnRatio = 0;
 
-    forAll(checkFaces, i)
+    for (const label facei : checkFaces)
     {
-        label facei = checkFaces[i];
-
         scalar ownVol = mag(cellVolumes[own[facei]]);
 
         scalar neiVol = -GREAT;
@@ -1381,28 +1343,27 @@ bool Foam::polyMeshGeometry::checkVolRatio
 
             if (ratio < warnRatio)
             {
+                ++nWarnRatio;
+
                 if (report)
                 {
                     Pout<< "Small ratio for face " << facei
                         << " ratio = " << ratio << endl;
                 }
-
                 if (setPtr)
                 {
                     setPtr->insert(facei);
                 }
-
-                nWarnRatio++;
             }
 
             minRatio = min(minRatio, ratio);
         }
     }
 
-    forAll(baffles, i)
+    for (const labelPair& baffle : baffles)
     {
-        label face0 = baffles[i].first();
-        label face1 = baffles[i].second();
+        const label face0 = baffle.first();
+        const label face1 = baffle.second();
 
         scalar ownVol = mag(cellVolumes[own[face0]]);
 
@@ -1414,18 +1375,17 @@ bool Foam::polyMeshGeometry::checkVolRatio
 
             if (ratio < warnRatio)
             {
+                ++nWarnRatio;
+
                 if (report)
                 {
                     Pout<< "Small ratio for face " << face0
                         << " ratio = " << ratio << endl;
                 }
-
                 if (setPtr)
                 {
                     setPtr->insert(face0);
                 }
-
-                nWarnRatio++;
             }
 
             minRatio = min(minRatio, ratio);
@@ -1492,10 +1452,8 @@ bool Foam::polyMeshGeometry::checkFaceAngles
 
     label errorFacei = -1;
 
-    forAll(checkFaces, i)
+    for (const label facei : checkFaces)
     {
-        label facei = checkFaces[i];
-
         const face& f = fcs[facei];
 
         const vector faceNormal = normalised(faceAreas[facei]);
@@ -1535,15 +1493,15 @@ bool Foam::polyMeshGeometry::checkFaceAngles
                         {
                             // Count only one error per face.
                             errorFacei = facei;
-                            nConcave++;
+                            ++nConcave;
                         }
 
+                        maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
+
                         if (setPtr)
                         {
                             setPtr->insert(facei);
                         }
-
-                        maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
                     }
                 }
             }
@@ -1619,53 +1577,8 @@ bool Foam::polyMeshGeometry::checkFaceTwist
 
     const faceList& fcs = mesh.faces();
 
-
     label nWarped = 0;
 
-//    forAll(checkFaces, i)
-//    {
-//        label facei = checkFaces[i];
-//
-//        const face& f = fcs[facei];
-//
-//        scalar magArea = mag(faceAreas[facei]);
-//
-//        if (f.size() > 3 && magArea > VSMALL)
-//        {
-//            const vector nf = faceAreas[facei] / magArea;
-//
-//            const point& fc = faceCentres[facei];
-//
-//            forAll(f, fpI)
-//            {
-//                vector triArea
-//                (
-//                    triPointRef
-//                    (
-//                        p[f[fpI]],
-//                        p[f.nextLabel(fpI)],
-//                        fc
-//                    ).areaNormal()
-//                );
-//
-//                scalar magTri = mag(triArea);
-//
-//                if (magTri > VSMALL && ((nf & triArea/magTri) < minTwist))
-//                {
-//                    nWarped++;
-//
-//                    if (setPtr)
-//                    {
-//                        setPtr->insert(facei);
-//                    }
-//
-//                    break;
-//                }
-//            }
-//        }
-//    }
-
-
     const labelList& own = mesh.faceOwner();
     const labelList& nei = mesh.faceNeighbour();
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
@@ -1673,16 +1586,14 @@ bool Foam::polyMeshGeometry::checkFaceTwist
     // Calculate coupled cell centre
     pointField neiCc(mesh.nBoundaryFaces());
 
-    for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
+    for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei)
     {
         neiCc[facei-mesh.nInternalFaces()] = cellCentres[own[facei]];
     }
     syncTools::swapBoundaryFacePositions(mesh, neiCc);
 
-    forAll(checkFaces, i)
+    for (const label facei : checkFaces)
     {
-        label facei = checkFaces[i];
-
         const face& f = fcs[facei];
 
         if (f.size() > 3)
@@ -1737,7 +1648,7 @@ bool Foam::polyMeshGeometry::checkFaceTwist
 
                     if (magTri > VSMALL && ((nf & triArea/magTri) < minTwist))
                     {
-                        nWarped++;
+                        ++nWarped;
 
                         if (setPtr)
                         {
@@ -1814,10 +1725,8 @@ bool Foam::polyMeshGeometry::checkTriangleTwist
 
     label nWarped = 0;
 
-    forAll(checkFaces, i)
+    for (const label facei : checkFaces)
     {
-        label facei = checkFaces[i];
-
         const face& f = fcs[facei];
 
         if (f.size() > 3)
@@ -1872,7 +1781,7 @@ bool Foam::polyMeshGeometry::checkTriangleTwist
 
                         if ((prevN & triN) < minTwist)
                         {
-                            nWarped++;
+                            ++nWarped;
 
                             if (setPtr)
                             {
@@ -1886,7 +1795,7 @@ bool Foam::polyMeshGeometry::checkTriangleTwist
                     }
                     else if (minTwist > 0)
                     {
-                        nWarped++;
+                        ++nWarped;
 
                         if (setPtr)
                         {
@@ -1964,10 +1873,8 @@ bool Foam::polyMeshGeometry::checkFaceFlatness
 
     label nWarped = 0;
 
-    forAll(checkFaces, i)
+    for (const label facei : checkFaces)
     {
-        label facei = checkFaces[i];
-
         const face& f = fcs[facei];
 
         if (f.size() > 3)
@@ -1989,8 +1896,7 @@ bool Foam::polyMeshGeometry::checkFaceFlatness
 
             if (sumArea/mag(faceAreas[facei]) < minFlatness)
             {
-                nWarped++;
-
+                ++nWarped;
                 if (setPtr)
                 {
                     setPtr->insert(facei);
@@ -2051,17 +1957,15 @@ bool Foam::polyMeshGeometry::checkFaceArea
 {
     label nZeroArea = 0;
 
-    forAll(checkFaces, i)
+    for (const label facei : checkFaces)
     {
-        label facei = checkFaces[i];
-
         if (mag(faceAreas[facei]) < minArea)
         {
+            ++nZeroArea;
             if (setPtr)
             {
                 setPtr->insert(facei);
             }
-            nZeroArea++;
         }
     }
 
@@ -2118,17 +2022,15 @@ bool Foam::polyMeshGeometry::checkCellDeterminant
     label nSumDet = 0;
     label nWarnDet = 0;
 
-    forAll(affectedCells, i)
+    for (const label celli : affectedCells)
     {
-        const cell& cFaces = cells[affectedCells[i]];
+        const cell& cFaces = cells[celli];
 
         tensor areaSum(Zero);
         scalar magAreaSum = 0;
 
-        forAll(cFaces, cFacei)
+        for (const label facei : cFaces)
         {
-            label facei = cFaces[cFacei];
-
             scalar magArea = mag(faceAreas[facei]);
 
             magAreaSum += magArea;
@@ -2139,20 +2041,15 @@ bool Foam::polyMeshGeometry::checkCellDeterminant
 
         minDet = min(minDet, scaledDet);
         sumDet += scaledDet;
-        nSumDet++;
+        ++nSumDet;
 
         if (scaledDet < warnDet)
         {
+            ++nWarnDet;
             if (setPtr)
             {
-                // Insert all faces of the cell.
-                forAll(cFaces, cFacei)
-                {
-                    label facei = cFaces[cFacei];
-                    setPtr->insert(facei);
-                }
+                setPtr->insert(cFaces);  // Insert all faces of the cell
             }
-            nWarnDet++;
         }
     }
 
-- 
GitLab