diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C
index 8c21017269..6bcc2af56f 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C
@@ -66,7 +66,6 @@ primitiveMesh::primitiveMesh()
     ppPtr_(NULL),
     cpPtr_(NULL),
 
-    allocSize_(0),
     labels_(0),
 
     cellCentresPtr_(NULL),
@@ -109,7 +108,6 @@ primitiveMesh::primitiveMesh
     ppPtr_(NULL),
     cpPtr_(NULL),
 
-    allocSize_(0),
     labels_(0),
 
     cellCentresPtr_(NULL),
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
index 3868ba2b25..6f24944cbc 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
@@ -54,6 +54,7 @@ SourceFiles
 #ifndef primitiveMesh_H
 #define primitiveMesh_H
 
+#include "DynamicList.H"
 #include "edgeList.H"
 #include "pointField.H"
 #include "SubField.H"
@@ -157,10 +158,8 @@ class primitiveMesh
 
         // On-the-fly edge addresing storage
 
-            //- Temporary storage for addressing. allocSize is the real size
-            //  of the labelList.
-            mutable label allocSize_;
-            mutable labelList labels_;
+            //- Temporary storage for addressing.
+            mutable DynamicList<label> labels_;
 
             //- Temporary storage for addressing
             mutable labelHashSet labelSet_;
@@ -705,31 +704,80 @@ public:
 
             // On-the-fly addressing calculation. These functions return either
             // a reference to the full addressing (if already calculated) or
-            // a reference to member data labels_ so be careful when not storing
+            // a reference to the supplied storage. The one-argument ones
+            // use member DynamicList labels_ so be careful when not storing
             // result.
 
-            //- cellCells using cells
+            //- cellCells using cells.
+            const labelList& cellCells
+            (
+                const label cellI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& cellCells(const label cellI) const;
 
             //- cellPoints using cells
+            const labelList& cellPoints
+            (
+                const label cellI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& cellPoints(const label cellI) const;
 
             //- pointCells using pointFaces
+            const labelList& pointCells
+            (
+                const label pointI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& pointCells(const label pointI) const;
 
             //- pointPoints using edges, pointEdges
+            const labelList& pointPoints
+            (
+                const label pointI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& pointPoints(const label pointI) const;
 
             //- faceEdges using pointFaces, edges, pointEdges
+            const labelList& faceEdges
+            (
+                const label faceI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& faceEdges(const label faceI) const;
 
             //- edgeFaces using pointFaces, edges, pointEdges
+            const labelList& edgeFaces
+            (
+                const label edgeI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& edgeFaces(const label edgeI) const;
 
             //- edgeCells using pointFaces, edges, pointEdges
+            const labelList& edgeCells
+            (
+                const label edgeI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& edgeCells(const label edgeI) const;
 
             //- cellEdges using cells, pointFaces, edges, pointEdges
+            const labelList& cellEdges
+            (
+                const label cellI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& cellEdges(const label cellI) const;
 
 
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCells.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCells.C
index 3c8f896bbe..e2eaf17a0f 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCells.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCells.C
@@ -41,6 +41,14 @@ void primitiveMesh::calcCellCells() const
     {
         Pout<< "primitiveMesh::calcCellCells() : calculating cellCells"
             << endl;
+
+        if (debug == -1)
+        {
+            // For checking calls:abort so we can quickly hunt down
+            // origin of call
+            FatalErrorIn("primitiveMesh::calcCellCells()")
+                << abort(FatalError);
+        }
     }
 
     // It is an error to attempt to recalculate cellCells
@@ -105,7 +113,11 @@ const labelListList& primitiveMesh::cellCells() const
 }
 
 
-const labelList& primitiveMesh::cellCells(const label cellI) const
+const labelList& primitiveMesh::cellCells
+(
+    const label cellI,
+    DynamicList<label>& storage
+) const
 {
     if (hasCellCells())
     {
@@ -117,16 +129,7 @@ const labelList& primitiveMesh::cellCells(const label cellI) const
         const labelList& nei = faceNeighbour();
         const cell& cFaces = cells()[cellI];
 
-        labels_.size() = allocSize_;
-
-        if (cFaces.size() > allocSize_)
-        {
-            labels_.clear();
-            allocSize_ = cFaces.size();
-            labels_.setSize(allocSize_);
-        }
-
-        label n = 0;
+        storage.clear();
 
         forAll(cFaces, i)
         {
@@ -136,22 +139,26 @@ const labelList& primitiveMesh::cellCells(const label cellI) const
             {
                 if (own[faceI] == cellI)
                 {
-                    labels_[n++] = nei[faceI];
+                    storage.append(nei[faceI]);
                 }
                 else
                 {
-                    labels_[n++] = own[faceI];
+                    storage.append(own[faceI]);
                 }
             }
         }
 
-        labels_.size() = n;
-
-        return labels_;
+        return storage;
     }
 }
 
 
+const labelList& primitiveMesh::cellCells(const label cellI) const
+{
+    return cellCells(cellI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellEdges.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellEdges.C
index 6588f4b721..47d6343f16 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellEdges.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellEdges.C
@@ -40,6 +40,14 @@ void Foam::primitiveMesh::calcCellEdges() const
         Pout<< "primitiveMesh::calcCellEdges() : "
             << "calculating cellEdges"
             << endl;
+
+        if (debug == -1)
+        {
+            // For checking calls:abort so we can quickly hunt down
+            // origin of call
+            FatalErrorIn("primitiveMesh::calcCellEdges()")
+                << abort(FatalError);
+        }
     }
 
     // It is an error to attempt to recalculate cellEdges
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellPoints.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellPoints.C
index c8f63e1eef..f813c3b15f 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellPoints.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellPoints.C
@@ -42,6 +42,14 @@ const labelListList& primitiveMesh::cellPoints() const
         {
             Pout<< "primitiveMesh::cellPoints() : "
                 << "calculating cellPoints" << endl;
+
+            if (debug == -1)
+            {
+                // For checking calls:abort so we can quickly hunt down
+                // origin of call
+                FatalErrorIn("primitiveMesh::cellPoints()")
+                    << abort(FatalError);
+            }
         }
 
         // Invert pointCells
@@ -53,7 +61,11 @@ const labelListList& primitiveMesh::cellPoints() const
 }
 
 
-const labelList& primitiveMesh::cellPoints(const label cellI) const
+const labelList& primitiveMesh::cellPoints
+(
+    const label cellI,
+    DynamicList<label>& storage
+) const
 {
     if (hasCellPoints())
     {
@@ -76,29 +88,28 @@ const labelList& primitiveMesh::cellPoints(const label cellI) const
             }
         }
 
-        labels_.size() = allocSize_;
-
-        if (labelSet_.size() > allocSize_)
+        storage.clear();
+        if (labelSet_.size() > storage.allocSize())
         {
-            labels_.clear();
-            allocSize_ = labelSet_.size();
-            labels_.setSize(allocSize_);
+            storage.setSize(labelSet_.size());
         }
 
-        label n = 0;
-
         forAllConstIter(labelHashSet, labelSet_, iter)
         {
-            labels_[n++] = iter.key();
+            storage.append(iter.key());
         }
 
-        labels_.size() = n;
-
-        return labels_;
+        return storage;
     }
 }
 
 
+const labelList& primitiveMesh::cellPoints(const label cellI) const
+{
+    return cellPoints(cellI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeCells.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeCells.C
index 6a2f309dd3..7b56c2e88f 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeCells.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeCells.C
@@ -41,6 +41,14 @@ const labelListList& primitiveMesh::edgeCells() const
         if (debug)
         {
             Pout<< "primitiveMesh::edgeCells() : calculating edgeCells" << endl;
+
+            if (debug == -1)
+            {
+                // For checking calls:abort so we can quickly hunt down
+                // origin of call
+                FatalErrorIn("primitiveMesh::edgeCells()")
+                    << abort(FatalError);
+            }
         }
         // Invert cellEdges
         ecPtr_ = new labelListList(nEdges());
@@ -51,7 +59,11 @@ const labelListList& primitiveMesh::edgeCells() const
 }
 
 
-const labelList& primitiveMesh::edgeCells(const label edgeI) const
+const labelList& primitiveMesh::edgeCells
+(
+    const label edgeI,
+    DynamicList<label>& storage
+) const
 {
     if (hasEdgeCells())
     {
@@ -62,24 +74,11 @@ const labelList& primitiveMesh::edgeCells(const label edgeI) const
         const labelList& own = faceOwner();
         const labelList& nei = faceNeighbour();
 
-        // edge faces can either return labels_ or reference in edgeLabels.
-        labelList labelsCopy;
-        if (!hasEdgeFaces())
-        {
-            labelsCopy = edgeFaces(edgeI);
-        }
+        // Construct edgeFaces
+        DynamicList<label> eFacesStorage;
+        const labelList& eFaces = edgeFaces(edgeI, eFacesStorage);
 
-        const labelList& eFaces =
-        (
-            hasEdgeFaces()
-          ? edgeFaces()[edgeI]
-          : labelsCopy
-        );
-
-        labels_.size() = allocSize_;
-
-        // labels_ should certainly be big enough for edge cells.
-        label n = 0;
+        storage.clear();
 
         // Do quadratic insertion.
         forAll(eFaces, i)
@@ -89,10 +88,10 @@ const labelList& primitiveMesh::edgeCells(const label edgeI) const
             {
                 label ownCellI = own[faceI];
 
-                // Check if not already in labels_
-                for (label j = 0; j < n; j++)
+                // Check if not already in storage
+                forAll(storage, j)
                 {
-                    if (labels_[j] == ownCellI)
+                    if (storage[j] == ownCellI)
                     {
                         ownCellI = -1;
                         break;
@@ -101,7 +100,7 @@ const labelList& primitiveMesh::edgeCells(const label edgeI) const
 
                 if (ownCellI != -1)
                 {
-                    labels_[n++] = ownCellI;
+                    storage.append(ownCellI);
                 }
             }
 
@@ -109,9 +108,9 @@ const labelList& primitiveMesh::edgeCells(const label edgeI) const
             {
                 label neiCellI = nei[faceI];
 
-                for (label j = 0; j < n; j++)
+                forAll(storage, j)
                 {
-                    if (labels_[j] == neiCellI)
+                    if (storage[j] == neiCellI)
                     {
                         neiCellI = -1;
                         break;
@@ -120,18 +119,22 @@ const labelList& primitiveMesh::edgeCells(const label edgeI) const
 
                 if (neiCellI != -1)
                 {
-                    labels_[n++] = neiCellI;
+                    storage.append(neiCellI);
                 }
             }
         }
 
-        labels_.size() = n;
-
-        return labels_;
+        return storage;
     }
 }
 
 
+const labelList& primitiveMesh::edgeCells(const label edgeI) const
+{
+    return edgeCells(edgeI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeFaces.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeFaces.C
index 99536bcd33..1281a713af 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeFaces.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeFaces.C
@@ -41,6 +41,14 @@ const labelListList& primitiveMesh::edgeFaces() const
         if (debug)
         {
             Pout<< "primitiveMesh::edgeFaces() : calculating edgeFaces" << endl;
+
+            if (debug == -1)
+            {
+                // For checking calls:abort so we can quickly hunt down
+                // origin of call
+                FatalErrorIn("primitiveMesh::edgeFaces()")
+                    << abort(FatalError);
+            }
         }
 
         // Invert faceEdges
@@ -52,7 +60,11 @@ const labelListList& primitiveMesh::edgeFaces() const
 }
 
 
-const labelList& primitiveMesh::edgeFaces(const label edgeI) const
+const labelList& primitiveMesh::edgeFaces
+(
+    const label edgeI,
+    DynamicList<label>& storage
+) const
 {
     if (hasEdgeFaces())
     {
@@ -67,9 +79,8 @@ const labelList& primitiveMesh::edgeFaces(const label edgeI) const
 
         label i0 = 0;
         label i1 = 0;
-        label n = 0;
 
-        labels_.size() = allocSize_;
+        storage.clear();
 
         while (i0 < pFaces0.size() && i1 < pFaces1.size())
         {
@@ -84,26 +95,23 @@ const labelList& primitiveMesh::edgeFaces(const label edgeI) const
             else
             {
                 // Equal. Append.
-                if (n == allocSize_)
-                {
-                    // Have setSize copy contents so far
-                    labels_.size() = n;
-                    allocSize_ = allocSize_*2 + 1;
-                    labels_.setSize(allocSize_);
-                }
-                labels_[n++] = pFaces0[i0];
+                storage.append(pFaces0[i0]);
                 ++i0;
                 ++i1;
             }
         }
 
-        labels_.size() = n;
-
-        return labels_;
+        return storage;
     }
 }
 
 
+const labelList& primitiveMesh::edgeFaces(const label edgeI) const
+{
+    return edgeFaces(edgeI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C
index b64eda7266..d7a79fb8fe 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C
@@ -582,11 +582,15 @@ void primitiveMesh::clearOutEdges()
     deleteDemandDrivenData(pePtr_);
     deleteDemandDrivenData(fePtr_);
     labels_.clear();
-    allocSize_ = 0;
+    labelSet_.clear();
 }
 
 
-const labelList& primitiveMesh::faceEdges(const label faceI) const
+const labelList& primitiveMesh::faceEdges
+(
+    const label faceI,
+    DynamicList<label>& storage
+) const
 {
     if (hasFaceEdges())
     {
@@ -597,34 +601,40 @@ const labelList& primitiveMesh::faceEdges(const label faceI) const
         const labelListList& pointEs = pointEdges();
         const face& f = faces()[faceI];
 
-        labels_.size() = allocSize_;
-
-        if (f.size() > allocSize_)
+        storage.clear();
+        if (f.size() > storage.allocSize())
         {
-            labels_.clear();
-            allocSize_ = f.size();
-            labels_.setSize(allocSize_);
+            storage.setSize(f.size());
         }
 
-        label n = 0;
-
         forAll(f, fp)
         {
-            labels_[n++] = findFirstCommonElementFromSortedLists
+            storage.append
             (
-                pointEs[f[fp]],
-                pointEs[f.nextLabel(fp)]
+                findFirstCommonElementFromSortedLists
+                (
+                    pointEs[f[fp]],
+                    pointEs[f.nextLabel(fp)]
+                )
             );
         }
 
-        labels_.size() = n;
-
-        return labels_;
+        return storage;
     }
 }
 
 
-const labelList& primitiveMesh::cellEdges(const label cellI) const
+const labelList& primitiveMesh::faceEdges(const label faceI) const
+{
+    return faceEdges(faceI, labels_);
+}
+
+
+const labelList& primitiveMesh::cellEdges
+(
+    const label cellI,
+    DynamicList<label>& storage
+) const
 {
     if (hasCellEdges())
     {
@@ -646,29 +656,29 @@ const labelList& primitiveMesh::cellEdges(const label cellI) const
             }
         }
 
-        labels_.size() = allocSize_;
+        storage.clear();
 
-        if (labelSet_.size() > allocSize_)
+        if (labelSet_.size() > storage.allocSize())
         {
-            labels_.clear();
-            allocSize_ = labelSet_.size();
-            labels_.setSize(allocSize_);
+            storage.setSize(labelSet_.size());
         }
 
-        label n =0;
-
         forAllConstIter(labelHashSet, labelSet_, iter)
         {
-            labels_[n++] = iter.key();
+            storage.append(iter.key());
         }
 
-        labels_.size() = n;
-
-        return labels_;
+        return storage;
     }
 }
 
 
+const labelList& primitiveMesh::cellEdges(const label cellI) const
+{
+    return cellEdges(cellI, labels_);;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointCells.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointCells.C
index 6d7af72e0b..2986863d97 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointCells.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointCells.C
@@ -43,6 +43,14 @@ void primitiveMesh::calcPointCells() const
         Pout<< "primitiveMesh::calcPointCells() : "
             << "calculating pointCells"
             << endl;
+
+        if (debug == -1)
+        {
+            // For checking calls:abort so we can quickly hunt down
+            // origin of call
+            FatalErrorIn("primitiveMesh::calcPointCells()")
+                << abort(FatalError);
+        }
     }
 
     // It is an error to attempt to recalculate pointCells
@@ -114,7 +122,11 @@ const labelListList& primitiveMesh::pointCells() const
 }
 
 
-const labelList& primitiveMesh::pointCells(const label pointI) const
+const labelList& primitiveMesh::pointCells
+(
+    const label pointI,
+    DynamicList<label>& storage
+) const
 {
     if (hasPointCells())
     {
@@ -126,58 +138,48 @@ const labelList& primitiveMesh::pointCells(const label pointI) const
         const labelList& nei = faceNeighbour();
         const labelList& pFaces = pointFaces()[pointI];
 
-        labels_.size() = allocSize_;
-
-        label n = 0;
+        storage.clear();
 
         forAll(pFaces, i)
         {
             const label faceI = pFaces[i];
 
             // Append owner
-            if (n == allocSize_)
-            {
-                labels_.size() = n;
-                allocSize_ = allocSize_*2 + 1;
-                labels_.setSize(allocSize_);
-            }
-            labels_[n++] = own[faceI];
+            storage.append(own[faceI]);
 
             // Append neighbour
             if (faceI < nInternalFaces())
             {
-                if (n == allocSize_)
-                {
-                    labels_.size() = n;
-                    allocSize_ = allocSize_*2 + 1;
-                    labels_.setSize(allocSize_);
-                }
-                labels_[n++] = nei[faceI];
+                storage.append(nei[faceI]);
             }
         }
-        labels_.size() = n;
-
 
         // Filter duplicates
-        sort(labels_);
+        sort(storage);
 
-        n = 1;
+        label n = 1;
 
-        for (label i = 1; i < labels_.size(); i++)
+        for (label i = 1; i < storage.size(); i++)
         {
-            if (labels_[i] != labels_[i-1])
+            if (storage[i] != storage[i-1])
             {
-                labels_[n++] = labels_[i];
+                storage[n++] = storage[i];
             }
         }
 
-        labels_.size() = n;
+        storage.size() = n;
 
-        return labels_;
+        return storage;
     }
 }
 
 
+const labelList& primitiveMesh::pointCells(const label pointI) const
+{
+    return pointCells(pointI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointPoints.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointPoints.C
index 099e266941..e1ffa07dff 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointPoints.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointPoints.C
@@ -40,6 +40,14 @@ void primitiveMesh::calcPointPoints() const
         Pout<< "primitiveMesh::calcPointPoints() : "
             << "calculating pointPoints"
             << endl;
+
+        if (debug == -1)
+        {
+            // For checking calls:abort so we can quickly hunt down
+            // origin of call
+            FatalErrorIn("primitiveMesh::calcPointPoints()")
+                << abort(FatalError);
+        }
     }
 
     // It is an error to attempt to recalculate pointPoints
@@ -97,7 +105,11 @@ const labelListList& primitiveMesh::pointPoints() const
 }
 
 
-const labelList& primitiveMesh::pointPoints(const label pointI) const
+const labelList& primitiveMesh::pointPoints
+(
+    const label pointI,
+    DynamicList<label>& storage
+) const
 {
     if (hasPointPoints())
     {
@@ -108,30 +120,29 @@ const labelList& primitiveMesh::pointPoints(const label pointI) const
         const edgeList& edges = this->edges();
         const labelList& pEdges = pointEdges()[pointI];
 
-        labels_.size() = allocSize_;
+        storage.clear();
 
-        if (pEdges.size() > allocSize_)
+        if (pEdges.size() > storage.allocSize())
         {
-            // Set size() so memory allocation behaves as normal.
-            labels_.clear();
-            allocSize_ = pEdges.size();
-            labels_.setSize(allocSize_);
+            storage.setSize(pEdges.size());
         }
 
-        label n = 0;
-
         forAll(pEdges, i)
         {
-            labels_[n++] = edges[pEdges[i]].otherVertex(pointI);
+            storage.append(edges[pEdges[i]].otherVertex(pointI));
         }
 
-        labels_.size() = n;
-
-        return labels_;
+        return storage;
     }
 }
 
 
+const labelList& primitiveMesh::pointPoints(const label pointI) const
+{
+    return pointPoints(pointI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index 520e6f69ef..d7dda063b9 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -1075,11 +1075,8 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
                 )
             );
 
-            pointField zonePoints(zonePointIndices.size());
-            forAll(zonePointIndices, i)
-            {
-                zonePoints[i] = localPoints[zonePointIndices[i]];
-            }
+            pointField zonePoints(localPoints, zonePointIndices);
+            scalarField zoneSnapDist(snapDist, zonePointIndices);
 
             // Find nearest for points both on faceZone and pp.
             List<pointIndexHit> hitInfo;
@@ -1088,17 +1085,18 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
             (
                 labelList(1, zoneSurfI),
                 zonePoints,
-                sqr(4*snapDist),
+                sqr(4*zoneSnapDist),
                 hitSurface,
                 hitInfo
             );
 
-            forAll(hitInfo, pointI)
+            forAll(hitInfo, i)
             {
-                if (hitInfo[pointI].hit())
+                label pointI = zonePointIndices[i]; 
+                if (hitInfo[i].hit())
                 {
                     patchDisp[pointI] =
-                        hitInfo[pointI].hitPoint()
+                        hitInfo[i].hitPoint()
                       - localPoints[pointI];
                 }
                 else
@@ -1107,7 +1105,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
                         << "For point:" << pointI
                         << " coordinate:" << localPoints[pointI]
                         << " did not find any surface within:"
-                        << 4*snapDist[pointI]
+                        << 4*zoneSnapDist[i]
                         << " meter." << endl;
                 }
             }
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 58c74d3fb4..27366a7a7a 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -514,7 +514,7 @@ void Foam::meshRefinement::markBoundaryFace
 {
     isBoundaryFace[faceI] = true;
 
-    const labelList& fEdges = mesh_.faceEdges()[faceI];
+    const labelList& fEdges = mesh_.faceEdges(faceI);
 
     forAll(fEdges, fp)
     {
@@ -623,12 +623,16 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     // If so what is the remaining non-boundary anchor point?
     labelHashSet nonBoundaryAnchors(mesh_.nCells()/10000);
 
+    // On-the-fly addressing storage.
+    DynamicList<label> dynFEdges;
+    DynamicList<label> dynCPoints;
+
     // Count of faces marked for baffling
     label nBaffleFaces = 0;
 
     forAll(cellLevel, cellI)
     {
-        const labelList cPoints(meshCutter_.cellPoints(cellI));
+        const labelList& cPoints = mesh_.cellPoints(cellI, dynCPoints);
 
         // Get number of anchor points (pointLevel == cellLevel)
 
@@ -714,11 +718,14 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     // Loop over all points. If a point is connected to 4 or more cells
     // with 7 anchor points on the boundary set those cell's non-boundary faces
     // to baffles
+
+    DynamicList<label> dynPCells;
+
     forAllConstIter(labelHashSet, nonBoundaryAnchors, iter)
     {
         label pointI = iter.key();
 
-        const labelList& pCells = mesh_.pointCells()[pointI];
+        const labelList& pCells = mesh_.pointCells(pointI, dynPCells);
 
         // Count number of 'hasSevenBoundaryAnchorPoints' cells.
         label n = 0;
@@ -806,7 +813,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     {
         if (facePatch[faceI] == -1)
         {
-            const labelList& fEdges = mesh_.faceEdges()[faceI];
+            const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
             label nFaceBoundaryEdges = 0;
 
             forAll(fEdges, fe)
@@ -840,7 +847,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
             {
                 if (facePatch[faceI] == -1)
                 {
-                    const labelList& fEdges = mesh_.faceEdges()[faceI];
+                    const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
                     label nFaceBoundaryEdges = 0;
 
                     forAll(fEdges, fe)
@@ -1239,6 +1246,7 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
     labelList nBafflesPerEdge(mesh_.nEdges(), 0);
 
 
+
     // Count number of boundary faces per edge
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -1255,7 +1263,7 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
 
             forAll(pp, i)
             {
-                const labelList& fEdges = mesh_.faceEdges()[faceI];
+                const labelList& fEdges = mesh_.faceEdges(faceI);
 
                 forAll(fEdges, fEdgeI)
                 {
@@ -1267,19 +1275,23 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
     }
 
 
+    DynamicList<label> fe0;
+    DynamicList<label> fe1;
+
+
     // Count number of duplicate boundary faces per edge
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     forAll(couples, i)
     {
-        const labelList& fEdges0 = mesh_.faceEdges()[couples[i].first()];
+        const labelList& fEdges0 = mesh_.faceEdges(couples[i].first(), fe0);
 
         forAll(fEdges0, fEdgeI)
         {
             nBafflesPerEdge[fEdges0[fEdgeI]] += 1000000;
         }
 
-        const labelList& fEdges1 = mesh_.faceEdges()[couples[i].second()];
+        const labelList& fEdges1 = mesh_.faceEdges(couples[i].second(), fe1);
 
         forAll(fEdges1, fEdgeI)
         {
@@ -1314,7 +1326,7 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
          == patches.whichPatch(couple.second())
         )
         {
-            const labelList& fEdges = mesh_.faceEdges()[couples[i].first()];
+            const labelList& fEdges = mesh_.faceEdges(couples[i].first());
 
             forAll(fEdges, fEdgeI)
             {
diff --git a/src/dynamicMesh/motionSmoother/motionSmoother.C b/src/dynamicMesh/motionSmoother/motionSmoother.C
index 5b7dd0c3a7..e82103bdb9 100644
--- a/src/dynamicMesh/motionSmoother/motionSmoother.C
+++ b/src/dynamicMesh/motionSmoother/motionSmoother.C
@@ -343,7 +343,7 @@ void Foam::motionSmoother::getAffectedFacesAndPoints
 
         forAllConstIter(pointSet, nbrPoints, iter)
         {
-            const labelList& pCells = mesh_.pointCells()[iter.key()];  
+            const labelList& pCells = mesh_.pointCells(iter.key());
 
             forAll(pCells, pCellI)
             {
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C
index 0742fd44b5..232e49e0e7 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C
@@ -330,15 +330,14 @@ Foam::label Foam::addPatchCellLayer::addSideFace
     const label meshEdgeI,              // corresponding mesh edge
     const label layerI,                 // layer
     const label numEdgeFaces,           // number of layers for edge
+    const labelList& meshFaces,         // precalculated edgeFaces
     polyTopoChange& meshMod
 ) const
 {
     // Edge to 'inflate' from
     label inflateEdgeI = -1;
 
-    // Mesh faces using edge
-    const labelList& meshFaces = mesh_.edgeFaces()[meshEdgeI];
-
+    // Check mesh faces using edge
     forAll(meshFaces, i)
     {
         if (mesh_.isInternalFace(meshFaces[i]))
@@ -620,6 +619,9 @@ void Foam::addPatchCellLayer::setRefinement
 
     const labelList& meshPoints = pp.meshPoints();
 
+    // Some storage for edge-face-addressing.
+    DynamicList<label> ef;
+
     // Precalculate mesh edges for pp.edges.
     labelList meshEdges(calcMeshEdges(mesh_, pp));
 
@@ -777,7 +779,9 @@ void Foam::addPatchCellLayer::setRefinement
                 label meshEdgeI = meshEdges[edgeI];
 
                 // Mesh faces using edge
-                const labelList& meshFaces = mesh_.edgeFaces()[meshEdgeI];
+
+                // Mesh faces using edge
+                const labelList& meshFaces = mesh_.edgeFaces(meshEdgeI, ef);
 
                 // Check that there is only one patchface using edge.
                 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
@@ -1353,6 +1357,12 @@ void Foam::addPatchCellLayer::setRefinement
                             patchFaceI
                         );
 
+                        const labelList& meshFaces = mesh_.edgeFaces
+                        (
+                            meshEdgeI,
+                            ef
+                        );
+
                         addSideFace
                         (
                             pp,
@@ -1365,6 +1375,7 @@ void Foam::addPatchCellLayer::setRefinement
                             meshEdgeI,      // corresponding mesh edge
                             i,
                             numEdgeSideFaces,
+                            meshFaces,
                             meshMod
                         );
                     }
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H
index 0b5caaba86..d1a1162e96 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H
@@ -232,6 +232,7 @@ class addPatchCellLayer
             const label meshEdgeI,
             const label layerI,
             const label numEdgeFaces,
+            const labelList& meshFaces,
             polyTopoChange&
         ) const;
 
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C
index ca059c2d1e..74e21817d2 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C
@@ -125,11 +125,11 @@ void Foam::combineFaces::regioniseFaces
 (
     const scalar minCos,
     const label cellI,
+    const labelList& cEdges,
     Map<label>& faceRegion
 ) const
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-    const labelList& cEdges = mesh_.cellEdges()[cellI];
 
     forAll(cEdges, i)
     {
@@ -220,9 +220,10 @@ bool Foam::combineFaces::faceNeighboursValid
         return true;
     }
 
-    const labelListList& faceEdges = mesh_.faceEdges();
     const cell& cFaces = mesh_.cells()[cellI];
 
+    DynamicList<label> storage;
+
     // Test for face collapsing to edge since too many neighbours merged.
     forAll(cFaces, cFaceI)
     {
@@ -230,7 +231,7 @@ bool Foam::combineFaces::faceNeighboursValid
 
         if (!faceRegion.found(faceI))
         {
-            const labelList& fEdges = faceEdges[faceI];
+            const labelList& fEdges = mesh_.faceEdges(faceI, storage);
 
             // Count number of remaining faces neighbouring faceI. This has
             // to be 3 or more.
@@ -299,6 +300,8 @@ Foam::labelListList Foam::combineFaces::getMergeSets
 {
     // Lists of faces that can be merged.
     DynamicList<labelList> allFaceSets(boundaryCells.size() / 10);
+    // Storage for on-the-fly cell-edge addressing.
+    DynamicList<label> storage;
 
     // On all cells regionise the faces
     forAllConstIter(labelHashSet, boundaryCells, iter)
@@ -307,9 +310,11 @@ Foam::labelListList Foam::combineFaces::getMergeSets
 
         const cell& cFaces = mesh_.cells()[cellI];
 
+        const labelList& cEdges = mesh_.cellEdges(cellI, storage);
+
         // Region per face
         Map<label> faceRegion(cFaces.size());
-        regioniseFaces(featureCos, cellI, faceRegion);
+        regioniseFaces(featureCos, cellI, cEdges, faceRegion);
 
         // Now we have in faceRegion for every face the region with planar
         // face sharing the same region. We now check whether the resulting
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.H
index 208407dff3..a65db41e44 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.H
@@ -103,6 +103,7 @@ class combineFaces
         (
             const scalar minCos,
             const label cellI,
+            const labelList& cEdges,
             Map<label>& faceRegion
         ) const;
 
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
index 2a38afad1a..c69f6c4533 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
@@ -372,7 +372,7 @@ Foam::scalar Foam::hexRef8::getLevel0EdgeLength() const
         {
             const label cLevel = cellLevel_[cellI];
 
-            const labelList& cEdges = mesh_.cellEdges()[cellI];
+            const labelList& cEdges = mesh_.cellEdges(cellI);
 
             forAll(cEdges, i)
             {
@@ -447,7 +447,7 @@ Foam::scalar Foam::hexRef8::getLevel0EdgeLength() const
     {
         const label cLevel = cellLevel_[cellI];
 
-        const labelList& cEdges = mesh_.cellEdges()[cellI];
+        const labelList& cEdges = mesh_.cellEdges(cellI);
 
         forAll(cEdges, i)
         {
@@ -1190,6 +1190,10 @@ void Foam::hexRef8::createInternalFaces
     // From edge mid to face mids
     Map<edge> midPointToFaceMids(24);
 
+    // Storage for on-the-fly addressing
+    DynamicList<label> storage;
+
+
     // Running count of number of internal faces added so far.
     label nFacesAdded = 0;
 
@@ -1198,7 +1202,7 @@ void Foam::hexRef8::createInternalFaces
         label faceI = cFaces[i];
 
         const face& f = mesh_.faces()[faceI];
-        const labelList& fEdges = mesh_.faceEdges()[faceI];
+        const labelList& fEdges = mesh_.faceEdges(faceI, storage);
 
         // We are on the cellI side of face f. The face will have 1 or 4
         // cLevel points and lots of higher numbered ones.
@@ -1299,7 +1303,7 @@ void Foam::hexRef8::createInternalFaces
                     {
                         dumpCell(cellI);
 
-                        const labelList cPoints(cellPoints(cellI));
+                        const labelList& cPoints = mesh_.cellPoints(cellI);
 
                         FatalErrorIn("createInternalFaces(..)")
                             << "cell:" << cellI << " cLevel:" << cLevel
@@ -1372,7 +1376,7 @@ void Foam::hexRef8::createInternalFaces
                     {
                         dumpCell(cellI);
 
-                        const labelList cPoints(cellPoints(cellI));
+                        const labelList& cPoints = mesh_.cellPoints(cellI);
 
                         FatalErrorIn("createInternalFaces(..)")
                             << "cell:" << cellI << " cLevel:" << cLevel
@@ -1454,7 +1458,7 @@ void Foam::hexRef8::walkFaceToMid
 ) const
 {
     const face& f = mesh_.faces()[faceI];
-    const labelList& fEdges = mesh_.faceEdges()[faceI];
+    const labelList& fEdges = mesh_.faceEdges(faceI);
 
     label fp = startFp;
 
@@ -1503,7 +1507,7 @@ void Foam::hexRef8::walkFaceFromMid
 ) const
 {
     const face& f = mesh_.faces()[faceI];
-    const labelList& fEdges = mesh_.faceEdges()[faceI];
+    const labelList& fEdges = mesh_.faceEdges(faceI);
 
     label fp = f.rcIndex(startFp);
 
@@ -2013,27 +2017,6 @@ Foam::hexRef8::hexRef8
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-//- Get points of a cell (without using cellPoints addressing)
-Foam::labelList Foam::hexRef8::cellPoints(const label cellI) const
-{
-    // Pick up points of the cell
-    const cell& cFaces = mesh_.cells()[cellI];
-
-    labelHashSet cPoints(4*cFaces.size());
-
-    forAll(cFaces, i)
-    {
-        const face& f = mesh_.faces()[cFaces[i]];
-
-        forAll(f, fp)
-        {
-            cPoints.insert(f[fp]);
-        }
-    }
-    return cPoints.toc();
-}
-
-
 Foam::labelList Foam::hexRef8::consistentRefinement
 (
     const labelList& cellsToRefine,
@@ -2358,13 +2341,11 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
         // as cell level purely for ease)
         labelList maxPointCount(mesh_.nPoints(), 0);
 
-        const labelListList& pointCells = mesh_.pointCells();
-
-        forAll(pointCells, pointI)
+        forAll(maxPointCount, pointI)
         {
             label& pLevel = maxPointCount[pointI];
 
-            const labelList& pCells = pointCells[pointI];
+            const labelList& pCells = mesh_.pointCells(pointI);
 
             forAll(pCells, i)
             {
@@ -2395,7 +2376,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
             // Loop over all cells using the point and check whether their
             // refinement level is much less than the maximum.
 
-            const labelList& pCells = pointCells[pointI];
+            const labelList& pCells = mesh_.pointCells(pointI);
 
             forAll(pCells, pCellI)
             {
@@ -3121,7 +3102,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
     {
         if (cellMidPoint[cellI] >= 0)
         {
-            const labelList& cEdges = mesh_.cellEdges()[cellI];
+            const labelList& cEdges = mesh_.cellEdges(cellI);
 
             forAll(cEdges, i)
             {
@@ -3458,7 +3439,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
 
         forAll(pointLevel_, pointI)
         {
-            const labelList& pCells = mesh_.pointCells()[pointI];
+            const labelList& pCells = mesh_.pointCells(pointI);
 
             forAll(pCells, pCellI)
             {
@@ -3498,7 +3479,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
                 {
                     dumpCell(cellI);
 
-                    const labelList cPoints(cellPoints(cellI));
+                    const labelList& cPoints = mesh_.cellPoints(cellI);
 
                     FatalErrorIn
                     (
@@ -3610,7 +3591,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
         {
             if (edgeMidPoint[edgeI] >= 0)
             {
-                const labelList& eFaces = mesh_.edgeFaces()[edgeI];
+                const labelList& eFaces = mesh_.edgeFaces(edgeI);
 
                 forAll(eFaces, i)
                 {
@@ -3768,13 +3749,16 @@ Foam::labelListList Foam::hexRef8::setRefinement
             << endl;
     }
 
+    DynamicList<label> eFacesStorage;
+    DynamicList<label> fEdgesStorage;
+
     forAll(edgeMidPoint, edgeI)
     {
         if (edgeMidPoint[edgeI] >= 0)
         {
             // Split edge. Check that face not already handled above.
 
-            const labelList& eFaces = mesh_.edgeFaces()[edgeI];
+            const labelList& eFaces = mesh_.edgeFaces(edgeI, eFacesStorage);
 
             forAll(eFaces, i)
             {
@@ -3785,7 +3769,11 @@ Foam::labelListList Foam::hexRef8::setRefinement
                     // Unsplit face. Add edge splits to face.
 
                     const face& f = mesh_.faces()[faceI];
-                    const labelList& fEdges = mesh_.faceEdges()[faceI];
+                    const labelList& fEdges = mesh_.faceEdges
+                    (
+                        faceI,
+                        fEdgesStorage
+                    );
 
                     DynamicList<label> newFaceVerts(f.size());
 
@@ -4715,14 +4703,12 @@ void Foam::hexRef8::checkRefinementLevels
     // Check 2:1 across points (instead of faces)
     if (maxPointDiff != -1)
     {
-        const labelListList& pointCells = mesh_.pointCells();
-
         // Determine per point the max cell level.
         labelList maxPointLevel(mesh_.nPoints(), 0);
 
-        forAll(pointCells, pointI)
+        forAll(maxPointLevel, pointI)
         {
-            const labelList& pCells = pointCells[pointI];
+            const labelList& pCells = mesh_.pointCells(pointI);
 
             label& pLevel = maxPointLevel[pointI];
 
@@ -4747,7 +4733,7 @@ void Foam::hexRef8::checkRefinementLevels
         {
             label pointI = pointsToCheck[i];
 
-            const labelList& pCells = pointCells[pointI];
+            const labelList& pCells = mesh_.pointCells(pointI);
 
             forAll(pCells, i)
             {
@@ -4881,11 +4867,11 @@ Foam::labelList Foam::hexRef8::getSplitPoints() const
     labelList splitMasterLevel(mesh_.nPoints(), 0);
 
     // Unmark all with not 8 cells
-    const labelListList& pointCells = mesh_.pointCells();
+    //const labelListList& pointCells = mesh_.pointCells();
 
-    forAll(pointCells, pointI)
+    for (label pointI = 0; pointI < mesh_.nPoints(); pointI++)
     {
-        const labelList& pCells = pointCells[pointI];
+        const labelList& pCells = mesh_.pointCells(pointI);
 
         if (pCells.size() != 8)
         {
@@ -4898,8 +4884,7 @@ Foam::labelList Foam::hexRef8::getSplitPoints() const
 
     forAll(visibleCells, cellI)
     {
-        //const labelList& cPoints = mesh_.cellPoints()[cellI];
-        const labelList cPoints(cellPoints(cellI));
+        const labelList& cPoints = mesh_.cellPoints(cellI);
 
         if (visibleCells[cellI] != -1 && history_.parentIndex(cellI) >= 0)
         {
@@ -5104,7 +5089,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
         {
             if (unrefinePoint.get(pointI) == 1)
             {
-                const labelList& pCells = mesh_.pointCells()[pointI];
+                const labelList& pCells = mesh_.pointCells(pointI);
 
                 forAll(pCells, j)
                 {
@@ -5244,7 +5229,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
         {
             if (unrefinePoint.get(pointI) == 1)
             {
-                const labelList& pCells = mesh_.pointCells()[pointI];
+                const labelList& pCells = mesh_.pointCells(pointI);
 
                 forAll(pCells, j)
                 {
@@ -5329,7 +5314,7 @@ void Foam::hexRef8::setUnrefinement
 
         forAll(splitPointLabels, i)
         {
-            const labelList& pCells = mesh_.pointCells()[splitPointLabels[i]];
+            const labelList& pCells = mesh_.pointCells(splitPointLabels[i]);
 
             forAll(pCells, j)
             {
@@ -5395,7 +5380,7 @@ void Foam::hexRef8::setUnrefinement
 
         // Get original cell label
 
-        const labelList& pCells = mesh_.pointCells()[pointI];
+        const labelList& pCells = mesh_.pointCells(pointI);
 
         // Check
         if (pCells.size() != 8)
@@ -5463,7 +5448,7 @@ void Foam::hexRef8::setUnrefinement
     {
         label pointI = splitPointLabels[i];
 
-        const labelList& pCells = mesh_.pointCells()[pointI];
+        const labelList& pCells = mesh_.pointCells(pointI);
 
         label masterCellI = min(pCells);
 
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H
index faccfdc271..ea67874474 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H
@@ -370,9 +370,6 @@ public:
 
         // Refinement
 
-            //- Helper:get points of a cell without using cellPoints addressing
-            labelList cellPoints(const label cellI) const;
-
             //- Given valid mesh and current cell level and proposed
             //  cells to refine calculate any clashes (due to 2:1) and return
             //  ok list of cells to refine.
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C
index fab2b74609..32cb3c5b82 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C
@@ -83,6 +83,7 @@ Foam::label Foam::removeFaces::changeFaceRegion
     const labelList& nFacesPerEdge,
     const label faceI,
     const label newRegion,
+    const labelList& fEdges,
     labelList& faceRegion
 ) const
 {
@@ -94,27 +95,33 @@ Foam::label Foam::removeFaces::changeFaceRegion
 
         nChanged = 1;
 
+        // Storage for on-the-fly addressing
+        DynamicList<label> fe;
+        DynamicList<label> ef;
+
         // Step to neighbouring faces across edges that will get removed
-
-        const labelList& fEdges = mesh_.faceEdges()[faceI];
-
         forAll(fEdges, i)
         {
             label edgeI = fEdges[i];
 
             if (nFacesPerEdge[edgeI] >= 0 && nFacesPerEdge[edgeI] <= 2)
             {
-                const labelList& eFaces = mesh_.edgeFaces()[edgeI];
+                const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
 
                 forAll(eFaces, j)
                 {
+                    label nbrFaceI = eFaces[j];
+
+                    const labelList& fEdges1 = mesh_.faceEdges(nbrFaceI, fe);
+
                     nChanged += changeFaceRegion
                     (
                         cellRegion,
                         removedFace,
                         nFacesPerEdge,
-                        eFaces[j],
+                        nbrFaceI,
                         newRegion,
+                        fEdges1,
                         faceRegion
                     );
                 }
@@ -166,7 +173,7 @@ Foam::boolList Foam::removeFaces::getFacesAffected
     //  Mark faces affected by removal of edges
     forAllConstIter(labelHashSet, edgesToRemove, iter)
     {
-        const labelList& eFaces = mesh_.edgeFaces()[iter.key()];
+        const labelList& eFaces = mesh_.edgeFaces(iter.key());
 
         forAll(eFaces, eFaceI)
         {
@@ -814,6 +821,10 @@ void Foam::removeFaces::setRefinement
     // Number of connected face regions
     label nRegions = 0;
 
+    // Storage for on-the-fly addressing
+    DynamicList<label> fe;
+    DynamicList<label> ef;
+
 
     {
         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
@@ -827,7 +838,7 @@ void Foam::removeFaces::setRefinement
         {
             label faceI = faceLabels[i];
 
-            const labelList& fEdges = mesh_.faceEdges()[faceI];
+            const labelList& fEdges = mesh_.faceEdges(faceI, fe);
 
             forAll(fEdges, i)
             {
@@ -835,8 +846,7 @@ void Foam::removeFaces::setRefinement
 
                 if (nFacesPerEdge[edgeI] == -1)
                 {
-                    nFacesPerEdge[edgeI] =
-                        mesh_.edgeFaces()[edgeI].size()-1;
+                    nFacesPerEdge[edgeI] = mesh_.edgeFaces(edgeI, ef).size()-1;
                 }
                 else
                 {
@@ -849,16 +859,15 @@ void Foam::removeFaces::setRefinement
         // Note that this only needs to be done for possibly coupled edges
         // so we could choose to loop only over boundary faces and use faceEdges
         // of those.
-        const labelListList& edgeFaces = mesh_.edgeFaces();
 
-        forAll(edgeFaces, edgeI)
+        forAll(mesh_.edges(), edgeI)
         {
             if (nFacesPerEdge[edgeI] == -1)
             {
                 // Edge not yet handled in loop above so is not used by any
                 // face to be removed.
 
-                const labelList& eFaces = edgeFaces[edgeI];
+                const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
 
                 if (eFaces.size() > 2)
                 {
@@ -922,7 +931,7 @@ void Foam::removeFaces::setRefinement
                 label f0 = -1;
                 label f1 = -1;
 
-                const labelList& eFaces = mesh_.edgeFaces()[edgeI];
+                const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
 
                 forAll(eFaces, i)
                 {
@@ -1152,6 +1161,7 @@ void Foam::removeFaces::setRefinement
                 nFacesPerEdge,
                 startFaceI,
                 nRegions,
+                mesh_.faceEdges(startFaceI, fe),
                 faceRegion
             );
 
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.H
index 83186bed8b..48ac6169e8 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.H
@@ -93,6 +93,7 @@ class removeFaces
             const labelList& nFacesPerEdge,
             const label faceI,
             const label newRegion,
+            const labelList& fEdges,
             labelList& faceRegion
         ) const;