diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/removeCells.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/removeCells.C
index dc8d7aaa6f..38269d63a6 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/removeCells.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/removeCells.C
@@ -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.
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "removeCells.H"
+#include "bitSet.H"
 #include "polyMesh.H"
 #include "polyTopoChange.H"
 #include "polyRemoveCell.H"
@@ -39,23 +40,43 @@ namespace Foam
     defineTypeNameAndDebug(removeCells, 0);
 }
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
-void Foam::removeCells::uncount
-(
-    const labelList& f,
-    labelList& nUsage
-)
+namespace
 {
-    forAll(f, fp)
+
+// Increase count (usage) of elements of list
+inline void incrCount(const Foam::labelUList& list, Foam::labelList& counter)
+{
+    for (auto idx : list)
     {
-        nUsage[f[fp]]--;
+        ++counter[idx];
     }
 }
 
+// Decrease count (usage) of elements of list
+inline void decrCount(const Foam::labelUList& list, Foam::labelList& counter)
+{
+    for (auto idx : list)
+    {
+        --counter[idx];
+    }
+}
+
+} // End anonymous namespace
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::removeCells::removeCells
+(
+    const polyMesh& mesh
+)
+:
+    removeCells(mesh, true)
+{}
+
+
 Foam::removeCells::removeCells
 (
     const polyMesh& mesh,
@@ -71,47 +92,37 @@ Foam::removeCells::removeCells
 
 Foam::labelList Foam::removeCells::getExposedFaces
 (
-    const labelList& cellLabels
+    const bitSet& removedCell
 ) const
 {
-    // Create list of cells to be removed
-    boolList removedCell(mesh_.nCells(), false);
-
-    // Go from labelList of cells-to-remove to a boolList.
-    forAll(cellLabels, i)
-    {
-        removedCell[cellLabels[i]] = true;
-    }
-
-
     const labelList& faceOwner = mesh_.faceOwner();
     const labelList& faceNeighbour = mesh_.faceNeighbour();
 
     // Count cells using face.
     labelList nCellsUsingFace(mesh_.nFaces(), 0);
 
-    for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
+    for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
     {
-        label own = faceOwner[facei];
-        label nei = faceNeighbour[facei];
+        const label own = faceOwner[facei];
+        const label nei = faceNeighbour[facei];
 
         if (!removedCell[own])
         {
-            nCellsUsingFace[facei]++;
+            ++nCellsUsingFace[facei];
         }
         if (!removedCell[nei])
         {
-            nCellsUsingFace[facei]++;
+            ++nCellsUsingFace[facei];
         }
     }
 
-    for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
+    for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); ++facei)
     {
-        label own = faceOwner[facei];
+        const label own = faceOwner[facei];
 
         if (!removedCell[own])
         {
-            nCellsUsingFace[facei]++;
+            ++nCellsUsingFace[facei];
         }
     }
 
@@ -147,7 +158,7 @@ Foam::labelList Foam::removeCells::getExposedFaces
 
     DynamicList<label> exposedFaces(mesh_.nFaces()/10);
 
-    for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
+    for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
     {
         if (nCellsUsingFace[facei] == 1)
         {
@@ -157,17 +168,15 @@ Foam::labelList Foam::removeCells::getExposedFaces
 
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    forAll(patches, patchi)
+    for (const polyPatch& pp : patches)
     {
-        const polyPatch& pp = patches[patchi];
-
         if (pp.coupled())
         {
             label facei = pp.start();
 
             forAll(pp, i)
             {
-                label own = faceOwner[facei];
+                const label own = faceOwner[facei];
 
                 if (nCellsUsingFace[facei] == 1 && !removedCell[own])
                 {
@@ -176,7 +185,7 @@ Foam::labelList Foam::removeCells::getExposedFaces
                     exposedFaces.append(facei);
                 }
 
-                facei++;
+                ++facei;
             }
         }
     }
@@ -187,9 +196,9 @@ Foam::labelList Foam::removeCells::getExposedFaces
 
 void Foam::removeCells::setRefinement
 (
-    const labelList& cellLabels,
-    const labelList& exposedFaceLabels,
-    const labelList& exposedPatchIDs,
+    const bitSet& removedCell,
+    const labelUList& exposedFaceLabels,
+    const labelUList& exposedPatchIDs,
     polyTopoChange& meshMod
 ) const
 {
@@ -209,13 +218,14 @@ void Foam::removeCells::setRefinement
 
     forAll(exposedFaceLabels, i)
     {
-        label patchi = exposedPatchIDs[i];
+        const label facei = exposedFaceLabels[i];
+        const label patchi = exposedPatchIDs[i];
 
         if (patchi < 0 || patchi >= patches.size())
         {
             FatalErrorInFunction
                 << "Invalid patch " << patchi
-                << " for exposed face " << exposedFaceLabels[i] << endl
+                << " for exposed face " << facei << nl
                 << "Valid patches 0.." << patches.size()-1
                 << abort(FatalError);
         }
@@ -223,28 +233,20 @@ void Foam::removeCells::setRefinement
         if (patches[patchi].coupled())
         {
             FatalErrorInFunction
-                << "Trying to put exposed face " << exposedFaceLabels[i]
+                << "Trying to put exposed face " << facei
                 << " into a coupled patch : " << patches[patchi].name()
-                << endl
+                << nl
                 << "This is illegal."
                 << abort(FatalError);
         }
 
-        newPatchID[exposedFaceLabels[i]] = patchi;
+        newPatchID[facei] = patchi;
     }
 
 
-    // Create list of cells to be removed
-    boolList removedCell(mesh_.nCells(), false);
-
-    // Go from labelList of cells-to-remove to a boolList and remove all
-    // cells mentioned.
-    forAll(cellLabels, i)
+    // Walk all the cells mentioned for removal
+    for (const label celli : removedCell)
     {
-        label celli = cellLabels[i];
-
-        removedCell[celli] = true;
-
         //Pout<< "Removing cell " << celli
         //    << " cc:" << mesh_.cellCentres()[celli] << endl;
 
@@ -260,26 +262,21 @@ void Foam::removeCells::setRefinement
     const labelList& faceNeighbour = mesh_.faceNeighbour();
     const faceZoneMesh& faceZones = mesh_.faceZones();
 
-    // Count starting number of faces using each point. Keep up to date whenever
-    // removing a face.
+    // Count starting number of faces using each point.
+    // Update whenever removing a face.
     labelList nFacesUsingPoint(mesh_.nPoints(), 0);
 
-    forAll(faces, facei)
+    for (const face& f : faces)
     {
-        const face& f = faces[facei];
-
-        forAll(f, fp)
-        {
-            nFacesUsingPoint[f[fp]]++;
-        }
+        incrCount(f, nFacesUsingPoint);
     }
 
 
-    for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
+    for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
     {
         const face& f = faces[facei];
-        label own = faceOwner[facei];
-        label nei = faceNeighbour[facei];
+        const label own = faceOwner[facei];
+        const label nei = faceNeighbour[facei];
 
         if (removedCell[own])
         {
@@ -290,7 +287,7 @@ void Foam::removeCells::setRefinement
                 //    << " fc:" << mesh_.faceCentres()[facei] << endl;
 
                 meshMod.setAction(polyRemoveFace(facei));
-                uncount(f, nFacesUsingPoint);
+                decrCount(f, nFacesUsingPoint);
             }
             else
             {
@@ -305,7 +302,7 @@ void Foam::removeCells::setRefinement
 
                 // nei is remaining cell. Facei becomes external cell
 
-                label zoneID = faceZones.whichZone(facei);
+                const label zoneID = faceZones.whichZone(facei);
                 bool zoneFlip = false;
 
                 if (zoneID >= 0)
@@ -353,7 +350,7 @@ void Foam::removeCells::setRefinement
             //    << " into patch " << newPatchID[facei] << endl;
 
             // own is remaining cell. Facei becomes external cell.
-            label zoneID = faceZones.whichZone(facei);
+            const label zoneID = faceZones.whichZone(facei);
             bool zoneFlip = false;
 
             if (zoneID >= 0)
@@ -380,10 +377,8 @@ void Foam::removeCells::setRefinement
         }
     }
 
-    forAll(patches, patchi)
+    for (const polyPatch& pp : patches)
     {
-        const polyPatch& pp = patches[patchi];
-
         if (pp.coupled())
         {
             label facei = pp.start();
@@ -396,7 +391,7 @@ void Foam::removeCells::setRefinement
                     //    << " fc:" << mesh_.faceCentres()[facei]
                     //    << " into patch " << newPatchID[facei] << endl;
 
-                    label zoneID = faceZones.whichZone(facei);
+                    const label zoneID = faceZones.whichZone(facei);
                     bool zoneFlip = false;
 
                     if (zoneID >= 0)
@@ -429,10 +424,10 @@ void Foam::removeCells::setRefinement
                     //    << endl;
 
                     meshMod.setAction(polyRemoveFace(facei));
-                    uncount(faces[facei], nFacesUsingPoint);
+                    decrCount(faces[facei], nFacesUsingPoint);
                 }
 
-                facei++;
+                ++facei;
             }
         }
         else
@@ -457,10 +452,10 @@ void Foam::removeCells::setRefinement
                     //    << endl;
 
                     meshMod.setAction(polyRemoveFace(facei));
-                    uncount(faces[facei], nFacesUsingPoint);
+                    decrCount(faces[facei], nFacesUsingPoint);
                 }
 
-                facei++;
+                ++facei;
             }
         }
     }
@@ -491,4 +486,35 @@ void Foam::removeCells::setRefinement
 }
 
 
+Foam::labelList Foam::removeCells::getExposedFaces
+(
+    const labelUList& cellsToRemove
+) const
+{
+    bitSet removeCell(mesh_.nCells(), cellsToRemove);
+
+    return getExposedFaces(removeCell);
+}
+
+
+void Foam::removeCells::setRefinement
+(
+    const labelUList& cellsToRemove,
+    const labelUList& exposedFaceLabels,
+    const labelUList& exposedPatchIDs,
+    polyTopoChange& meshMod
+) const
+{
+    bitSet removedCell(mesh_.nCells(), cellsToRemove);
+
+    setRefinement
+    (
+        removedCell,
+        exposedFaceLabels,
+        exposedPatchIDs,
+        meshMod
+    );
+}
+
+
 // ************************************************************************* //
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/removeCells.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/removeCells.H
index b27a9062ba..5088ffe03a 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/removeCells.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/removeCells.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,7 +25,7 @@ Class
     Foam::removeCells
 
 Description
-    Given list of cells to remove insert all the topology changes.
+    Given list of cells to remove, insert all the topology changes.
 
     Works in two passes:
     - get faces that will become boundary faces
@@ -48,7 +48,8 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
+// Forward declarations
+class bitSet;
 class polyMesh;
 class polyTopoChange;
 class mapPolyMesh;
@@ -67,17 +68,6 @@ class removeCells
         //- Whether or not to synchronize parallel case.
         const bool syncPar_;
 
-
-    // Private Member Functions
-
-        //- Decrease count of elements of f
-        static void uncount
-        (
-            const labelList& f,
-            labelList& nUsage
-        );
-
-
 public:
 
     //- Runtime type information
@@ -86,35 +76,58 @@ public:
 
     // Constructors
 
-        //- Construct from mesh. syncPar: do parallel synchronization.
-        removeCells(const polyMesh& mesh, const bool syncPar = true);
+        //- Construct from mesh. With parallel synchronization.
+        explicit removeCells(const polyMesh& mesh);
+
+        //- Construct from mesh, optionally with parallel synchronization.
+        removeCells(const polyMesh& mesh, const bool syncPar);
+
+
+    //- Destructor
+    ~removeCells() = default;
 
 
     // Member Functions
 
-        // Topology changes
+        //- Get labels of faces exposed after cells removal.
+        //  These are
+        //  - internal faces that become boundary faces
+        //  - coupled faces that become uncoupled (since one of the sides
+        //    gets deleted)
+        labelList getExposedFaces(const bitSet& removedCell) const;
 
-            //- Get labels of exposed faces.
-            //  These are
-            //  - internal faces that become boundary faces
-            //  - coupled faces that become uncoupled (since one of the sides
-            //    gets deleted)
-            labelList getExposedFaces(const labelList& cellsToRemove) const;
+        //- Get labels of faces exposed after cells removal.
+        //  These are
+        //  - internal faces that become boundary faces
+        //  - coupled faces that become uncoupled (since one of the sides
+        //    gets deleted)
+        labelList getExposedFaces(const labelUList& cellsToRemove) const;
 
-            //- Play commands into polyTopoChange to remove cells.
-            //  patchIDs is for every element in facesToExpose (see above) the
-            //  patch it has to go into. This cannot be a coupled patch!
-            void setRefinement
-            (
-                const labelList& cellsToRemove,
-                const labelList& facesToExpose,
-                const labelList& patchIDs,
-                polyTopoChange&
-            ) const;
+        //- Play commands into polyTopoChange to remove cells.
+        //  patchIDs is for every element in facesToExpose (see above) the
+        //  patch it has to go into. This cannot be a coupled patch!
+        void setRefinement
+        (
+            const bitSet& removedCell,
+            const labelUList& facesToExpose,
+            const labelUList& patchIDs,
+            polyTopoChange&
+        ) const;
 
-            //- Force recalculation of locally stored data on topological change
-            void updateMesh(const mapPolyMesh&)
-            {}
+        //- Play commands into polyTopoChange to remove cells.
+        //  patchIDs is for every element in facesToExpose (see above) the
+        //  patch it has to go into. This cannot be a coupled patch!
+        void setRefinement
+        (
+            const labelUList& cellsToRemove,
+            const labelUList& facesToExpose,
+            const labelUList& patchIDs,
+            polyTopoChange&
+        ) const;
+
+        //- Force recalculation of locally stored data on topological change
+        void updateMesh(const mapPolyMesh&)
+        {}
 };