From ad2baed5af13cd495864d80b0d628d5f162c07da Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Fri, 16 Nov 2018 10:30:02 +0100
Subject: [PATCH] ENH update use of bitSet in dynamicRefineFvMesh #1075

---
 .../dynamicRefineFvMesh/dynamicRefineFvMesh.C | 473 +++++++-----------
 .../dynamicRefineFvMesh/dynamicRefineFvMesh.H |  28 +-
 2 files changed, 206 insertions(+), 295 deletions(-)

diff --git a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
index 8a85e3d25d..abc72e686f 100644
--- a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
+++ b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
@@ -33,6 +33,7 @@ License
 #include "pointFields.H"
 #include "sigFpe.H"
 #include "cellSet.H"
+#include "HashOps.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -44,32 +45,6 @@ namespace Foam
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-Foam::label Foam::dynamicRefineFvMesh::count
-(
-    const bitSet& l,
-    const unsigned int val
-)
-{
-    label n = 0;
-    forAll(l, i)
-    {
-        if (l.get(i) == val)
-        {
-            n++;
-        }
-
-        // debug also serves to get-around Clang compiler trying to optimsie
-        // out this forAll loop under O3 optimisation
-        if (debug & 128)
-        {
-            Info<< "n=" << n << endl;
-        }
-    }
-
-    return n;
-}
-
-
 void Foam::dynamicRefineFvMesh::calculateProtectedCells
 (
     bitSet& unrefineableCell
@@ -88,81 +63,89 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
     // Get neighbouring cell level
     labelList neiLevel(nFaces()-nInternalFaces());
 
-    for (label facei = nInternalFaces(); facei < nFaces(); facei++)
+    for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
     {
         neiLevel[facei-nInternalFaces()] = cellLevel[faceOwner()[facei]];
     }
     syncTools::swapBoundaryFaceList(*this, neiLevel);
 
 
+    bitSet seedFace;
+
     while (true)
     {
         // Pick up faces on border of protected cells
-        boolList seedFace(nFaces(), false);
+        seedFace.reset();
+        seedFace.resize(nFaces());
 
-        forAll(faceNeighbour(), facei)
+        for (label facei = 0; facei < nInternalFaces(); ++facei)
         {
-            label own = faceOwner()[facei];
-            bool ownProtected = unrefineableCell.get(own);
-            label nei = faceNeighbour()[facei];
-            bool neiProtected = unrefineableCell.get(nei);
+            const label own = faceOwner()[facei];
+            const label nei = faceNeighbour()[facei];
 
-            if (ownProtected && (cellLevel[nei] > cellLevel[own]))
+            if
+            (
+                // Protected owner
+                (
+                    unrefineableCell.test(own)
+                 && (cellLevel[nei] > cellLevel[own])
+                )
+             ||
+                // Protected neighbour
+                (
+                    unrefineableCell.test(nei)
+                 && (cellLevel[own] > cellLevel[nei])
+                )
+            )
             {
-                seedFace[facei] = true;
-            }
-            else if (neiProtected && (cellLevel[own] > cellLevel[nei]))
-            {
-                seedFace[facei] = true;
+                seedFace.set(facei);
             }
         }
         for (label facei = nInternalFaces(); facei < nFaces(); facei++)
         {
-            label own = faceOwner()[facei];
-            bool ownProtected = unrefineableCell.get(own);
+            const label own = faceOwner()[facei];
+
             if
             (
-                ownProtected
-             && (neiLevel[facei-nInternalFaces()] > cellLevel[own])
+                // Protected owner
+                (
+                    unrefineableCell.test(own)
+                 && (neiLevel[facei-nInternalFaces()] > cellLevel[own])
+                )
             )
             {
-                seedFace[facei] = true;
+                seedFace.set(facei);
             }
         }
 
-        syncTools::syncFaceList(*this, seedFace, orEqOp<bool>());
+        syncTools::syncFaceList(*this, seedFace, orEqOp<unsigned int>());
 
 
         // Extend unrefineableCell
         bool hasExtended = false;
 
-        for (label facei = 0; facei < nInternalFaces(); facei++)
+        for (label facei = 0; facei < nInternalFaces(); ++facei)
         {
-            if (seedFace[facei])
+            if (seedFace.test(facei))
             {
-                label own = faceOwner()[facei];
-                if (unrefineableCell.get(own) == 0)
+                if (unrefineableCell.set(faceOwner()[facei]))
                 {
-                    unrefineableCell.set(own, 1);
                     hasExtended = true;
                 }
-
-                label nei = faceNeighbour()[facei];
-                if (unrefineableCell.get(nei) == 0)
+                if (unrefineableCell.set(faceNeighbour()[facei]))
                 {
-                    unrefineableCell.set(nei, 1);
                     hasExtended = true;
                 }
             }
         }
-        for (label facei = nInternalFaces(); facei < nFaces(); facei++)
+        for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
         {
-            if (seedFace[facei])
+            if (seedFace.test(facei))
             {
-                label own = faceOwner()[facei];
-                if (unrefineableCell.get(own) == 0)
+                const label own = faceOwner()[facei];
+
+                if (unrefineableCell.set(own))
                 {
-                    unrefineableCell.set(own, 1);
                     hasExtended = true;
                 }
             }
@@ -198,9 +181,9 @@ void Foam::dynamicRefineFvMesh::readDict()
 
     // Rework into hashtable.
     correctFluxes_.resize(fluxVelocities.size());
-    forAll(fluxVelocities, i)
+    for (const auto& pr : fluxVelocities)
     {
-        correctFluxes_.insert(fluxVelocities[i][0], fluxVelocities[i][1]);
+        correctFluxes_.insert(pr.first(), pr.second());
     }
 
     refineDict.readEntry("dumpLevel", dumpLevel_);
@@ -210,6 +193,7 @@ void Foam::dynamicRefineFvMesh::readDict()
 void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
 {
     dynamicFvMesh::mapFields(mpm);
+
     // Correct the flux for modified/added faces. All the faces which only
     // have been renumbered will already have been handled by the mapping.
     {
@@ -220,46 +204,42 @@ void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
         // on the coarse cell that get split into four (or rather the
         // master face gets modified and three faces get added from the master)
         // Estimate number of faces created
-        labelHashSet masterFaces
-        (
-            max
-            (
-                mag(nFaces()-mpm.nOldFaces())/4,
-                nFaces()/100
-            )
-        );
+
+        bitSet masterFaces(nFaces());
 
         forAll(faceMap, facei)
         {
-            label oldFacei = faceMap[facei];
+            const label oldFacei = faceMap[facei];
 
             if (oldFacei >= 0)
             {
-                label masterFacei = reverseFaceMap[oldFacei];
+                const label masterFacei = reverseFaceMap[oldFacei];
 
                 if (masterFacei < 0)
                 {
                     FatalErrorInFunction
                         << "Problem: should not have removed faces"
                         << " when refining."
-                        << nl << "face:" << facei << abort(FatalError);
+                        << nl << "face:" << facei << endl
+                        << abort(FatalError);
                 }
                 else if (masterFacei != facei)
                 {
-                    masterFaces.insert(masterFacei);
+                    masterFaces.set(masterFacei);
                 }
             }
         }
+
         if (debug)
         {
-            Pout<< "Found " << masterFaces.size() << " split faces " << endl;
+            Pout<< "Found " << masterFaces.count() << " split faces " << endl;
         }
 
         HashTable<surfaceScalarField*> fluxes
         (
             lookupClass<surfaceScalarField>()
         );
-        forAllIter(HashTable<surfaceScalarField*>, fluxes, iter)
+        forAllIters(fluxes, iter)
         {
             if (!correctFluxes_.found(iter.key()))
             {
@@ -282,13 +262,13 @@ void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
                 continue;
             }
 
+            surfaceScalarField& phi = *iter();
+
             if (UName == "NaN")
             {
                 Pout<< "Setting surfaceScalarField " << iter.key()
                     << " to NaN" << endl;
 
-                surfaceScalarField& phi = *iter();
-
                 sigFpe::fillNan(phi.primitiveFieldRef());
 
                 continue;
@@ -301,7 +281,6 @@ void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
                     << endl;
             }
 
-            surfaceScalarField& phi = *iter();
             const surfaceScalarField phiU
             (
                 fvc::interpolate
@@ -312,9 +291,9 @@ void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
             );
 
             // Recalculate new internal faces.
-            for (label facei = 0; facei < nInternalFaces(); facei++)
+            for (label facei = 0; facei < nInternalFaces(); ++facei)
             {
-                label oldFacei = faceMap[facei];
+                const label oldFacei = faceMap[facei];
 
                 if (oldFacei == -1)
                 {
@@ -329,8 +308,8 @@ void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
             }
 
             // Recalculate new boundary faces.
-            surfaceScalarField::Boundary& phiBf =
-                phi.boundaryFieldRef();
+            surfaceScalarField::Boundary& phiBf = phi.boundaryFieldRef();
+
             forAll(phiBf, patchi)
             {
                 fvsPatchScalarField& patchPhi = phiBf[patchi];
@@ -341,7 +320,7 @@ void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
 
                 forAll(patchPhi, i)
                 {
-                    label oldFacei = faceMap[facei];
+                    const label oldFacei = faceMap[facei];
 
                     if (oldFacei == -1)
                     {
@@ -354,23 +333,21 @@ void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
                         patchPhi[i] = patchPhiU[i];
                     }
 
-                    facei++;
+                    ++facei;
                 }
             }
 
             // Update master faces
-            forAllConstIter(labelHashSet, masterFaces, iter)
+            for (const label facei : masterFaces)
             {
-                label facei = iter.key();
-
                 if (isInternalFace(facei))
                 {
                     phi[facei] = phiU[facei];
                 }
                 else
                 {
-                    label patchi = boundaryMesh().whichPatch(facei);
-                    label i = facei - boundaryMesh()[patchi].start();
+                    const label patchi = boundaryMesh().whichPatch(facei);
+                    const label i = facei - boundaryMesh()[patchi].start();
 
                     const fvsPatchScalarField& patchPhiU =
                         phiU.boundaryField()[patchi];
@@ -424,9 +401,9 @@ Foam::dynamicRefineFvMesh::refine
     if (debug)
     {
         // Check map.
-        for (label facei = 0; facei < nInternalFaces(); facei++)
+        for (label facei = 0; facei < nInternalFaces(); ++facei)
         {
-            label oldFacei = map().faceMap()[facei];
+            const label oldFacei = map().faceMap()[facei];
 
             if (oldFacei >= nInternalFaces())
             {
@@ -474,14 +451,17 @@ Foam::dynamicRefineFvMesh::refine
 
         forAll(newProtectedCell, celli)
         {
-            label oldCelli = map().cellMap()[celli];
-            newProtectedCell.set(celli, protectedCell_.get(oldCelli));
+            const label oldCelli = map().cellMap()[celli];
+            if (protectedCell_.test(oldCelli))
+            {
+                newProtectedCell.set(celli);
+            }
         }
         protectedCell_.transfer(newProtectedCell);
     }
 
     // Debug: Check refinement levels (across faces only)
-    meshCutter_.checkRefinementLevels(-1, labelList(0));
+    meshCutter_.checkRefinementLevels(-1, labelList());
 
     return map;
 }
@@ -508,21 +488,19 @@ Foam::dynamicRefineFvMesh::unrefine
     Map<label> faceToSplitPoint(3*splitPoints.size());
 
     {
-        forAll(splitPoints, i)
+        for (const label pointi : splitPoints)
         {
-            label pointi = splitPoints[i];
-
             const labelList& pEdges = pointEdges()[pointi];
 
-            forAll(pEdges, j)
+            for (const label edgei : pEdges)
             {
-                label otherPointi = edges()[pEdges[j]].otherVertex(pointi);
+                const label otherPointi = edges()[edgei].otherVertex(pointi);
 
                 const labelList& pFaces = pointFaces()[otherPointi];
 
-                forAll(pFaces, pFacei)
+                for (const label facei : pFaces)
                 {
-                    faceToSplitPoint.insert(pFaces[pFacei], otherPointi);
+                    faceToSplitPoint.insert(facei, otherPointi);
                 }
             }
         }
@@ -565,7 +543,7 @@ Foam::dynamicRefineFvMesh::unrefine
         (
             lookupClass<surfaceScalarField>()
         );
-        forAllIter(HashTable<surfaceScalarField*>, fluxes, iter)
+        forAllIters(fluxes, iter)
         {
             if (!correctFluxes_.found(iter.key()))
             {
@@ -588,12 +566,11 @@ Foam::dynamicRefineFvMesh::unrefine
                 continue;
             }
 
-            if (debug)
-            {
-                Info<< "Mapping flux " << iter.key()
-                    << " using interpolated flux " << UName
-                    << endl;
-            }
+            DebugInfo
+                << "Mapping flux " << iter.key()
+                << " using interpolated flux " << UName
+                << endl;
+
 
             surfaceScalarField& phi = *iter();
             surfaceScalarField::Boundary& phiBf =
@@ -609,15 +586,15 @@ Foam::dynamicRefineFvMesh::unrefine
             );
 
 
-            forAllConstIter(Map<label>, faceToSplitPoint, iter)
+            forAllConstIters(faceToSplitPoint, iter)
             {
-                label oldFacei = iter.key();
-                label oldPointi = iter();
+                const label oldFacei = iter.key();
+                const label oldPointi = iter.object();
 
                 if (reversePointMap[oldPointi] < 0)
                 {
                     // midpoint was removed. See if face still exists.
-                    label facei = reverseFaceMap[oldFacei];
+                    const label facei = reverseFaceMap[oldFacei];
 
                     if (facei >= 0)
                     {
@@ -652,17 +629,17 @@ Foam::dynamicRefineFvMesh::unrefine
 
         forAll(newProtectedCell, celli)
         {
-            label oldCelli = map().cellMap()[celli];
-            if (oldCelli >= 0)
+            const label oldCelli = map().cellMap()[celli];
+            if (protectedCell_.test(oldCelli))
             {
-                newProtectedCell.set(celli, protectedCell_.get(oldCelli));
+                newProtectedCell.set(celli);
             }
         }
         protectedCell_.transfer(newProtectedCell);
     }
 
     // Debug: Check refinement levels (across faces only)
-    meshCutter_.checkRefinementLevels(-1, labelList(0));
+    meshCutter_.checkRefinementLevels(-1, labelList());
 
     return map;
 }
@@ -677,9 +654,9 @@ Foam::dynamicRefineFvMesh::maxPointField(const scalarField& pFld) const
     {
         const labelList& pCells = pointCells()[pointi];
 
-        forAll(pCells, i)
+        for (const label celli : pCells)
         {
-            vFld[pCells[i]] = max(vFld[pCells[i]], pFld[pointi]);
+            vFld[celli] = max(vFld[celli], pFld[pointi]);
         }
     }
     return vFld;
@@ -695,9 +672,9 @@ Foam::dynamicRefineFvMesh::maxCellField(const volScalarField& vFld) const
     {
         const labelList& pCells = pointCells()[pointi];
 
-        forAll(pCells, i)
+        for (const label celli : pCells)
         {
-            pFld[pointi] = max(pFld[pointi], vFld[pCells[i]]);
+            pFld[pointi] = max(pFld[pointi], vFld[celli]);
         }
     }
     return pFld;
@@ -714,9 +691,9 @@ Foam::dynamicRefineFvMesh::cellToPoint(const scalarField& vFld) const
         const labelList& pCells = pointCells()[pointi];
 
         scalar sum = 0.0;
-        forAll(pCells, i)
+        for (const label celli : pCells)
         {
-            sum += vFld[pCells[i]];
+            sum += vFld[celli];
         }
         pFld[pointi] = sum/pCells.size();
     }
@@ -731,7 +708,7 @@ Foam::scalarField Foam::dynamicRefineFvMesh::error
     const scalar maxLevel
 ) const
 {
-    scalarField c(fld.size(), -1);
+    scalarField c(fld.size(), scalar(-1));
 
     forAll(fld, i)
     {
@@ -774,7 +751,7 @@ void Foam::dynamicRefineFvMesh::selectRefineCandidates
     {
         if (cellError[celli] > 0)
         {
-            candidateCell.set(celli, 1);
+            candidateCell.set(celli);
         }
     }
 }
@@ -798,7 +775,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
     calculateProtectedCells(unrefineableCell);
 
     // Count current selection
-    label nLocalCandidates = count(candidateCell, 1);
+    label nLocalCandidates = candidateCell.count();
     label nCandidates = returnReduce(nLocalCandidates, sumOp<label>());
 
     // Collect all cells
@@ -806,16 +783,12 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
 
     if (nCandidates < nTotToRefine)
     {
-        forAll(candidateCell, celli)
+        for (const label celli : candidateCell)
         {
             if
             (
-                cellLevel[celli] < maxRefinement
-             && candidateCell.get(celli)
-             && (
-                    unrefineableCell.empty()
-                 || !unrefineableCell.get(celli)
-                )
+                (!unrefineableCell.test(celli))
+             && cellLevel[celli] < maxRefinement
             )
             {
                 candidates.append(celli);
@@ -825,18 +798,14 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
     else
     {
         // Sort by error? For now just truncate.
-        for (label level = 0; level < maxRefinement; level++)
+        for (label level = 0; level < maxRefinement; ++level)
         {
-            forAll(candidateCell, celli)
+            for (const label celli : candidateCell)
             {
                 if
                 (
-                    cellLevel[celli] == level
-                 && candidateCell.get(celli)
-                 && (
-                        unrefineableCell.empty()
-                     || !unrefineableCell.get(celli)
-                    )
+                    (!unrefineableCell.test(celli))
+                 && cellLevel[celli] == level
                 )
                 {
                     candidates.append(celli);
@@ -889,17 +858,13 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
     if (protectedCell_.size())
     {
         // Get all points on a protected cell
-        forAll(pointCells, pointI)
+        forAll(pointCells, pointi)
         {
-            const labelList& pCells = pointCells[pointI];
-
-            forAll(pCells, pCellI)
+            for (const label celli : pointCells[pointi])
             {
-                label cellI = pCells[pCellI];
-
-                if (protectedCell_[cellI])
+                if (protectedCell_.test(celli))
                 {
-                    protectedPoint[pointI] = true;
+                    protectedPoint.set(pointi);
                     break;
                 }
             }
@@ -910,36 +875,29 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
             *this,
             protectedPoint,
             orEqOp<unsigned int>(),
-            0U
+            0u
         );
 
-        if (debug)
-        {
-            Info<< "From "
-                << returnReduce(protectedCell_.count(), sumOp<label>())
-                << " protected cells found "
-                << returnReduce(protectedPoint.count(), sumOp<label>())
-                << " protected points." << endl;
-        }
+        DebugInfo<< "From "
+            << returnReduce(protectedCell_.count(), sumOp<label>())
+            << " protected cells found "
+            << returnReduce(protectedPoint.count(), sumOp<label>())
+            << " protected points." << endl;
     }
 
 
     DynamicList<label> newSplitPoints(splitPoints.size());
 
-    forAll(splitPoints, i)
+    for (const label pointi : splitPoints)
     {
-        label pointi = splitPoints[i];
-
         if (!protectedPoint[pointi] && pFld[pointi] < unrefineLevel)
         {
             // Check that all cells are not marked
-            const labelList& pCells = pointCells[pointi];
-
             bool hasMarked = false;
 
-            forAll(pCells, pCelli)
+            for (const label celli : pointCells[pointi])
             {
-                if (markedCell.get(pCells[pCelli]))
+                if (markedCell.test(celli))
                 {
                     hasMarked = true;
                     break;
@@ -980,37 +938,29 @@ void Foam::dynamicRefineFvMesh::extendMarkedCells
 ) const
 {
     // Mark faces using any marked cell
-    boolList markedFace(nFaces(), false);
+    bitSet markedFace(nFaces());
 
-    forAll(markedCell, celli)
+    for (const label celli : markedCell)
     {
-        if (markedCell.get(celli))
-        {
-            const cell& cFaces = cells()[celli];
-
-            forAll(cFaces, i)
-            {
-                markedFace[cFaces[i]] = true;
-            }
-        }
+        markedFace.set(cells()[celli]);  // set multiple faces
     }
 
-    syncTools::syncFaceList(*this, markedFace, orEqOp<bool>());
+    syncTools::syncFaceList(*this, markedFace, orEqOp<unsigned int>());
 
     // Update cells using any markedFace
-    for (label facei = 0; facei < nInternalFaces(); facei++)
+    for (label facei = 0; facei < nInternalFaces(); ++facei)
     {
-        if (markedFace[facei])
+        if (markedFace.test(facei))
         {
-            markedCell.set(faceOwner()[facei], 1);
-            markedCell.set(faceNeighbour()[facei], 1);
+            markedCell.set(faceOwner()[facei]);
+            markedCell.set(faceNeighbour()[facei]);
         }
     }
-    for (label facei = nInternalFaces(); facei < nFaces(); facei++)
+    for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
     {
-        if (markedFace[facei])
+        if (markedFace.test(facei))
         {
-            markedCell.set(faceOwner()[facei], 1);
+            markedCell.set(faceOwner()[facei]);
         }
     }
 }
@@ -1018,37 +968,31 @@ void Foam::dynamicRefineFvMesh::extendMarkedCells
 
 void Foam::dynamicRefineFvMesh::checkEightAnchorPoints
 (
-    bitSet& protectedCell,
-    label& nProtected
+    bitSet& protectedCell
 ) const
 {
     const labelList& cellLevel = meshCutter_.cellLevel();
     const labelList& pointLevel = meshCutter_.pointLevel();
 
-    labelList nAnchorPoints(nCells(), 0);
+    labelList nAnchorPoints(nCells(), Zero);
 
     forAll(pointLevel, pointi)
     {
         const labelList& pCells = pointCells(pointi);
 
-        forAll(pCells, pCelli)
+        for (const label celli : pCells)
         {
-            label celli = pCells[pCelli];
-
             if (pointLevel[pointi] <= cellLevel[celli])
             {
                 // Check if cell has already 8 anchor points -> protect cell
                 if (nAnchorPoints[celli] == 8)
                 {
-                    if (protectedCell.set(celli, true))
-                    {
-                        nProtected++;
-                    }
+                    protectedCell.set(celli);
                 }
 
-                if (!protectedCell[celli])
+                if (!protectedCell.test(celli))
                 {
-                    nAnchorPoints[celli]++;
+                    ++nAnchorPoints[celli];
                 }
             }
         }
@@ -1057,10 +1001,9 @@ void Foam::dynamicRefineFvMesh::checkEightAnchorPoints
 
     forAll(protectedCell, celli)
     {
-        if (!protectedCell[celli] && nAnchorPoints[celli] != 8)
+        if (nAnchorPoints[celli] != 8)
         {
-            protectedCell.set(celli, true);
-            nProtected++;
+            protectedCell.set(celli);
         }
     }
 }
@@ -1072,9 +1015,9 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
 :
     dynamicFvMesh(io),
     meshCutter_(*this),
-    dumpLevel_(false),
+    protectedCell_(nCells()),
     nRefinementIterations_(0),
-    protectedCell_(nCells(), 0)
+    dumpLevel_(false)
 {
     // Read static part of dictionary
     readDict();
@@ -1091,28 +1034,23 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
     // Count number of points <= cellLevel
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    labelList nAnchors(nCells(), 0);
-
-    label nProtected = 0;
+    labelList nAnchors(nCells(), Zero);
 
     forAll(pointCells(), pointi)
     {
         const labelList& pCells = pointCells()[pointi];
 
-        forAll(pCells, i)
+        for (const label celli : pCells)
         {
-            label celli = pCells[i];
-
-            if (!protectedCell_.get(celli))
+            if (!protectedCell_.test(celli))
             {
                 if (pointLevel[pointi] <= cellLevel[celli])
                 {
-                    nAnchors[celli]++;
+                    ++nAnchors[celli];
 
                     if (nAnchors[celli] > 8)
                     {
-                        protectedCell_.set(celli, 1);
-                        nProtected++;
+                        protectedCell_.set(celli);
                     }
                 }
             }
@@ -1128,22 +1066,22 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
     {
         labelList neiLevel(nFaces());
 
-        for (label facei = 0; facei < nInternalFaces(); facei++)
+        for (label facei = 0; facei < nInternalFaces(); ++facei)
         {
             neiLevel[facei] = cellLevel[faceNeighbour()[facei]];
         }
-        for (label facei = nInternalFaces(); facei < nFaces(); facei++)
+        for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
         {
             neiLevel[facei] = cellLevel[faceOwner()[facei]];
         }
         syncTools::swapFaceList(*this, neiLevel);
 
 
-        boolList protectedFace(nFaces(), false);
+        bitSet protectedFace(nFaces());
 
         forAll(faceOwner(), facei)
         {
-            label faceLevel = max
+            const label faceLevel = max
             (
                 cellLevel[faceOwner()[facei]],
                 neiLevel[facei]
@@ -1153,39 +1091,36 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
 
             label nAnchors = 0;
 
-            forAll(f, fp)
+            for (const label pointi : f)
             {
-                if (pointLevel[f[fp]] <= faceLevel)
+                if (pointLevel[pointi] <= faceLevel)
                 {
-                    nAnchors++;
+                    ++nAnchors;
 
                     if (nAnchors > 4)
                     {
-                        protectedFace[facei] = true;
+                        protectedFace.set(facei);
                         break;
                     }
                 }
             }
         }
 
-        syncTools::syncFaceList(*this, protectedFace, orEqOp<bool>());
+        syncTools::syncFaceList(*this, protectedFace, orEqOp<unsigned int>());
 
-        for (label facei = 0; facei < nInternalFaces(); facei++)
+        for (label facei = 0; facei < nInternalFaces(); ++facei)
         {
-            if (protectedFace[facei])
+            if (protectedFace.test(facei))
             {
-                protectedCell_.set(faceOwner()[facei], 1);
-                nProtected++;
-                protectedCell_.set(faceNeighbour()[facei], 1);
-                nProtected++;
+                protectedCell_.set(faceOwner()[facei]);
+                protectedCell_.set(faceNeighbour()[facei]);
             }
         }
-        for (label facei = nInternalFaces(); facei < nFaces(); facei++)
+        for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
         {
-            if (protectedFace[facei])
+            if (protectedFace.test(facei))
             {
-                protectedCell_.set(faceOwner()[facei], 1);
-                nProtected++;
+                protectedCell_.set(faceOwner()[facei]);
             }
         }
 
@@ -1196,21 +1131,15 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
 
             if (cFaces.size() < 6)
             {
-                if (protectedCell_.set(celli, 1))
-                {
-                    nProtected++;
-                }
+                protectedCell_.set(celli);
             }
             else
             {
-                forAll(cFaces, cFacei)
+                for (const label cfacei : cFaces)
                 {
-                    if (faces()[cFaces[cFacei]].size() < 4)
+                    if (faces()[cfacei].size() < 4)
                     {
-                        if (protectedCell_.set(celli, 1))
-                        {
-                            nProtected++;
-                        }
+                        protectedCell_.set(celli);
                         break;
                     }
                 }
@@ -1218,26 +1147,24 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
         }
 
         // Check cells for 8 corner points
-        checkEightAnchorPoints(protectedCell_, nProtected);
+        checkEightAnchorPoints(protectedCell_);
     }
 
-    if (returnReduce(nProtected, sumOp<label>()) == 0)
+    if (!returnReduce(protectedCell_.any(), orOp<bool>()))
     {
         protectedCell_.clear();
     }
     else
     {
+        cellSet protectedCells
+        (
+            *this,
+            "protectedCells",
+            HashSetOps::used(protectedCell_)
+        );
 
-        cellSet protectedCells(*this, "protectedCells", nProtected);
-        forAll(protectedCell_, celli)
-        {
-            if (protectedCell_[celli])
-            {
-                protectedCells.insert(celli);
-            }
-        }
-
-        Info<< "Detected " << returnReduce(nProtected, sumOp<label>())
+        Info<< "Detected "
+            << returnReduce(protectedCells.size(), sumOp<label>())
             << " cells that are protected from refinement."
             << " Writing these to cellSet "
             << protectedCells.name()
@@ -1248,12 +1175,6 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::dynamicRefineFvMesh::~dynamicRefineFvMesh()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::dynamicRefineFvMesh::update()
@@ -1277,7 +1198,7 @@ bool Foam::dynamicRefineFvMesh::update()
         ).optionalSubDict(typeName + "Coeffs")
     );
 
-    label refineInterval = refineDict.get<label>("refineInterval");
+    const label refineInterval = refineDict.get<label>("refineInterval");
 
     bool hasChanged = false;
 
@@ -1297,14 +1218,12 @@ bool Foam::dynamicRefineFvMesh::update()
     }
 
 
-
-
     // Note: cannot refine at time 0 since no V0 present since mesh not
     //       moved yet.
 
     if (time().timeIndex() > 0 && time().timeIndex() % refineInterval == 0)
     {
-        label maxCells = refineDict.get<label>("maxCells");
+        const label maxCells = refineDict.get<label>("maxCells");
 
         if (maxCells <= 0)
         {
@@ -1315,7 +1234,7 @@ bool Foam::dynamicRefineFvMesh::update()
                 << exit(FatalError);
         }
 
-        label maxRefinement = refineDict.get<label>("maxRefinement");
+        const label maxRefinement = refineDict.get<label>("maxRefinement");
 
         if (maxRefinement <= 0)
         {
@@ -1339,8 +1258,7 @@ bool Foam::dynamicRefineFvMesh::update()
             "unrefineLevel",
             GREAT
         );
-        const label nBufferLayers =
-            refineDict.get<label>("nBufferLayers");
+        const label nBufferLayers = refineDict.get<label>("nBufferLayers");
 
         // Cells marked for refinement or otherwise protected from unrefinement.
         bitSet refineCell(nCells());
@@ -1368,7 +1286,7 @@ bool Foam::dynamicRefineFvMesh::update()
                 )
             );
 
-            label nCellsToRefine = returnReduce
+            const label nCellsToRefine = returnReduce
             (
                 cellsToRefine.size(), sumOp<label>()
             );
@@ -1388,19 +1306,16 @@ bool Foam::dynamicRefineFvMesh::update()
 
                     forAll(cellMap, celli)
                     {
-                        label oldCelli = cellMap[celli];
+                        const label oldCelli = cellMap[celli];
 
-                        if (oldCelli < 0)
+                        if
+                        (
+                            (oldCelli < 0)
+                         || (reverseCellMap[oldCelli] != celli)
+                         || (refineCell.test(oldCelli))
+                        )
                         {
-                            newRefineCell.set(celli, 1);
-                        }
-                        else if (reverseCellMap[oldCelli] != celli)
-                        {
-                            newRefineCell.set(celli, 1);
-                        }
-                        else
-                        {
-                            newRefineCell.set(celli, refineCell.get(oldCelli));
+                            newRefineCell.set(celli);
                         }
                     }
                     refineCell.transfer(newRefineCell);
@@ -1408,7 +1323,7 @@ bool Foam::dynamicRefineFvMesh::update()
 
                 // Extend with a buffer layer to prevent neighbouring points
                 // being unrefined.
-                for (label i = 0; i < nBufferLayers; i++)
+                for (label i = 0; i < nBufferLayers; ++i)
                 {
                     extendMarkedCells(refineCell);
                 }
@@ -1430,7 +1345,7 @@ bool Foam::dynamicRefineFvMesh::update()
                 )
             );
 
-            label nSplitPoints = returnReduce
+            const label nSplitPoints = returnReduce
             (
                 pointsToUnrefine.size(),
                 sumOp<label>()
@@ -1498,7 +1413,7 @@ bool Foam::dynamicRefineFvMesh::writeObject
                 false
             ),
             *this,
-            dimensionedScalar("level", dimless, 0)
+            dimensionedScalar(dimless, Zero)
         );
 
         const labelList& cellLevel = meshCutter_.cellLevel();
diff --git a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.H b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.H
index 2ccd324486..84d9fa7fe9 100644
--- a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.H
+++ b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.H
@@ -29,7 +29,8 @@ Description
 
     Determines which cells to refine/unrefine and does all in update().
 
-
+    \verbatim
+    {
         // How often to refine
         refineInterval  1;
         // Field to be refinement on
@@ -59,6 +60,8 @@ Description
 
         // Write the refinement level as a volScalarField
         dumpLevel       true;
+    }
+    \endverbatim
 
 
 SourceFiles
@@ -91,24 +94,21 @@ protected:
         //- Mesh cutting engine
         hexRef8 meshCutter_;
 
-        //- Dump cellLevel for post-processing
-        bool dumpLevel_;
-
         //- Fluxes to map
         HashTable<word> correctFluxes_;
 
-        //- Number of refinement/unrefinement steps done so far.
-        label nRefinementIterations_;
-
         //- Protected cells (usually since not hexes)
         bitSet protectedCell_;
 
+        //- Number of refinement/unrefinement steps done so far.
+        label nRefinementIterations_;
+
+        //- Dump cellLevel for post-processing
+        bool dumpLevel_;
+
 
     // Protected Member Functions
 
-        //- Count set/unset elements in packedlist.
-        static label count(const bitSet&, const unsigned int);
-
         //- Calculate cells that cannot be refined since would trigger
         //  refinement of protectedCell_ (since 2:1 refinement cascade)
         void calculateProtectedCells(bitSet& unrefineableCell) const;
@@ -180,11 +180,7 @@ protected:
             void extendMarkedCells(bitSet& markedCell) const;
 
             //- Check all cells have 8 anchor points
-            void checkEightAnchorPoints
-            (
-                bitSet& protectedCell,
-                label& nProtected
-            ) const;
+            void checkEightAnchorPoints(bitSet& protectedCell) const;
 
             //- Map single non-flux surface<Type>Field
             //  for new internal faces (e.g. AMR refine). This currently
@@ -233,7 +229,7 @@ public:
 
 
     //- Destructor
-    virtual ~dynamicRefineFvMesh();
+    virtual ~dynamicRefineFvMesh() = default;
 
 
     // Member Functions