From 00b94c4a406e2c0d589f682e2934f1b57a1651f8 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Fri, 6 Oct 2017 14:18:54 +0200
Subject: [PATCH] ENH: handle partially attached master/slave faces (issue
 #608)

- addresses problems if faces are already partly attached
  (eg, on a single vertex or along an edge).
---
 .../polyTopoChange/polyTopoChange.C           | 521 ++++++++--------
 .../polyTopoChange/polyTopoChange.H           | 155 +++--
 .../polyTopoChange/polyTopoChangeTemplates.C  |  38 +-
 .../slidingInterface/coupleSlidingInterface.C | 571 +++++++++---------
 .../decoupleSlidingInterface.C                | 132 ++--
 .../enrichedPatch/enrichedPatch.C             |  39 +-
 .../enrichedPatch/enrichedPatch.H             |  37 +-
 .../enrichedPatch/enrichedPatchCutFaces.C     |  64 +-
 .../enrichedPatch/enrichedPatchFaces.C        | 214 +++----
 .../enrichedPatch/enrichedPatchI.H            |  63 ++
 .../enrichedPatch/enrichedPatchMasterPoints.C |  72 +--
 .../enrichedPatch/enrichedPatchPointMap.C     |  29 +-
 .../enrichedPatch/enrichedPatchPointPoints.C  |  11 +-
 .../slidingInterfaceAttachedAddressing.C      |  84 ++-
 .../slidingInterfaceProjectPoints.C           | 295 +++++----
 .../mesh/stitchMesh/simple-cube1/Allclean     |   8 +
 .../mesh/stitchMesh/simple-cube1/Allmesh      |  19 +
 tutorials/mesh/stitchMesh/simple-cube1/Allrun |  12 +
 .../stitchMesh/simple-cube1/scripts/stitch    |  50 ++
 .../simple-cube1/system/blockMeshDict         | 132 ++++
 .../simple-cube1/system/controlDict           |  45 ++
 .../simple-cube1/system/createPatchDict       |  53 ++
 .../stitchMesh/simple-cube1/system/fvSchemes  |  55 ++
 .../stitchMesh/simple-cube1/system/fvSolution |  43 ++
 .../simple-cube1/system/topoSetDict           | 106 ++++
 .../simple-cube1/system/topoSetDict.patches   |  88 +++
 26 files changed, 1739 insertions(+), 1197 deletions(-)
 create mode 100644 src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchI.H
 create mode 100755 tutorials/mesh/stitchMesh/simple-cube1/Allclean
 create mode 100755 tutorials/mesh/stitchMesh/simple-cube1/Allmesh
 create mode 100755 tutorials/mesh/stitchMesh/simple-cube1/Allrun
 create mode 100755 tutorials/mesh/stitchMesh/simple-cube1/scripts/stitch
 create mode 100644 tutorials/mesh/stitchMesh/simple-cube1/system/blockMeshDict
 create mode 100644 tutorials/mesh/stitchMesh/simple-cube1/system/controlDict
 create mode 100644 tutorials/mesh/stitchMesh/simple-cube1/system/createPatchDict
 create mode 100644 tutorials/mesh/stitchMesh/simple-cube1/system/fvSchemes
 create mode 100644 tutorials/mesh/stitchMesh/simple-cube1/system/fvSolution
 create mode 100644 tutorials/mesh/stitchMesh/simple-cube1/system/topoSetDict
 create mode 100644 tutorials/mesh/stitchMesh/simple-cube1/system/topoSetDict.patches

diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
index 78e9b7b309..90c3c63520 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,22 +54,22 @@ namespace Foam
 // Renumber with special handling for merged items (marked with <-1)
 void Foam::polyTopoChange::renumberReverseMap
 (
-    const labelList& map,
+    const labelUList& oldToNew,
     DynamicList<label>& elems
 )
 {
     forAll(elems, elemI)
     {
-        label val = elems[elemI];
+        const label val = elems[elemI];
 
         if (val >= 0)
         {
-            elems[elemI] = map[val];
+            elems[elemI] = oldToNew[val];
         }
         else if (val < -1)
         {
-            label mergedVal = -val-2;
-            elems[elemI] = -map[mergedVal]-2;
+            const label mergedVal = -val-2;
+            elems[elemI] = -oldToNew[mergedVal]-2;
         }
     }
 }
@@ -77,52 +77,52 @@ void Foam::polyTopoChange::renumberReverseMap
 
 void Foam::polyTopoChange::renumber
 (
-    const labelList& map,
-    labelHashSet& elems
+    const labelUList& oldToNew,
+    labelHashSet& labels
 )
 {
-    labelHashSet newElems(elems.size());
+    labelHashSet newSet(labels.capacity());
 
-    forAllConstIter(labelHashSet, elems, iter)
+    for (const label val : labels)
     {
-        label newElem = map[iter.key()];
+        const label newVal = oldToNew[val];
 
-        if (newElem >= 0)
+        if (newVal >= 0)
         {
-            newElems.insert(newElem);
+            newSet.insert(newVal);
         }
     }
 
-    elems.transfer(newElems);
+    labels.transfer(newSet);
 }
 
 
 // Renumber and remove -1 elements.
 void Foam::polyTopoChange::renumberCompact
 (
-    const labelList& map,
+    const labelUList& oldToNew,
     labelList& elems
 )
 {
-    label newElemI = 0;
+    label nElem = 0;
 
-    forAll(elems, elemI)
+    for (const label val : elems)
     {
-        label newVal = map[elems[elemI]];
+        const label newVal = oldToNew[val];
 
         if (newVal != -1)
         {
-            elems[newElemI++] = newVal;
+            elems[nElem++] = newVal;
         }
     }
-    elems.setSize(newElemI);
+    elems.setSize(nElem);
 }
 
 
 void Foam::polyTopoChange::countMap
 (
-    const labelList& map,
-    const labelList& reverseMap,
+    const labelUList& map,
+    const labelUList& reverseMap,
     label& nAdd,
     label& nInflate,
     label& nMerge,
@@ -136,7 +136,7 @@ void Foam::polyTopoChange::countMap
 
     forAll(map, newCelli)
     {
-        label oldCelli = map[newCelli];
+        const label oldCelli = map[newCelli];
 
         if (oldCelli >= 0)
         {
@@ -164,7 +164,7 @@ void Foam::polyTopoChange::countMap
 
     forAll(reverseMap, oldCelli)
     {
-        label newCelli = reverseMap[oldCelli];
+        const label newCelli = reverseMap[oldCelli];
 
         if (newCelli >= 0)
         {
@@ -224,8 +224,8 @@ void Foam::polyTopoChange::writeMeshStats(const polyMesh& mesh, Ostream& os)
 
 void Foam::polyTopoChange::getMergeSets
 (
-    const labelList& reverseCellMap,
-    const labelList& cellMap,
+    const labelUList& reverseCellMap,
+    const labelUList& cellMap,
     List<objectMap>& cellsFromCells
 )
 {
@@ -234,7 +234,7 @@ void Foam::polyTopoChange::getMergeSets
 
     forAll(reverseCellMap, oldCelli)
     {
-        label newCelli = reverseCellMap[oldCelli];
+        const label newCelli = reverseCellMap[oldCelli];
 
         if (newCelli < -1)
         {
@@ -267,15 +267,15 @@ void Foam::polyTopoChange::getMergeSets
 
     forAll(reverseCellMap, oldCelli)
     {
-        label newCelli = reverseCellMap[oldCelli];
+        const label newCelli = reverseCellMap[oldCelli];
 
         if (newCelli < -1)
         {
-            label mergeCelli = -newCelli-2;
+            const label mergeCelli = -newCelli-2;
 
             // oldCelli was merged into mergeCelli
 
-            label setI = cellToMergeSet[mergeCelli];
+            const label setI = cellToMergeSet[mergeCelli];
 
             objectMap& mergeSet = cellsFromCells[setI];
 
@@ -305,9 +305,9 @@ void Foam::polyTopoChange::getMergeSets
 
 bool Foam::polyTopoChange::hasValidPoints(const face& f) const
 {
-    forAll(f, fp)
+    for (const label fp : f)
     {
-        if (f[fp] < 0 || f[fp] >= points_.size())
+        if (fp < 0 || fp >= points_.size())
         {
             return false;
         }
@@ -595,7 +595,7 @@ void Foam::polyTopoChange::makeCellCells
 // Handles removed cells. Returns number of remaining cells.
 Foam::label Foam::polyTopoChange::getCellOrder
 (
-    const CompactListList<label>& cellCellAddressing,
+    const CompactListList<label, labelList>& cellCellAddressing,
     labelList& oldToNew
 ) const
 {
@@ -716,8 +716,8 @@ Foam::label Foam::polyTopoChange::getCellOrder
 void Foam::polyTopoChange::getFaceOrder
 (
     const label nActiveFaces,
-    const labelList& cellFaces,
-    const labelList& cellFaceOffsets,
+    const labelUList& cellFaces,
+    const labelUList& cellFaceOffsets,
 
     labelList& oldToNew,
     labelList& patchSizes,
@@ -882,7 +882,7 @@ void Foam::polyTopoChange::getFaceOrder
 void Foam::polyTopoChange::reorderCompactFaces
 (
     const label newSize,
-    const labelList& oldToNew
+    const labelUList& oldToNew
 )
 {
     reorder(oldToNew, faces_);
@@ -1069,9 +1069,9 @@ void Foam::polyTopoChange::compact
             newPointi = nActivePoints;
         }
 
-        forAllConstIter(labelHashSet, retiredPoints_, iter)
+        for (const label pointi : retiredPoints_)
         {
-            localPointMap[iter.key()] = newPointi++;
+            localPointMap[pointi] = newPointi++;
         }
 
 
@@ -1282,7 +1282,7 @@ void Foam::polyTopoChange::compact
 Foam::labelList Foam::polyTopoChange::selectFaces
 (
     const primitiveMesh& mesh,
-    const labelList& faceLabels,
+    const labelUList& faceLabels,
     const bool internalFacesOnly
 )
 {
@@ -1331,7 +1331,7 @@ Foam::labelList Foam::polyTopoChange::selectFaces
 // label)
 void Foam::polyTopoChange::calcPatchPointMap
 (
-    const List<Map<label>>& oldPatchMeshPointMaps,
+    const UList<Map<label>>& oldPatchMeshPointMaps,
     const polyBoundaryMesh& boundary,
     labelListList& patchPointMap
 ) const
@@ -1353,19 +1353,11 @@ void Foam::polyTopoChange::calcPatchPointMap
             if (meshPoints[i] < pointMap_.size())
             {
                 // Check if old point was part of same patch
-                Map<label>::const_iterator ozmpmIter = oldMeshPointMap.find
+                curPatchPointRnb[i] = oldMeshPointMap.lookup
                 (
-                    pointMap_[meshPoints[i]]
+                    pointMap_[meshPoints[i]],
+                    -1
                 );
-
-                if (ozmpmIter != oldMeshPointMap.end())
-                {
-                    curPatchPointRnb[i] = ozmpmIter();
-                }
-                else
-                {
-                    curPatchPointRnb[i] = -1;
-                }
             }
             else
             {
@@ -1394,38 +1386,24 @@ void Foam::polyTopoChange::calcFaceInflationMaps
         label nFacesFromPoints = 0;
 
         // Collect all still existing faces connected to this point.
-        forAllConstIter(Map<label>, faceFromPoint_, iter)
+        forAllConstIters(faceFromPoint_, iter)
         {
-            label newFacei = iter.key();
+            const label facei = iter.key();
+            const label pointi = iter.object();
 
-            if (region_[newFacei] == -1)
-            {
-                // Get internal faces using point on old mesh
-                facesFromPoints[nFacesFromPoints++] = objectMap
+            // Get internal or patch faces using point on old mesh
+            const bool internal = (region_[facei] == -1);
+
+            facesFromPoints[nFacesFromPoints++] = objectMap
+            (
+                facei,
+                selectFaces
                 (
-                    newFacei,
-                    selectFaces
-                    (
-                        mesh,
-                        mesh.pointFaces()[iter()],
-                        true
-                    )
-                );
-            }
-            else
-            {
-                // Get patch faces using point on old mesh
-                facesFromPoints[nFacesFromPoints++] = objectMap
-                (
-                    newFacei,
-                    selectFaces
-                    (
-                        mesh,
-                        mesh.pointFaces()[iter()],
-                        false
-                    )
-                );
-            }
+                    mesh,
+                    mesh.pointFaces()[pointi],
+                    internal
+                )
+            );
         }
     }
 
@@ -1440,38 +1418,24 @@ void Foam::polyTopoChange::calcFaceInflationMaps
         label nFacesFromEdges = 0;
 
         // Collect all still existing faces connected to this edge.
-        forAllConstIter(Map<label>, faceFromEdge_, iter)
+        forAllConstIters(faceFromEdge_, iter)
         {
-            label newFacei = iter.key();
+            const label facei = iter.key();
+            const label edgei = iter.object();
 
-            if (region_[newFacei] == -1)
-            {
-                // Get internal faces using edge on old mesh
-                facesFromEdges[nFacesFromEdges++] = objectMap
+            // Get internal or patch faces using edge on old mesh
+            const bool internal = (region_[facei] == -1);
+
+            facesFromEdges[nFacesFromEdges++] = objectMap
+            (
+                facei,
+                selectFaces
                 (
-                    newFacei,
-                    selectFaces
-                    (
-                        mesh,
-                        mesh.edgeFaces(iter()),
-                        true
-                    )
-                );
-            }
-            else
-            {
-                // Get patch faces using edge on old mesh
-                facesFromEdges[nFacesFromEdges++] = objectMap
-                (
-                    newFacei,
-                    selectFaces
-                    (
-                        mesh,
-                        mesh.edgeFaces(iter()),
-                        false
-                    )
-                );
-            }
+                    mesh,
+                    mesh.edgeFaces(edgei),
+                    internal
+                )
+            );
         }
     }
 
@@ -1503,13 +1467,16 @@ void Foam::polyTopoChange::calcCellInflationMaps
     {
         label nCellsFromPoints = 0;
 
-        // Collect all still existing faces connected to this point.
-        forAllConstIter(Map<label>, cellFromPoint_, iter)
+        // Collect all still existing cells connected to this point.
+        forAllConstIters(cellFromPoint_, iter)
         {
+            const label celli = iter.key();
+            const label pointi = iter.object();
+
             cellsFromPoints[nCellsFromPoints++] = objectMap
             (
-                iter.key(),
-                mesh.pointCells()[iter()]
+                celli,
+                mesh.pointCells()[pointi]
             );
         }
     }
@@ -1521,13 +1488,16 @@ void Foam::polyTopoChange::calcCellInflationMaps
     {
         label nCellsFromEdges = 0;
 
-        // Collect all still existing faces connected to this point.
-        forAllConstIter(Map<label>, cellFromEdge_, iter)
+        // Collect all still existing cells connected to this edge.
+        forAllConstIters(cellFromEdge_, iter)
         {
+            const label celli = iter.key();
+            const label edgei = iter.object();
+
             cellsFromEdges[nCellsFromEdges++] = objectMap
             (
-                iter.key(),
-                mesh.edgeCells()[iter()]
+                celli,
+                mesh.edgeCells()[edgei]
             );
         }
     }
@@ -1542,9 +1512,10 @@ void Foam::polyTopoChange::calcCellInflationMaps
         labelList twoCells(2);
 
         // Collect all still existing faces connected to this point.
-        forAllConstIter(Map<label>, cellFromFace_, iter)
+        forAllConstIters(cellFromFace_, iter)
         {
-            label oldFacei = iter();
+            const label celli = iter.key();
+            const label oldFacei = iter.object();
 
             if (mesh.isInternalFace(oldFacei))
             {
@@ -1552,7 +1523,7 @@ void Foam::polyTopoChange::calcCellInflationMaps
                 twoCells[1] = mesh.faceNeighbour()[oldFacei];
                 cellsFromFaces[nCellsFromFaces++] = objectMap
                 (
-                    iter.key(),
+                    celli,
                     twoCells
                 );
             }
@@ -1560,7 +1531,7 @@ void Foam::polyTopoChange::calcCellInflationMaps
             {
                 cellsFromFaces[nCellsFromFaces++] = objectMap
                 (
-                    iter.key(),
+                    celli,
                     labelList(1, mesh.faceOwner()[oldFacei])
                 );
             }
@@ -1600,49 +1571,51 @@ void Foam::polyTopoChange::resetZones
 
         labelList nPoints(pointZones.size(), 0);
 
-        forAllConstIter(Map<label>, pointZone_, iter)
+        forAllConstIters(pointZone_, iter)
         {
-            label zoneI = iter();
+            const label pointi = iter.key();
+            const label zonei = iter.object();
 
-            if (zoneI < 0 || zoneI >= pointZones.size())
+            if (zonei < 0 || zonei >= pointZones.size())
             {
                 FatalErrorInFunction
-                    << "Illegal zoneID " << zoneI << " for point "
-                    << iter.key() << " coord " << mesh.points()[iter.key()]
+                    << "Illegal zoneID " << zonei << " for point "
+                    << pointi << " coord " << mesh.points()[pointi]
                     << abort(FatalError);
             }
-            nPoints[zoneI]++;
+            nPoints[zonei]++;
         }
 
         // Distribute points per zone
 
         labelListList addressing(pointZones.size());
-        forAll(addressing, zoneI)
+        forAll(addressing, zonei)
         {
-            addressing[zoneI].setSize(nPoints[zoneI]);
+            addressing[zonei].setSize(nPoints[zonei]);
         }
         nPoints = 0;
 
-        forAllConstIter(Map<label>, pointZone_, iter)
+        forAllConstIters(pointZone_, iter)
         {
-            label zoneI = iter();
+            const label pointi = iter.key();
+            const label zonei = iter.object();
 
-            addressing[zoneI][nPoints[zoneI]++] = iter.key();
+            addressing[zonei][nPoints[zonei]++] = pointi;
         }
         // Sort the addressing
-        forAll(addressing, zoneI)
+        forAll(addressing, zonei)
         {
-            stableSort(addressing[zoneI]);
+            stableSort(addressing[zonei]);
         }
 
         // So now we both have old zones and the new addressing.
         // Invert the addressing to get pointZoneMap.
-        forAll(addressing, zoneI)
+        forAll(addressing, zonei)
         {
-            const pointZone& oldZone = pointZones[zoneI];
-            const labelList& newZoneAddr = addressing[zoneI];
+            const pointZone& oldZone = pointZones[zonei];
+            const labelList& newZoneAddr = addressing[zonei];
 
-            labelList& curPzRnb = pointZoneMap[zoneI];
+            labelList& curPzRnb = pointZoneMap[zonei];
             curPzRnb.setSize(newZoneAddr.size());
 
             forAll(newZoneAddr, i)
@@ -1660,17 +1633,17 @@ void Foam::polyTopoChange::resetZones
 
         // Reset the addresing on the zone
         newMesh.pointZones().clearAddressing();
-        forAll(newMesh.pointZones(), zoneI)
+        forAll(newMesh.pointZones(), zonei)
         {
             if (debug)
             {
-                Pout<< "pointZone:" << zoneI
-                    << "  name:" << newMesh.pointZones()[zoneI].name()
-                    << "  size:" << addressing[zoneI].size()
+                Pout<< "pointZone:" << zonei
+                    << "  name:" << newMesh.pointZones()[zonei].name()
+                    << "  size:" << addressing[zonei].size()
                     << endl;
             }
 
-            newMesh.pointZones()[zoneI] = addressing[zoneI];
+            newMesh.pointZones()[zonei] = addressing[zonei];
         }
     }
 
@@ -1684,71 +1657,73 @@ void Foam::polyTopoChange::resetZones
 
         labelList nFaces(faceZones.size(), 0);
 
-        forAllConstIter(Map<label>, faceZone_, iter)
+        forAllConstIters(faceZone_, iter)
         {
-            label zoneI = iter();
+            const label facei = iter.key();
+            const label zonei = iter.object();
 
-            if (zoneI < 0 || zoneI >= faceZones.size())
+            if (zonei < 0 || zonei >= faceZones.size())
             {
                 FatalErrorInFunction
-                    << "Illegal zoneID " << zoneI << " for face "
-                    << iter.key()
+                    << "Illegal zoneID " << zonei << " for face "
+                    << facei
                     << abort(FatalError);
             }
-            nFaces[zoneI]++;
+            nFaces[zonei]++;
         }
 
         labelListList addressing(faceZones.size());
         boolListList flipMode(faceZones.size());
 
-        forAll(addressing, zoneI)
+        forAll(addressing, zonei)
         {
-            addressing[zoneI].setSize(nFaces[zoneI]);
-            flipMode[zoneI].setSize(nFaces[zoneI]);
+            addressing[zonei].setSize(nFaces[zonei]);
+            flipMode[zonei].setSize(nFaces[zonei]);
         }
         nFaces = 0;
 
-        forAllConstIter(Map<label>, faceZone_, iter)
+        forAllConstIters(faceZone_, iter)
         {
-            label zoneI = iter();
-            label facei = iter.key();
+            const label facei = iter.key();
+            const label zonei = iter.object();
 
-            label index = nFaces[zoneI]++;
+            const label index = nFaces[zonei]++;
 
-            addressing[zoneI][index] = facei;
-            flipMode[zoneI][index] = faceZoneFlip_[facei];
+            addressing[zonei][index] = facei;
+            flipMode[zonei][index] = faceZoneFlip_[facei];
         }
+
         // Sort the addressing
-        forAll(addressing, zoneI)
+        forAll(addressing, zonei)
         {
             labelList newToOld;
-            sortedOrder(addressing[zoneI], newToOld);
+            sortedOrder(addressing[zonei], newToOld);
             {
-                labelList newAddressing(addressing[zoneI].size());
+                labelList newAddressing(addressing[zonei].size());
                 forAll(newAddressing, i)
                 {
-                    newAddressing[i] = addressing[zoneI][newToOld[i]];
+                    newAddressing[i] = addressing[zonei][newToOld[i]];
                 }
-                addressing[zoneI].transfer(newAddressing);
+                addressing[zonei].transfer(newAddressing);
             }
             {
-                boolList newFlipMode(flipMode[zoneI].size());
+                boolList newFlipMode(flipMode[zonei].size());
                 forAll(newFlipMode, i)
                 {
-                    newFlipMode[i] = flipMode[zoneI][newToOld[i]];
+                    newFlipMode[i] = flipMode[zonei][newToOld[i]];
                 }
-                flipMode[zoneI].transfer(newFlipMode);
+                flipMode[zonei].transfer(newFlipMode);
             }
         }
 
         // So now we both have old zones and the new addressing.
         // Invert the addressing to get faceZoneFaceMap.
-        forAll(addressing, zoneI)
+        forAll(addressing, zonei)
         {
-            const faceZone& oldZone = faceZones[zoneI];
-            const labelList& newZoneAddr = addressing[zoneI];
+            const faceZone& oldZone = faceZones[zonei];
+            const labelList& newZoneAddr = addressing[zonei];
 
-            labelList& curFzFaceRnb = faceZoneFaceMap[zoneI];
+            labelList& curFzFaceRnb = faceZoneFaceMap[zonei];
 
             curFzFaceRnb.setSize(newZoneAddr.size());
 
@@ -1769,20 +1744,20 @@ void Foam::polyTopoChange::resetZones
 
         // Reset the addresing on the zone
         newMesh.faceZones().clearAddressing();
-        forAll(newMesh.faceZones(), zoneI)
+        forAll(newMesh.faceZones(), zonei)
         {
             if (debug)
             {
-                Pout<< "faceZone:" << zoneI
-                    << "  name:" << newMesh.faceZones()[zoneI].name()
-                    << "  size:" << addressing[zoneI].size()
+                Pout<< "faceZone:" << zonei
+                    << "  name:" << newMesh.faceZones()[zonei].name()
+                    << "  size:" << addressing[zonei].size()
                     << endl;
             }
 
-            newMesh.faceZones()[zoneI].resetAddressing
+            newMesh.faceZones()[zonei].resetAddressing
             (
-                addressing[zoneI],
-                flipMode[zoneI]
+                addressing[zonei],
+                flipMode[zonei]
             );
         }
     }
@@ -1799,51 +1774,51 @@ void Foam::polyTopoChange::resetZones
 
         forAll(cellZone_, celli)
         {
-            label zoneI = cellZone_[celli];
+            const label zonei = cellZone_[celli];
 
-            if (zoneI >= cellZones.size())
+            if (zonei >= cellZones.size())
             {
                 FatalErrorInFunction
-                    << "Illegal zoneID " << zoneI << " for cell "
+                    << "Illegal zoneID " << zonei << " for cell "
                     << celli << abort(FatalError);
             }
 
-            if (zoneI >= 0)
+            if (zonei >= 0)
             {
-                nCells[zoneI]++;
+                nCells[zonei]++;
             }
         }
 
         labelListList addressing(cellZones.size());
-        forAll(addressing, zoneI)
+        forAll(addressing, zonei)
         {
-            addressing[zoneI].setSize(nCells[zoneI]);
+            addressing[zonei].setSize(nCells[zonei]);
         }
         nCells = 0;
 
         forAll(cellZone_, celli)
         {
-            label zoneI = cellZone_[celli];
+            const label zonei = cellZone_[celli];
 
-            if (zoneI >= 0)
+            if (zonei >= 0)
             {
-                addressing[zoneI][nCells[zoneI]++] = celli;
+                addressing[zonei][nCells[zonei]++] = celli;
             }
         }
         // Sort the addressing
-        forAll(addressing, zoneI)
+        forAll(addressing, zonei)
         {
-            stableSort(addressing[zoneI]);
+            stableSort(addressing[zonei]);
         }
 
         // So now we both have old zones and the new addressing.
         // Invert the addressing to get cellZoneMap.
-        forAll(addressing, zoneI)
+        forAll(addressing, zonei)
         {
-            const cellZone& oldZone = cellZones[zoneI];
-            const labelList& newZoneAddr = addressing[zoneI];
+            const cellZone& oldZone = cellZones[zonei];
+            const labelList& newZoneAddr = addressing[zonei];
 
-            labelList& curCellRnb = cellZoneMap[zoneI];
+            labelList& curCellRnb = cellZoneMap[zonei];
 
             curCellRnb.setSize(newZoneAddr.size());
 
@@ -1863,17 +1838,17 @@ void Foam::polyTopoChange::resetZones
 
         // Reset the addresing on the zone
         newMesh.cellZones().clearAddressing();
-        forAll(newMesh.cellZones(), zoneI)
+        forAll(newMesh.cellZones(), zonei)
         {
             if (debug)
             {
-                Pout<< "cellZone:" << zoneI
-                    << "  name:" << newMesh.cellZones()[zoneI].name()
-                    << "  size:" << addressing[zoneI].size()
+                Pout<< "cellZone:" << zonei
+                    << "  name:" << newMesh.cellZones()[zonei].name()
+                    << "  size:" << addressing[zonei].size()
                     << endl;
             }
 
-            newMesh.cellZones()[zoneI] = addressing[zoneI];
+            newMesh.cellZones()[zonei] = addressing[zonei];
         }
     }
 }
@@ -1882,7 +1857,7 @@ void Foam::polyTopoChange::resetZones
 void Foam::polyTopoChange::calcFaceZonePointMap
 (
     const polyMesh& mesh,
-    const List<Map<label>>& oldFaceZoneMeshPointMaps,
+    const UList<Map<label>>& oldFaceZoneMeshPointMaps,
     labelListList& faceZonePointMap
 ) const
 {
@@ -1890,15 +1865,15 @@ void Foam::polyTopoChange::calcFaceZonePointMap
 
     faceZonePointMap.setSize(faceZones.size());
 
-    forAll(faceZones, zoneI)
+    forAll(faceZones, zonei)
     {
-        const faceZone& newZone = faceZones[zoneI];
+        const faceZone& newZone = faceZones[zonei];
 
         const labelList& newZoneMeshPoints = newZone().meshPoints();
 
-        const Map<label>& oldZoneMeshPointMap = oldFaceZoneMeshPointMaps[zoneI];
+        const Map<label>& oldZoneMeshPointMap = oldFaceZoneMeshPointMaps[zonei];
 
-        labelList& curFzPointRnb = faceZonePointMap[zoneI];
+        labelList& curFzPointRnb = faceZonePointMap[zonei];
 
         curFzPointRnb.setSize(newZoneMeshPoints.size());
 
@@ -1906,20 +1881,12 @@ void Foam::polyTopoChange::calcFaceZonePointMap
         {
             if (newZoneMeshPoints[pointi] < pointMap_.size())
             {
-                Map<label>::const_iterator ozmpmIter =
-                    oldZoneMeshPointMap.find
+                curFzPointRnb[pointi] =
+                    oldZoneMeshPointMap.lookup
                     (
-                        pointMap_[newZoneMeshPoints[pointi]]
+                        pointMap_[newZoneMeshPoints[pointi]],
+                        -1
                     );
-
-                if (ozmpmIter != oldZoneMeshPointMap.end())
-                {
-                    curFzPointRnb[pointi] = ozmpmIter();
-                }
-                else
-                {
-                    curFzPointRnb[pointi] = -1;
-                }
             }
             else
             {
@@ -1934,8 +1901,8 @@ void Foam::polyTopoChange::reorderCoupledFaces
 (
     const bool syncParallel,
     const polyBoundaryMesh& boundary,
-    const labelList& patchStarts,
-    const labelList& patchSizes,
+    const labelUList& patchStarts,
+    const labelUList& patchSizes,
     const pointField& points
 )
 {
@@ -1986,7 +1953,7 @@ void Foam::polyTopoChange::reorderCoupledFaces
             labelList patchFaceMap(patchSizes[patchi], -1);
             labelList patchFaceRotation(patchSizes[patchi], 0);
 
-            bool changed = boundary[patchi].order
+            const bool changed = boundary[patchi].order
             (
                 pBufs,
                 primitivePatch
@@ -2006,7 +1973,7 @@ void Foam::polyTopoChange::reorderCoupledFaces
             if (changed)
             {
                 // Merge patch face reordering into mesh face reordering table
-                label start = patchStarts[patchi];
+                const label start = patchStarts[patchi];
 
                 forAll(patchFaceMap, patchFacei)
                 {
@@ -2169,11 +2136,11 @@ void Foam::polyTopoChange::compactAndReorder
     // later on to calculate the faceZone pointMaps.
     oldFaceZoneMeshPointMaps.setSize(mesh.faceZones().size());
 
-    forAll(mesh.faceZones(), zoneI)
+    forAll(mesh.faceZones(), zonei)
     {
-        const faceZone& oldZone = mesh.faceZones()[zoneI];
+        const faceZone& oldZone = mesh.faceZones()[zonei];
 
-        oldFaceZoneMeshPointMaps[zoneI] = oldZone().meshPointMap();
+        oldFaceZoneMeshPointMaps[zonei] = oldZone().meshPointMap();
     }
 }
 
@@ -2290,16 +2257,16 @@ void Foam::polyTopoChange::clear()
 void Foam::polyTopoChange::addMesh
 (
     const polyMesh& mesh,
-    const labelList& patchMap,
-    const labelList& pointZoneMap,
-    const labelList& faceZoneMap,
-    const labelList& cellZoneMap
+    const labelUList& patchMap,
+    const labelUList& pointZoneMap,
+    const labelUList& faceZoneMap,
+    const labelUList& cellZoneMap
 )
 {
     label maxRegion = nPatches_ - 1;
-    forAll(patchMap, i)
+    for (const label regioni : patchMap)
     {
-        maxRegion = max(maxRegion, patchMap[i]);
+        maxRegion = max(maxRegion, regioni);
     }
     nPatches_ = maxRegion + 1;
 
@@ -2318,13 +2285,13 @@ void Foam::polyTopoChange::addMesh
         // Precalc offset zones
         labelList newZoneID(points.size(), -1);
 
-        forAll(pointZones, zoneI)
+        forAll(pointZones, zonei)
         {
-            const labelList& pointLabels = pointZones[zoneI];
+            const labelList& pointLabels = pointZones[zonei];
 
-            forAll(pointLabels, j)
+            for (const label pointi : pointLabels)
             {
-                newZoneID[pointLabels[j]] = pointZoneMap[zoneI];
+                newZoneID[pointi] = pointZoneMap[zonei];
             }
         }
 
@@ -2362,14 +2329,12 @@ void Foam::polyTopoChange::addMesh
         // Precalc offset zones
         labelList newZoneID(nAllCells, -1);
 
-        forAll(cellZones, zoneI)
+        forAll(cellZones, zonei)
         {
-            const labelList& cellLabels = cellZones[zoneI];
+            const labelList& cellLabels = cellZones[zonei];
 
-            forAll(cellLabels, j)
+            for (const label celli : cellLabels)
             {
-                label celli = cellLabels[j];
-
                 if (newZoneID[celli] != -1)
                 {
                     WarningInFunction
@@ -2377,13 +2342,13 @@ void Foam::polyTopoChange::addMesh
                         << " centre:" << mesh.cellCentres()[celli]
                         << " is in two zones:"
                         << cellZones[newZoneID[celli]].name()
-                        << " and " << cellZones[zoneI].name() << endl
+                        << " and " << cellZones[zonei].name() << endl
                         << "    This is not supported."
                         << " Continuing with first zone only." << endl;
                 }
                 else
                 {
-                    newZoneID[celli] = cellZoneMap[zoneI];
+                    newZoneID[celli] = cellZoneMap[zonei];
                 }
             }
         }
@@ -2424,15 +2389,15 @@ void Foam::polyTopoChange::addMesh
         labelList newZoneID(nAllFaces, -1);
         boolList zoneFlip(nAllFaces, false);
 
-        forAll(faceZones, zoneI)
+        forAll(faceZones, zonei)
         {
-            const labelList& faceLabels = faceZones[zoneI];
-            const boolList& flipMap = faceZones[zoneI].flipMap();
+            const labelList& faceLabels = faceZones[zonei];
+            const boolList& flipMap = faceZones[zonei].flipMap();
 
-            forAll(faceLabels, j)
+            forAll(faceLabels, facei)
             {
-                newZoneID[faceLabels[j]] = faceZoneMap[zoneI];
-                zoneFlip[faceLabels[j]] = flipMap[j];
+                newZoneID[faceLabels[facei]] = faceZoneMap[zonei];
+                zoneFlip[faceLabels[facei]] = flipMap[facei];
             }
         }
 
@@ -2473,7 +2438,7 @@ void Foam::polyTopoChange::addMesh
             }
             forAll(pp, patchFacei)
             {
-                label facei = pp.start() + patchFacei;
+                const label facei = pp.start() + patchFacei;
 
                 addFace
                 (
@@ -2663,7 +2628,7 @@ Foam::label Foam::polyTopoChange::addPoint
     const bool inCell
 )
 {
-    label pointi = points_.size();
+    const label pointi = points_.size();
 
     points_.append(pt);
     pointMap_.append(masterPointID);
@@ -2687,7 +2652,7 @@ void Foam::polyTopoChange::modifyPoint
 (
     const label pointi,
     const point& pt,
-    const label newZoneID,
+    const label zoneID,
     const bool inCell
 )
 {
@@ -2706,22 +2671,13 @@ void Foam::polyTopoChange::modifyPoint
     }
     points_[pointi] = pt;
 
-    Map<label>::iterator pointFnd = pointZone_.find(pointi);
-
-    if (pointFnd != pointZone_.end())
+    if (zoneID >= 0)
     {
-        if (newZoneID >= 0)
-        {
-            pointFnd() = newZoneID;
-        }
-        else
-        {
-            pointZone_.erase(pointFnd);
-        }
+        pointZone_.set(pointi, zoneID);
     }
-    else if (newZoneID >= 0)
+    else
     {
-        pointZone_.insert(pointi, newZoneID);
+        pointZone_.erase(pointi);
     }
 
     if (inCell)
@@ -2889,29 +2845,24 @@ void Foam::polyTopoChange::modifyFace
     region_[facei] = patchID;
 
     flipFaceFlux_[facei] = (flipFaceFlux ? 1 : 0);
-
-    Map<label>::iterator faceFnd = faceZone_.find(facei);
-
-    if (faceFnd != faceZone_.end())
-    {
-        if (zoneID >= 0)
-        {
-            faceFnd() = zoneID;
-        }
-        else
-        {
-            faceZone_.erase(faceFnd);
-        }
-    }
-    else if (zoneID >= 0)
-    {
-        faceZone_.insert(facei, zoneID);
-    }
     faceZoneFlip_[facei] = (zoneFlip ? 1 : 0);
+
+    if (zoneID >= 0)
+    {
+        faceZone_.set(facei, zoneID);
+    }
+    else
+    {
+        faceZone_.erase(facei);
+    }
 }
 
 
-void Foam::polyTopoChange::removeFace(const label facei, const label mergeFacei)
+void Foam::polyTopoChange::removeFace
+(
+    const label facei,
+    const label mergeFacei
+)
 {
     if (facei < 0 || facei >= faces_.size())
     {
@@ -2949,8 +2900,8 @@ void Foam::polyTopoChange::removeFace(const label facei, const label mergeFacei)
     faceFromEdge_.erase(facei);
     faceFromPoint_.erase(facei);
     flipFaceFlux_[facei] = 0;
-    faceZone_.erase(facei);
     faceZoneFlip_[facei] = 0;
+    faceZone_.erase(facei);
 }
 
 
@@ -3001,7 +2952,11 @@ void Foam::polyTopoChange::modifyCell
 }
 
 
-void Foam::polyTopoChange::removeCell(const label celli, const label mergeCelli)
+void Foam::polyTopoChange::removeCell
+(
+    const label celli,
+    const label mergeCelli
+)
 {
     if (celli < 0 || celli >= cellMap_.size())
     {
@@ -3124,7 +3079,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::changeMesh
 
         forAll(pointMap_, newPointi)
         {
-            label oldPointi = pointMap_[newPointi];
+            const label oldPointi = pointMap_[newPointi];
 
             if (oldPointi >= 0)
             {
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.H
index e97c6ef7db..f9d7640583 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.H
@@ -197,30 +197,57 @@ class polyTopoChange
 
     // Private Member Functions
 
-        //- Reorder contents of container according to map
+        //- Reorder contents of container according to oldToNew map
         template<class T>
-        static void reorder(const labelList& map, DynamicList<T>&);
-        template<class T>
-        static void reorder(const labelList& map, List<DynamicList<T>>&);
-        template<class T>
-        static void renumberKey(const labelList& map, Map<T>&);
+        static void reorder
+        (
+            const labelUList& oldToNew,
+            DynamicList<T>& lst
+        );
+
+        template<class T>
+        static void reorder
+        (
+            const labelUList& oldToNew,
+            List<DynamicList<T>>& lst
+        );
+
+        template<class T>
+        static void renumberKey
+        (
+            const labelUList& oldToNew,
+            Map<T>& map
+        );
+
+        //- Renumber elements of container according to oldToNew map
+        static void renumber
+        (
+            const labelUList& oldToNew,
+            labelHashSet& labels
+        );
 
-        //- Renumber elements of container according to map
-        static void renumber(const labelList&, labelHashSet&);
         //- Special handling of reverse maps which have <-1 in them
-        static void renumberReverseMap(const labelList&, DynamicList<label>&);
+        static void renumberReverseMap
+        (
+            const labelUList& oldToNew,
+            DynamicList<label>& elems
+        );
 
         //- Renumber & compact elements of list according to map
-        static void renumberCompact(const labelList&, labelList&);
+        static void renumberCompact
+        (
+            const labelUList& oldToNew,
+            labelList& elems
+        );
 
         //- Get all set elements as a labelHashSet
-        static labelHashSet getSetIndices(const PackedBoolList&);
+        static labelHashSet getSetIndices(const PackedBoolList& lst);
 
         //- Count number of added and removed quantities from maps.
         static void countMap
         (
-            const labelList& map,
-            const labelList& reverseMap,
+            const labelUList& map,
+            const labelUList& reverseMap,
             label& nAdd,
             label& nInflate,
             label& nMerge,
@@ -228,19 +255,19 @@ class polyTopoChange
         );
 
         //- Print some stats about mesh
-        static void writeMeshStats(const polyMesh& mesh, Ostream&);
+        static void writeMeshStats(const polyMesh& mesh, Ostream& os);
 
         //- Calculate object maps. Requires reverseMap to have destination
         //  to be marked with <-1.
         static void getMergeSets
         (
-            const labelList& reverseCellMap,
-            const labelList& cellMap,
+            const labelUList& reverseCellMap,
+            const labelUList& cellMap,
             List<objectMap>& cellsFromCells
         );
 
         //- Are all face vertices valid
-        bool hasValidPoints(const face&) const;
+        bool hasValidPoints(const face& f) const;
 
         //- Return face points
         pointField facePoints(const face& f) const;
@@ -248,7 +275,7 @@ class polyTopoChange
         //- Check inputs to modFace or addFace
         void checkFace
         (
-            const face&,
+            const face& f,
             const label facei,
             const label own,
             const label nei,
@@ -274,16 +301,16 @@ class polyTopoChange
         //- Cell ordering (bandCompression). Returns number of remaining cells.
         label getCellOrder
         (
-            const CompactListList<label, labelList>&,
-            labelList&
+            const CompactListList<label, labelList>& cellCellAddressing,
+            labelList& oldToNew
         ) const;
 
         //- Do upper-triangular ordering and patch ordering.
         void getFaceOrder
         (
             const label nActiveFaces,
-            const labelList& cellFaces,
-            const labelList& cellFaceOffsets,
+            const labelUList& cellFaces,
+            const labelUList& cellFaceOffsets,
 
             labelList& oldToNew,
             labelList& patchSizes,
@@ -294,7 +321,7 @@ class polyTopoChange
         void reorderCompactFaces
         (
             const label newSize,
-            const labelList& oldToNew
+            const labelUList& oldToNew
         );
 
         //- Remove all unused/removed points/faces/cells and update
@@ -314,50 +341,50 @@ class polyTopoChange
         //- Select either internal or external faces out of faceLabels
         static labelList selectFaces
         (
-            const primitiveMesh&,
-            const labelList& faceLabels,
+            const primitiveMesh& mesh,
+            const labelUList& faceLabels,
             const bool internalFacesOnly
         );
 
         //- Calculate mapping for patchpoints only
         void calcPatchPointMap
         (
-            const List<Map<label>>&,
-            const polyBoundaryMesh&,
-            labelListList&
+            const UList<Map<label>>& oldPatchMeshPointMaps,
+            const polyBoundaryMesh& boundary,
+            labelListList& patchPointMap
         ) const;
 
         void calcFaceInflationMaps
         (
-            const polyMesh&,
-            List<objectMap>&,
-            List<objectMap>&,
-            List<objectMap>&
+            const polyMesh& mesh,
+            List<objectMap>& facesFromPoints,
+            List<objectMap>& facesFromEdges,
+            List<objectMap>& facesFromFaces
         ) const;
 
         void calcCellInflationMaps
         (
-            const polyMesh&,
-            List<objectMap>&,
-            List<objectMap>&,
-            List<objectMap>&,
-            List<objectMap>&
+            const polyMesh& mesh,
+            List<objectMap>& cellsFromPoints,
+            List<objectMap>& cellsFromEdges,
+            List<objectMap>& cellsFromFaces,
+            List<objectMap>& cellsFromCells
         ) const;
 
         void resetZones
         (
-            const polyMesh&,        // mesh to get existing info from
-            polyMesh&,              // mesh to change zones on
-            labelListList&,
-            labelListList&,
-            labelListList&
+            const polyMesh& mesh,   // mesh to get existing info from
+            polyMesh& newMesh,      // mesh to change zones on
+            labelListList& pointZoneMap,
+            labelListList& faceZoneFaceMap,
+            labelListList& cellZoneMap
         ) const;
 
         void calcFaceZonePointMap
         (
-            const polyMesh&,
-            const List<Map<label>>&,
-            labelListList&
+            const polyMesh& mesh,
+            const UList<Map<label>>& oldFaceZoneMeshPointMaps,
+            labelListList& faceZonePointMap
         ) const;
 
 
@@ -367,15 +394,15 @@ class polyTopoChange
             void reorderCoupledFaces
             (
                 const bool syncParallel,
-                const polyBoundaryMesh&,
-                const labelList& patchStarts,
-                const labelList& patchSizes,
+                const polyBoundaryMesh& boundary,
+                const labelUList& patchStarts,
+                const labelUList& patchSizes,
                 const pointField& points
             );
 
         void compactAndReorder
         (
-            const polyMesh&,
+            const polyMesh& mesh,
             const bool syncParallel,
             const bool orderCells,
             const bool orderPoints,
@@ -445,12 +472,15 @@ public:
             }
 
             //- Is point removed?
+            //  Considered removed if point is GREAT.
             inline bool pointRemoved(const label pointi) const;
 
             //- Is face removed?
+            //  Considered removed if face is empty
             inline bool faceRemoved(const label facei) const;
 
             //- Is cell removed?
+            //  Considered removed if the cellMap is -2
             inline bool cellRemoved(const label celli) const;
 
 
@@ -463,11 +493,11 @@ public:
             //  or zone ids.
             void addMesh
             (
-                const polyMesh&,
-                const labelList& patchMap,
-                const labelList& pointZoneMap,
-                const labelList& faceZoneMap,
-                const labelList& cellZoneMap
+                const polyMesh& mesh,
+                const labelUList& patchMap,
+                const labelUList& pointZoneMap,
+                const labelUList& faceZoneMap,
+                const labelUList& cellZoneMap
             );
 
             //- Explicitly pre-size the dynamic storage for expected mesh
@@ -491,7 +521,7 @@ public:
             //  - inCell = false: add retired point (to end of point list)
             label addPoint
             (
-                const point&,
+                const point& pt,
                 const label masterPointID,
                 const label zoneID,
                 const bool inCell
@@ -499,17 +529,18 @@ public:
 
             //- Modify coordinate.
             //  Notes:
+            //  - zoneID = +ve (add to zoneID), -ve (remove from zones)
             //  - inCell = false: add retired point (to end of point list)
             void modifyPoint
             (
-                const label,
-                const point&,
-                const label newZoneID,
+                const label pointi,
+                const point& pt,
+                const label zoneID,
                 const bool inCell
             );
 
             //- Remove/merge point.
-            void removePoint(const label, const label);
+            void removePoint(const label pointi, const label mergePointi);
 
             //- Add face to cells. Return new face label.
             //  own,nei<0, zoneID>=0 : add inactive face (to end of face list)
@@ -541,7 +572,7 @@ public:
             );
 
             //- Remove/merge face.
-            void removeFace(const label, const label);
+            void removeFace(const label facei, const label mergeFacei);
 
             //- Add cell. Return new cell label.
             label addCell
@@ -554,10 +585,10 @@ public:
             );
 
             //- Modify zone of cell
-            void modifyCell(const label, const label zoneID);
+            void modifyCell(const label celli, const label zoneID);
 
             //- Remove/merge cell.
-            void removeCell(const label, const label);
+            void removeCell(const label celli, const label mergeCelli);
 
             //- Explicitly set the number of patches if construct-without-mesh
             //  used.
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChangeTemplates.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChangeTemplates.C
index ddca640d04..731d0beb01 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChangeTemplates.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChangeTemplates.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  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,20 +30,20 @@ License
 template<class T>
 void Foam::polyTopoChange::reorder
 (
-    const labelList& oldToNew,
+    const labelUList& oldToNew,
     DynamicList<T>& lst
 )
 {
     // Create copy
     DynamicList<T> oldLst(lst);
 
-    forAll(oldToNew, elemI)
+    forAll(oldToNew, i)
     {
-        label newElemI = oldToNew[elemI];
+        const label newIdx = oldToNew[i];
 
-        if (newElemI != -1)
+        if (newIdx >= 0)
         {
-            lst[newElemI] = oldLst[elemI];
+            lst[newIdx] = oldLst[i];
         }
     }
 }
@@ -52,20 +52,20 @@ void Foam::polyTopoChange::reorder
 template<class T>
 void Foam::polyTopoChange::reorder
 (
-    const labelList& oldToNew,
+    const labelUList& oldToNew,
     List<DynamicList<T>>& lst
 )
 {
     // Create copy
     List<DynamicList<T>> oldLst(lst);
 
-    forAll(oldToNew, elemI)
+    forAll(oldToNew, i)
     {
-        label newElemI = oldToNew[elemI];
+        const label newIdx = oldToNew[i];
 
-        if (newElemI != -1)
+        if (newIdx >= 0)
         {
-            lst[newElemI].transfer(oldLst[elemI]);
+            lst[newIdx].transfer(oldLst[i]);
         }
     }
 }
@@ -74,23 +74,23 @@ void Foam::polyTopoChange::reorder
 template<class T>
 void Foam::polyTopoChange::renumberKey
 (
-    const labelList& oldToNew,
-    Map<T>& elems
+    const labelUList& oldToNew,
+    Map<T>& map
 )
 {
-    Map<T> newElems(elems.size());
+    Map<T> newMap(map.capacity());
 
-    forAllConstIter(typename Map<T>, elems, iter)
+    forAllConstIters(map, iter)
     {
-        label newElem = oldToNew[iter.key()];
+        const label newKey = oldToNew[iter.key()];
 
-        if (newElem >= 0)
+        if (newKey >= 0)
         {
-            newElems.insert(newElem, iter());
+            newMap.insert(newKey, iter.object());
         }
     }
 
-    elems.transfer(newElems);
+    map.transfer(newMap);
 }
 
 
diff --git a/src/dynamicMesh/slidingInterface/coupleSlidingInterface.C b/src/dynamicMesh/slidingInterface/coupleSlidingInterface.C
index 1d5794a507..9257fa57d2 100644
--- a/src/dynamicMesh/slidingInterface/coupleSlidingInterface.C
+++ b/src/dynamicMesh/slidingInterface/coupleSlidingInterface.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  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -44,6 +44,7 @@ License
 
 const Foam::scalar Foam::slidingInterface::edgeCoPlanarTolDefault_ = 0.8;
 
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 // Index of debug signs:
@@ -69,16 +70,13 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
 {
     if (debug)
     {
-        Pout<< "void slidingInterface::coupleInterface"
-            << "(polyTopoChange& ref) : "
-            << "Coupling sliding interface " << name() << endl;
+        Pout<< FUNCTION_NAME << nl
+            << ": Coupling sliding interface " << name() << endl;
     }
 
     const polyMesh& mesh = topoChanger().mesh();
-
     const pointField& points = mesh.points();
     const faceList& faces = mesh.faces();
-
     const labelList& own = mesh.faceOwner();
     const labelList& nei = mesh.faceNeighbour();
     const faceZoneMesh& faceZones = mesh.faceZones();
@@ -156,20 +154,19 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
     // Create mapping for every merged point of the slave patch
     forAll(slavePointPointHits, pointi)
     {
-        if (slavePointPointHits[pointi] >= 0)
-        {
-            // Pout<< "Inserting point merge pair: " << slaveMeshPoints[pointi]
-            //     << " : " << masterMeshPoints[slavePointPointHits[pointi]]
-            //     << endl;
+        const label slaveHitPti = slavePointPointHits[pointi];
 
+        if (slaveHitPti >= 0)
+        {
             pointMergeMap.insert
             (
                 slaveMeshPoints[pointi],
-                masterMeshPoints[slavePointPointHits[pointi]]
+                masterMeshPoints[slaveHitPti]
             );
         }
     }
 
+
     // Collect the list of used edges for every slave edge
 
     List<labelHashSet> usedMasterEdges(slaveEdges.size());
@@ -187,22 +184,20 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
 
             // Mark all current master edges as used for all the current slave
             // edges
-            forAll(curSlaveEdges, slaveEdgeI)
+            for (const label slaveEdgei : curSlaveEdges)
             {
-                labelHashSet& sm = usedMasterEdges[curSlaveEdges[slaveEdgeI]];
-
-                forAll(curMasterEdges, masterEdgeI)
-                {
-                    sm.insert(curMasterEdges[masterEdgeI]);
-                }
+                usedMasterEdges[slaveEdgei].insert
+                (
+                    curMasterEdges
+                );
             }
         }
         else if (slavePointEdgeHits[pointi] > -1)
         {
             // For edge hits, add the master edge
-            forAll(curSlaveEdges, slaveEdgeI)
+            for (const label slaveEdgei : curSlaveEdges)
             {
-                usedMasterEdges[curSlaveEdges[slaveEdgeI]].insert
+                usedMasterEdges[slaveEdgei].insert
                 (
                     slavePointEdgeHits[pointi]
                 );
@@ -219,21 +214,18 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
         {
             const labelList& curMasterEdges = masterPointEdges[masterPointi];
 
-            labelHashSet& sm =
-                usedMasterEdges[masterPointEdgeHits[masterPointi]];
-
-            forAll(curMasterEdges, masterEdgeI)
-            {
-                sm.insert(curMasterEdges[masterEdgeI]);
-            }
+            usedMasterEdges[masterPointEdgeHits[masterPointi]].insert
+            (
+                curMasterEdges
+            );
         }
     }
 
     // Pout<< "used edges: " << endl;
-    // forAll(usedMasterEdges, edgeI)
+    // forAll(usedMasterEdges, edgei)
     // {
-    //     Pout<< "edge: " << edgeI
-    //         << " used: " << usedMasterEdges[edgeI].toc()
+    //     Pout<< "edge: " << edgei
+    //         << " used: " << usedMasterEdges[edgei].toc()
     //         << endl;
     // }
 
@@ -249,13 +241,13 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
         {
             // Create a new point on the master edge
 
-            point edgeCutPoint =
+            const point edgeCutPoint =
                 masterEdges[slavePointEdgeHits[pointi]].line
                 (
                     masterLocalPoints
                 ).nearestDist(projectedSlavePoints[pointi]).hitPoint();
 
-            label newPoint =
+            const label newPointi =
                 ref.setAction
                 (
                     polyAddPoint
@@ -267,32 +259,37 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
                     )
                 );
 
-            // Pout<< "Inserting merge pair off edge: "
-            //     << slaveMeshPoints[pointi] << " " << newPoint
-            //     << " cut point: " << edgeCutPoint
-            //     << " orig: " << slaveLocalPoints[pointi]
-            //     << " proj: " << projectedSlavePoints[pointi]
-            //     << endl;
-
             // Add the new edge point into the merge map
-            pointMergeMap.insert(slaveMeshPoints[pointi], newPoint);
+            pointMergeMap.insert
+            (
+                slaveMeshPoints[pointi],
+                newPointi
+            );
 
             pointsIntoMasterEdges[slavePointEdgeHits[pointi]].append
             (
-                newPoint
+                newPointi
             );
 
             // Add the point into the enriched patch map
             pointMap.insert
             (
-                newPoint,
+                newPointi,
                 edgeCutPoint
             );
 
             if (debug)
             {
                 Pout<< "e";
-                // Pout<< newPoint << " = " << edgeCutPoint << endl;
+
+                // Pout<< "Inserting merge pair off edge: "
+                //     << slaveMeshPoints[pointi] << " " << newPointi
+                //     << " cut point: " << edgeCutPoint
+                //     << " orig: " << slaveLocalPoints[pointi]
+                //     << " proj: " << projectedSlavePoints[pointi]
+                //     << endl;
+
+                // Pout<< newPointi << " = " << edgeCutPoint << endl;
             }
         }
     }
@@ -312,7 +309,7 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
          && slavePointFaceHits[pointi].hit()
         )
         {
-            label newPoint =
+            const label newPointi =
                 ref.setAction
                 (
                     polyAddPoint
@@ -326,22 +323,26 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
 
             // Pout<< "Inserting merge pair off face: "
             //     << slaveMeshPoints[pointi]
-            //     << " " << newPoint
+            //     << " " << newPointi
             //     << endl;
 
             // Add the new edge point into the merge map
-            pointMergeMap.insert(slaveMeshPoints[pointi], newPoint);
+            pointMergeMap.insert
+            (
+                slaveMeshPoints[pointi],
+                newPointi
+            );
 
             // Add the point into the enriched patch map
             pointMap.insert
             (
-                newPoint,
+                newPointi,
                 projectedSlavePoints[pointi]
             );
 
             if (debug)
             {
-                Pout<< "f: " << newPoint << " = "
+                Pout<< "f: " << newPointi << " = "
                     << projectedSlavePoints[pointi] << endl;
             }
         }
@@ -390,7 +391,7 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
     //
     if (debug)
     {
-        Pout<< "Processing slave edges " << endl;
+        Pout<< "Processing slave edges" << endl;
     }
 
     if (!cutPointEdgePairMapPtr_)
@@ -430,65 +431,72 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
             //     << " curUme: " << curUme
             //     << endl;
 
-            // Clear the maps
-            curFaceMap.clear();
-            addedFaces.clear();
-
             // Grab the faces for start and end points.
             const label startFace =
                 slavePointFaceHits[curEdge.start()].hitObject();
-            const label endFace = slavePointFaceHits[curEdge.end()].hitObject();
+
+            const label endFace =
+                slavePointFaceHits[curEdge.end()].hitObject();
 
             // Pout<< "startFace: " << slavePointFaceHits[curEdge.start()]
             //     << " endFace: " << slavePointFaceHits[curEdge.end()]
             //     << endl;
 
-            // Insert the start face into the list
-            curFaceMap.insert(startFace);
-            addedFaces.insert(startFace);
-
-            // Pout<< "curFaceMap: " << curFaceMap.toc() << endl;
-
-            label nSweeps = 0;
             bool completed = false;
 
-            while (nSweeps < edgeFaceEscapeLimit_)
+            if (!completed)
             {
-                nSweeps++;
+                // Forward sweep
 
-                if (addedFaces.found(endFace))
-                {
-                    completed = true;
-                }
-
-                // Add all face neighbours of face in the map
-                const labelList cf = addedFaces.toc();
+                // Clear the maps
+                curFaceMap.clear();
                 addedFaces.clear();
 
-                forAll(cf, cfI)
-                {
-                    const labelList& curNbrs = masterFaceFaces[cf[cfI]];
+                // Insert the start face into the list
+                curFaceMap.insert(startFace);
+                addedFaces.insert(startFace);
 
-                    forAll(curNbrs, nbrI)
+                // Pout<< "curFaceMap: " << curFaceMap.toc() << endl;
+
+                for
+                (
+                    label nSweeps = 0;
+                    nSweeps < edgeFaceEscapeLimit_;
+                    ++nSweeps
+                )
+                {
+                    completed = addedFaces.found(endFace);
+
+                    // Add all face neighbours of face in the map
+                    const labelList cf(std::move(addedFaces.toc()));
+                    addedFaces.clear();
+
+                    for (const label cfi : cf)
                     {
-                        if (!curFaceMap.found(curNbrs[nbrI]))
+                        const labelList& curNbrs = masterFaceFaces[cfi];
+
+                        for (const label nbri : curNbrs)
                         {
-                            curFaceMap.insert(curNbrs[nbrI]);
-                            addedFaces.insert(curNbrs[nbrI]);
+                            if (curFaceMap.insert(nbri))
+                            {
+                                addedFaces.insert(nbri);
+                            }
                         }
                     }
-                }
 
-                if (completed) break;
+                    if (completed) break;
 
-                if (debug)
-                {
-                    Pout<< ".";
+                    if (debug)
+                    {
+                        Pout<< ".";
+                    }
                 }
             }
 
             if (!completed)
             {
+                // Reverse sweep
+
                 if (debug)
                 {
                     Pout<< "x";
@@ -497,34 +505,31 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
                 // It is impossible to reach the end from the start, probably
                 // due to disconnected domain.  Do search in opposite direction
 
-                label nReverseSweeps = 0;
-
                 addedFaces.clear();
                 addedFaces.insert(endFace);
 
-                while (nReverseSweeps < edgeFaceEscapeLimit_)
+                for
+                (
+                    label nSweeps = 0;
+                    nSweeps < edgeFaceEscapeLimit_;
+                    ++nSweeps
+                )
                 {
-                    nReverseSweeps++;
-
-                    if (addedFaces.found(startFace))
-                    {
-                        completed = true;
-                    }
+                    completed = addedFaces.found(startFace);
 
                     // Add all face neighbours of face in the map
-                    const labelList cf = addedFaces.toc();
+                    const labelList cf(std::move(addedFaces.toc()));
                     addedFaces.clear();
 
-                    forAll(cf, cfI)
+                    for (const label cfi : cf)
                     {
-                        const labelList& curNbrs = masterFaceFaces[cf[cfI]];
+                        const labelList& curNbrs = masterFaceFaces[cfi];
 
-                        forAll(curNbrs, nbrI)
+                        for (const label nbri : curNbrs)
                         {
-                            if (!curFaceMap.found(curNbrs[nbrI]))
+                            if (curFaceMap.insert(nbri))
                             {
-                                curFaceMap.insert(curNbrs[nbrI]);
-                                addedFaces.insert(curNbrs[nbrI]);
+                                addedFaces.insert(nbri);
                             }
                         }
                     }
@@ -561,27 +566,22 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
                 nFacesPerSlaveEdge_*primitiveMesh::edgesPerFace_
             );
 
-            const labelList curFaces = curFaceMap.toc();
-
-            // Pout<< "curFaces: " << curFaces << endl;
-
-            forAll(curFaces, facei)
+            for (const label facei : curFaceMap)
             {
-                // Pout<< "face: " << curFaces[facei] << " "
-                //     << masterPatch[curFaces[facei]]
+                // Pout<< "face: " << facei << " "
+                //     << masterPatch[facei]
                 //     << " local: "
-                //     << masterPatch.localFaces()[curFaces[facei]]
+                //     << masterPatch.localFaces()[facei]
                 //     << endl;
 
-                const labelList& me = masterFaceEdges[curFaces[facei]];
-
-                forAll(me, meI)
-                {
-                    curMasterEdgesMap.insert(me[meI]);
-                }
+                curMasterEdgesMap.insert
+                (
+                    masterFaceEdges[facei]
+                );
             }
 
-            const labelList curMasterEdges = curMasterEdgesMap.toc();
+            const labelList curMasterEdges(std::move(curMasterEdgesMap.toc()));
+            curMasterEdgesMap.clear();
 
             // For all master edges to intersect, skip the ones
             // already used and cut the rest with a cutting plane.  If
@@ -597,7 +597,7 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
             const point& a = projectedSlavePoints[curEdge.start()];
             const point& b = projectedSlavePoints[curEdge.end()];
 
-            point c =
+            const point c =
                 0.5*
                 (
                     slaveLocalPoints[curEdge.start()]
@@ -607,7 +607,7 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
                 );
 
             // Create the plane
-            plane cutPlane(a, b, c);
+            const plane cutPlane(a, b, c);
 
             // Pout<< "a: " << a
             //     << " b: " << b
@@ -615,14 +615,16 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
             //     << " plane: " << cutPlane
             //     << endl;
 
-            linePointRef curSlaveLine = curEdge.line(projectedSlavePoints);
+            const linePointRef curSlaveLine =
+                curEdge.line(projectedSlavePoints);
+
             const scalar curSlaveLineMag = curSlaveLine.mag();
 
             // Pout<< "curSlaveLine: " << curSlaveLine << endl;
 
-            forAll(curMasterEdges, masterEdgeI)
+            for (const label cmeIndex : curMasterEdges)
             {
-                if (!curUme.found(curMasterEdges[masterEdgeI]))
+                if (!curUme.found(cmeIndex))
                 {
                     // New edge
                     if (debug)
@@ -630,7 +632,6 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
                         Pout<< "n";
                     }
 
-                    const label cmeIndex = curMasterEdges[masterEdgeI];
                     const edge& cme = masterEdges[cmeIndex];
 
                     // Pout<< "Edge " << cmeIndex
@@ -638,7 +639,7 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
                     //     << " line: " << cme.line(masterLocalPoints)
                     //     << endl;
 
-                    scalar cutOnMaster =
+                    const scalar cutOnMaster =
                         cutPlane.lineIntersect
                         (
                             cme.line(masterLocalPoints)
@@ -651,18 +652,18 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
                     )
                     {
                         // Master is cut, check the slave
-                        point masterCutPoint =
+                        const point masterCutPoint =
                             masterLocalPoints[cme.start()]
                           + cutOnMaster*cme.vec(masterLocalPoints);
 
-                        pointHit slaveCut =
+                        const pointHit slaveCut =
                             curSlaveLine.nearestDist(masterCutPoint);
 
                         if (slaveCut.hit())
                         {
                             // Strict checking of slave cut to avoid capturing
                             // end points.
-                            scalar cutOnSlave =
+                            const scalar cutOnSlave =
                                 (
                                     (
                                         slaveCut.hitPoint()
@@ -672,7 +673,7 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
 
                             // Calculate merge tolerance from the
                             // target edge length
-                            scalar mergeTol = edgeCoPlanarTol_*mag(b - a);
+                            const scalar mergeTol = edgeCoPlanarTol_*mag(b - a);
 
                             // Pout<< "cutOnMaster: " << cutOnMaster
                             //     << " masterCutPoint: " << masterCutPoint
@@ -697,7 +698,7 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
                                 // to edge points The point is nominally
                                 // added from the start of the master edge
                                 // and added to the cut point zone
-                                label newPoint =
+                                const label newPointi =
                                     ref.setAction
                                     (
                                         polyAddPoint
@@ -709,7 +710,7 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
                                         )
                                     );
 
-                                // Pout<< "Inserting point: " << newPoint
+                                // Pout<< "Inserting point: " << newPointi
                                 //     << " as edge to edge intersection.  "
                                 //     << "Slave edge: "
                                 //     << edgeI << " " << curEdge
@@ -717,16 +718,19 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
                                 //     << cmeIndex << " " << cme
                                 //     << endl;
 
-                                pointsIntoSlaveEdges[edgeI].append(newPoint);
+                                pointsIntoSlaveEdges[edgeI].append
+                                (
+                                    newPointi
+                                );
                                 pointsIntoMasterEdges[cmeIndex].append
                                 (
-                                    newPoint
+                                    newPointi
                                 );
 
                                 // Add the point into the enriched patch map
                                 pointMap.insert
                                 (
-                                    newPoint,
+                                    newPointi,
                                     masterCutPoint
                                 );
 
@@ -734,25 +738,25 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
                                 // create cut point
                                 addToCpepm.insert
                                 (
-                                    newPoint,    // Cut point index
+                                    newPointi,   // Cut point index
                                     Pair<edge>
                                     (
                                         edge
                                         (
                                             masterMeshPoints[cme.start()],
                                             masterMeshPoints[cme.end()]
-                                        ),    // Master edge
+                                        ),  // Master edge
                                         edge
                                         (
                                             slaveMeshPoints[curEdge.start()],
                                             slaveMeshPoints[curEdge.end()]
-                                        )// Slave edge
+                                        )   // Slave edge
                                     )
                                 );
 
                                 if (debug)
                                 {
-                                    Pout<< " " << newPoint << " = "
+                                    Pout<< " " << newPointi << " = "
                                         << masterCutPoint << " ";
                                 }
                             }
@@ -955,12 +959,7 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
 
             forAll(rsf, i)
             {
-                Map<label>::const_iterator mpIter = pointMergeMap.find(rsf[i]);
-
-                if (mpIter != pointMergeMap.end())
-                {
-                    rsf[i] = mpIter();
-                }
+                rsf[i] = pointMergeMap.lookup(rsf[i], rsf[i]);
             }
 
             if (curCutFace == rsf)
@@ -1240,10 +1239,10 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
 
     if (debug)
     {
-        Pout<< "Number of orphaned faces: "
-            << "master = " << nOrphanedMasters << " out of "
+        Pout<< "Orphaned faces: "
+            << "master = " << nOrphanedMasters << "/"
             << orphanedMaster.size()
-            << " slave = " << nOrphanedSlaves << " out of "
+            << " slave = " << nOrphanedSlaves << "/"
             << orphanedSlave.size() << endl;
     }
 
@@ -1268,12 +1267,9 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
     // Pout<< "masterStickOuts: " << masterStickOuts << endl;
 
     // Re-create the master stick-out faces
-    forAll(masterStickOuts, facei)
+    for (const label curFaceID : masterStickOuts)
     {
         // Renumber the face and remove additional points
-
-        const label curFaceID = masterStickOuts[facei];
-
         const face& oldRichFace = faces[curFaceID];
 
         bool changed = false;
@@ -1282,47 +1278,46 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
         face oldFace(oldRichFace.size());
         label nOldFace = 0;
 
-        forAll(oldRichFace, pointi)
+        for (const label pointi : oldRichFace)
         {
-            if (ref.pointRemoved(oldRichFace[pointi]))
+            if (ref.pointRemoved(pointi))
             {
                 changed = true;
             }
             else
             {
                 // Point off patch
-                oldFace[nOldFace] = oldRichFace[pointi];
+                oldFace[nOldFace] = pointi;
                 nOldFace++;
             }
         }
 
         oldFace.setSize(nOldFace);
 
-        // Pout<< "old rich master face: " << oldRichFace
-        //     << " old face: " << oldFace
-        //     << endl;
+        // Pout<< "old rich face[" << curFaceID << "]: " << oldRichFace
+        //     << " old face: " << oldFace << endl;
 
         DynamicList<label> newFaceLabels(2*oldFace.size());
 
         forAll(oldFace, pointi)
         {
-            if (masterMeshPointMap.found(oldFace[pointi]))
+            const label localFirstLabel =
+                masterMeshPointMap.lookup(oldFace[pointi], -1);
+
+            if (localFirstLabel != -1)
             {
                 // Point is in master patch. Add it
 
-                // If the point is a direct hit, grab its label; otherwise
-                // keep the original
-                if (pointMergeMap.found(oldFace[pointi]))
+                // If the point is a direct hit, grab its label;
+                // otherwise keep the original
+                newFaceLabels.append
+                (
+                    pointMergeMap.lookup(oldFace[pointi], oldFace[pointi])
+                );
+
+                if (newFaceLabels.last() != oldFace[pointi])
                 {
                     changed = true;
-                    newFaceLabels.append
-                    (
-                        pointMergeMap.find(oldFace[pointi])()
-                    );
-                }
-                else
-                {
-                    newFaceLabels.append(oldFace[pointi]);
                 }
 
                 // Find if there are additional points inserted onto the edge
@@ -1331,39 +1326,35 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
                 // 1) Find all the edges in the master patch coming
                 //    out of the current point.
                 // 2) If the next point in the face to pick the right edge
-                const label localFirstLabel =
-                    masterMeshPointMap.find(oldFace[pointi])();
 
                 const labelList& curEdges = masterPointEdges[localFirstLabel];
 
-                const label  nextLabel = oldFace.nextLabel(pointi);
+                const label nextLabel = oldFace.nextLabel(pointi);
 
-                Map<label>::const_iterator mmpmIter =
-                    masterMeshPointMap.find(nextLabel);
+                const label localNextLabel =
+                    masterMeshPointMap.lookup(nextLabel, -1);
 
-                if (mmpmIter != masterMeshPointMap.end())
+                if (localNextLabel != -1)
                 {
                     // Pout<< "found label pair " << oldFace[pointi]
                     //     << " and " << nextLabel;
                     // Find the points on the edge between them
-                    const label localNextLabel = mmpmIter();
 
-                    forAll(curEdges, curEdgeI)
+                    for (const label curEdgei : curEdges)
                     {
                         if
                         (
-                            masterEdges[curEdges[curEdgeI]].otherVertex
+                            masterEdges[curEdgei].otherVertex
                             (
                                 localFirstLabel
                             )
                          == localNextLabel
                         )
                         {
-                            // Pout<< " found edge: " << curEdges[curEdgeI]
-                            //     << endl;
+                            // Pout<< " found edge: " << curEdgei << endl;
 
                             // Get points on current edge
-                            const labelList& curPime = pime[curEdges[curEdgeI]];
+                            const labelList& curPime = pime[curEdgei];
 
                             if (curPime.size())
                             {
@@ -1374,9 +1365,10 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
                                 const point& startPoint =
                                     masterLocalPoints[localFirstLabel];
 
-                                vector e =
-                                    masterLocalPoints[localNextLabel]
-                                  - startPoint;
+                                const point& endPoint =
+                                    masterLocalPoints[localNextLabel];
+
+                                vector e = (endPoint - startPoint);
 
                                 e /= magSqr(e);
 
@@ -1385,13 +1377,13 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
                                 forAll(curPime, curPimeI)
                                 {
                                     edgePointWeights[curPimeI] =
-                                        (
-                                            e
-                                          & (
-                                              pointMap.find(curPime[curPimeI])()
-                                            - startPoint
-                                            )
-                                        );
+                                    (
+                                        e
+                                      & (
+                                           pointMap[curPime[curPimeI]]
+                                         - startPoint
+                                        )
+                                    );
                                 }
 
                                 if (debug)
@@ -1466,25 +1458,24 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
             }
 
             // Get face zone and its flip
-            label modifiedFaceZone = faceZones.whichZone(curFaceID);
-            bool modifiedFaceZoneFlip = false;
+            const label modifiedFaceZone = faceZones.whichZone(curFaceID);
 
-            if (modifiedFaceZone >= 0)
-            {
-                modifiedFaceZoneFlip =
-                    faceZones[modifiedFaceZone].flipMap()
-                    [
-                        faceZones[modifiedFaceZone].whichFace(curFaceID)
-                    ];
-            }
+            const bool modifiedFaceZoneFlip =
+            (
+                modifiedFaceZone >= 0
+              ?
+                faceZones[modifiedFaceZone].flipMap()
+                [
+                    faceZones[modifiedFaceZone].whichFace(curFaceID)
+                ]
+              : false
+            );
 
             face newFace;
             newFace.transfer(newFaceLabels);
 
-            // Pout<< "Modifying master stick-out face " << curFaceID
-            //     << " old face: " << oldFace
-            //     << " new face: " << newFace
-            //     << endl;
+            // Pout<< "Modifying master stick-out face[" << curFaceID
+            //     << "]: old: " << oldFace << " new: " << newFace << endl;
 
             // Modify the face
             if (mesh.isInternalFace(curFaceID))
@@ -1540,11 +1531,9 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
 
     // Re-create the slave stick-out faces
 
-    forAll(slaveStickOuts, facei)
+    for (const label curFaceID : slaveStickOuts)
     {
         // Renumber the face and remove additional points
-        const label curFaceID = slaveStickOuts[facei];
-
         const face& oldRichFace = faces[curFaceID];
 
         bool changed = false;
@@ -1553,33 +1542,28 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
         face oldFace(oldRichFace.size());
         label nOldFace = 0;
 
-        forAll(oldRichFace, pointi)
+        for (const label pointi : oldRichFace)
         {
             if
             (
-                rpm.found(oldRichFace[pointi])
-             || slaveMeshPointMap.found(oldRichFace[pointi])
+                rpm.found(pointi)
+             || slaveMeshPointMap.found(pointi)
             )
             {
                 // Point definitely live. Add it
-                oldFace[nOldFace] = oldRichFace[pointi];
+                oldFace[nOldFace] = pointi;
                 nOldFace++;
             }
-            else if
-            (
-                ref.pointRemoved(oldRichFace[pointi])
-             || masterMeshPointMap.found(oldRichFace[pointi])
-            )
+            else if (ref.pointRemoved(pointi))
             {
-                // Point removed and not on slave patch
-                // (first if takes care of that!) or
-                // point belonging to master patch
+                // Point removed, not on slave patch and not retired
+                // (first if takes care of that!)
                 changed = true;
             }
             else
             {
-                // Point off patch
-                oldFace[nOldFace] = oldRichFace[pointi];
+                // Point on master or slave patch
+                oldFace[nOldFace] = pointi;
                 nOldFace++;
             }
         }
@@ -1595,33 +1579,29 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
         forAll(oldFace, pointi)
         {
             // Try to find the point in retired points
-            label curP = oldFace[pointi];
+            const label curP = rpm.lookup(oldFace[pointi], oldFace[pointi]);
 
-            Map<label>::const_iterator rpmIter = rpm.find(oldFace[pointi]);
-
-            if (rpmIter != rpm.end())
+            if (curP != oldFace[pointi])
             {
                 changed = true;
-                curP = rpmIter();
             }
 
-            if (slaveMeshPointMap.found(curP))
+            const label localFirstLabel = slaveMeshPointMap.lookup(curP, -1);
+
+            if (localFirstLabel != -1)
             {
                 // Point is in slave patch. Add it
 
-                // If the point is a direct hit, grab its label; otherwise
-                // keep the original
-                if (pointMergeMap.found(curP))
+                // If the point is a direct hit, grab its label;
+                // otherwise keep the original
+                newFaceLabels.append
+                (
+                    pointMergeMap.lookup(curP, curP)
+                );
+
+                if (newFaceLabels.last() != curP)
                 {
                     changed = true;
-                    newFaceLabels.append
-                    (
-                        pointMergeMap.find(curP)()
-                    );
-                }
-                else
-                {
-                    newFaceLabels.append(curP);
                 }
 
                 // Find if there are additional points inserted onto the edge
@@ -1631,59 +1611,51 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
                 //    out of the current point.
                 // 2) Use the next point in the face to pick the right edge
 
-                const label localFirstLabel =
-                    slaveMeshPointMap.find(curP)();
-
                 const labelList& curEdges = slavePointEdges[localFirstLabel];
 
                 label nextLabel = oldFace.nextLabel(pointi);
 
-                Map<label>::const_iterator rpmNextIter =
-                    rpm.find(nextLabel);
+                nextLabel = rpm.lookup(nextLabel, nextLabel);
 
-                if (rpmNextIter != rpm.end())
-                {
-                    nextLabel = rpmNextIter();
-                }
+                const label localNextLabel =
+                    slaveMeshPointMap.lookup(nextLabel, -1);
 
-                Map<label>::const_iterator mmpmIter =
-                    slaveMeshPointMap.find(nextLabel);
-
-                if (mmpmIter != slaveMeshPointMap.end())
+                if (localNextLabel != -1)
                 {
                     // Both points on the slave patch.
                     // Find the points on the edge between them
-                    const label localNextLabel = mmpmIter();
 
-                    forAll(curEdges, curEdgeI)
+                    for (const label curEdgei : curEdges)
                     {
                         if
                         (
-                            slaveEdges[curEdges[curEdgeI]].otherVertex
+                            slaveEdges[curEdgei].otherVertex
                             (
                                 localFirstLabel
                             )
                          == localNextLabel
                         )
                         {
-                            // Pout<< " found edge: " << curEdges[curEdgeI]
+                            // Pout<< " found edge: " << curEdgei
                             //     << endl;
 
                             // Get points on current edge
-                            const labelList& curPise = pise[curEdges[curEdgeI]];
+                            const labelList& curPise = pise[curEdgei];
 
                             if (curPise.size())
                             {
                                 changed = true;
+
                                 // Pout<< "curPise: " << curPise << endl;
                                 // Insert the edge points into the face
                                 // in the correct order
                                 const point& startPoint =
                                     projectedSlavePoints[localFirstLabel];
 
-                                vector e =
-                                    projectedSlavePoints[localNextLabel]
-                                  - startPoint;
+                                const point& endPoint =
+                                    projectedSlavePoints[localNextLabel];
+
+                                vector e = endPoint - startPoint;
 
                                 e /= magSqr(e);
 
@@ -1695,7 +1667,7 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
                                     (
                                         e
                                       & (
-                                            pointMap.find(curPise[curPiseI])()
+                                            pointMap[curPise[curPiseI]]
                                           - startPoint
                                         )
                                     );
@@ -1773,24 +1745,25 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
             }
 
             // Get face zone and its flip
-            label modifiedFaceZone = faceZones.whichZone(curFaceID);
-            bool modifiedFaceZoneFlip = false;
+            const label modifiedFaceZone =
+                faceZones.whichZone(curFaceID);
 
-            if (modifiedFaceZone >= 0)
-            {
-                modifiedFaceZoneFlip =
-                    faceZones[modifiedFaceZone].flipMap()
-                    [
-                        faceZones[modifiedFaceZone].whichFace(curFaceID)
-                    ];
-            }
+            const bool modifiedFaceZoneFlip
+            (
+                modifiedFaceZone >= 0
+              ?
+                faceZones[modifiedFaceZone].flipMap()
+                [
+                    faceZones[modifiedFaceZone].whichFace(curFaceID)
+                ]
+              : false
+            );
 
             face newFace;
             newFace.transfer(newFaceLabels);
 
-            // Pout<< "Modifying slave stick-out face " << curFaceID
-            //     << " old face: " << oldFace
-            //     << " new face: " << newFace
+            // Pout<< "Modifying slave stick-out face[" << curFaceID
+            //     << "]: old: " << oldFace << " new: " << newFace
             //     << endl;
 
             // Modify the face
@@ -1852,9 +1825,16 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
 
     label nRetiredPoints = 0;
 
-    forAll(slaveMeshPoints, pointi)
+    for (const label slavePointi : slaveMeshPoints)
     {
-        if (pointMergeMap.found(slaveMeshPoints[pointi]))
+        const label masterPointi = pointMergeMap.lookup(slavePointi, -1);
+
+        if (slavePointi == masterPointi)
+        {
+            // Identity mapping (ie, slave point already exists on master patch)
+            continue;
+        }
+        else if (masterPointi != -1)
         {
             // Retire the point - only used for supporting the detached
             // slave patch
@@ -1862,32 +1842,28 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
 
             // ref.setAction
             // (
-            //    polyModifyPoint
-            //    (
-            //        slaveMeshPoints[pointi],             // point ID
-            //        points[slaveMeshPoints[pointi]],     // point
-            //        false,                               // remove from zone
-            //        mesh.pointZones().whichZone(slaveMeshPoints[pointi]),
-            //                                             // zone
-            //        false                                // in a cell
-            //    )
+            //     polyModifyPoint
+            //     (
+            //         slavePointi,                         // point ID
+            //         points[slavePointi],                 // point
+            //         false,                               // remove from zone
+            //         mesh.pointZones().whichZone(slavePointi), // zone
+            //         false                                // in a cell
+            //     )
             // );
-            //Pout<< "MJ retire slave point " << slaveMeshPoints[pointi]
-            //    << " coord " << points[slaveMeshPoints[pointi]]
+            //
+            //Pout<< "MJ retire slave point " << slavePointi
+            //    << " coord " << points[slavePointi]
             //    << endl;
             ref.setAction
             (
                 polyRemovePoint
                 (
-                    slaveMeshPoints[pointi]
+                    slavePointi
                 )
             );
 
-            addToRpm.insert
-            (
-                pointMergeMap.find(slaveMeshPoints[pointi])(),
-                slaveMeshPoints[pointi]
-            );
+            addToRpm.insert(masterPointi, slavePointi);
         }
         else
         {
@@ -1895,10 +1871,10 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
             (
                 polyModifyPoint
                 (
-                    slaveMeshPoints[pointi],             // point ID
-                    points[slaveMeshPoints[pointi]],     // point
+                    slavePointi,                         // point ID
+                    points[slavePointi],                 // point
                     false,                               // remove from zone
-                    mesh.pointZones().whichZone(slaveMeshPoints[pointi]),// zone
+                    mesh.pointZones().whichZone(slavePointi),// zone
                     true                                 // in a cell
                 )
             );
@@ -1907,15 +1883,15 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
 
     if (debug)
     {
-        Pout<< "Retired " << nRetiredPoints << " out of "
+        Pout<< "Retired " << nRetiredPoints << "/"
             << slaveMeshPoints.size() << " points." << endl;
     }
 
     // Grab cut face master and slave addressing
-    if (cutFaceMasterPtr_) deleteDemandDrivenData(cutFaceMasterPtr_);
+    deleteDemandDrivenData(cutFaceMasterPtr_);
     cutFaceMasterPtr_ = new labelList(cutPatch.cutFaceMaster());
 
-    if (cutFaceSlavePtr_) deleteDemandDrivenData(cutFaceSlavePtr_);
+    deleteDemandDrivenData(cutFaceSlavePtr_);
     cutFaceSlavePtr_ = new labelList(cutPatch.cutFaceSlave());
 
     // Finished coupling
@@ -1923,9 +1899,8 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
 
     if (debug)
     {
-        Pout<< "void slidingInterface::coupleInterface("
-            << "polyTopoChange& ref) : "
-            << "Finished coupling sliding interface " << name() << endl;
+        Pout<< FUNCTION_NAME << nl
+            << ": Finished coupling sliding interface " << name() << endl;
     }
 }
 
diff --git a/src/dynamicMesh/slidingInterface/decoupleSlidingInterface.C b/src/dynamicMesh/slidingInterface/decoupleSlidingInterface.C
index e454e31981..95c3f91f2e 100644
--- a/src/dynamicMesh/slidingInterface/decoupleSlidingInterface.C
+++ b/src/dynamicMesh/slidingInterface/decoupleSlidingInterface.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  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,18 +40,16 @@ void Foam::slidingInterface::decoupleInterface
 {
     if (debug)
     {
-        Pout<< "void slidingInterface::decoupleInterface("
-            << "polyTopoChange& ref) const : "
-            << "Decoupling sliding interface " << name() << endl;
+        Pout<< FUNCTION_NAME << nl
+            << ": Decoupling sliding interface " << name() << endl;
     }
 
     if (!attached_)
     {
         if (debug)
         {
-            Pout<< "void slidingInterface::decoupleInterface("
-                << "polyTopoChange& ref) const : "
-                << "Interface already decoupled." << endl;
+            Pout<< FUNCTION_NAME << nl
+                << ": Interface already decoupled." << endl;
         }
 
         return;
@@ -61,22 +59,23 @@ void Foam::slidingInterface::decoupleInterface
     clearCouple(ref);
 
     const polyMesh& mesh = topoChanger().mesh();
+    const pointField& points = mesh.points();
     const faceList& faces = mesh.faces();
     const cellList& cells = mesh.cells();
-
     const labelList& own = mesh.faceOwner();
     const labelList& nei = mesh.faceNeighbour();
+    const faceZoneMesh& faceZones = mesh.faceZones();
 
     // Master side
 
     const primitiveFacePatch& masterPatch =
-        mesh.faceZones()[masterFaceZoneID_.index()]();
+        faceZones[masterFaceZoneID_.index()]();
 
     const labelList& masterPatchAddr =
-        mesh.faceZones()[masterFaceZoneID_.index()];
+        faceZones[masterFaceZoneID_.index()];
 
     const boolList& masterPatchFlip =
-        mesh.faceZones()[masterFaceZoneID_.index()].flipMap();
+        faceZones[masterFaceZoneID_.index()].flipMap();
 
     const labelList& masterFc = masterFaceCells();
 
@@ -119,13 +118,13 @@ void Foam::slidingInterface::decoupleInterface
     // Slave side
 
     const primitiveFacePatch& slavePatch =
-        mesh.faceZones()[slaveFaceZoneID_.index()]();
+        faceZones[slaveFaceZoneID_.index()]();
 
     const labelList& slavePatchAddr =
-        mesh.faceZones()[slaveFaceZoneID_.index()];
+        faceZones[slaveFaceZoneID_.index()];
 
     const boolList& slavePatchFlip =
-        mesh.faceZones()[slaveFaceZoneID_.index()].flipMap();
+        faceZones[slaveFaceZoneID_.index()].flipMap();
 
     const labelList& slaveFc = slaveFaceCells();
 
@@ -147,16 +146,7 @@ void Foam::slidingInterface::decoupleInterface
         // Recover retired points on the slave side
         forAll(newFace, pointi)
         {
-            Map<label>::const_iterator rpmIter = rpm.find(newFace[pointi]);
-            if (rpmIter != rpm.end())
-            {
-                // Master of retired point; grab its original
-                // Pout<< "Reinstating retired point: " << newFace[pointi]
-                //     << " with old: " << rpm.find(newFace[pointi])()
-                //     << endl;
-
-                newFace[pointi] = rpmIter();
-            }
+            newFace[pointi] = rpm.lookup(newFace[pointi], newFace[pointi]);
         }
 
         ref.setAction
@@ -181,12 +171,10 @@ void Foam::slidingInterface::decoupleInterface
     // Grab the list of faces in the layer
     const labelList& masterStickOuts = masterStickOutFaces();
 
-    forAll(masterStickOuts, facei)
+    for (const label curFaceID : masterStickOuts)
     {
         // Renumber the face and remove additional points
 
-        const label curFaceID = masterStickOuts[facei];
-
         const face& oldFace = faces[curFaceID];
 
         DynamicList<label> newFaceLabels(oldFace.size());
@@ -219,17 +207,18 @@ void Foam::slidingInterface::decoupleInterface
             }
 
             // Get face zone and its flip
-            label modifiedFaceZone = mesh.faceZones().whichZone(curFaceID);
-            bool modifiedFaceZoneFlip = false;
+            const label modifiedFaceZone = faceZones.whichZone(curFaceID);
 
-            if (modifiedFaceZone >= 0)
-            {
-                modifiedFaceZoneFlip =
-                    mesh.faceZones()[modifiedFaceZone].flipMap()
-                    [
-                        mesh.faceZones()[modifiedFaceZone].whichFace(curFaceID)
-                    ];
-            }
+            const bool modifiedFaceZoneFlip =
+            (
+                modifiedFaceZone >= 0
+              ?
+                faceZones[modifiedFaceZone].flipMap()
+                [
+                    faceZones[modifiedFaceZone].whichFace(curFaceID)
+                ]
+              : false
+            );
 
             face newFace;
             newFace.transfer(newFaceLabels);
@@ -265,23 +254,22 @@ void Foam::slidingInterface::decoupleInterface
         primitiveMesh::facesPerCell_*(masterPatch.size() + slavePatch.size())
     );
 
-    forAll(slaveFc, facei)
+    for (const label slaveFci : slaveFc)
     {
-        const labelList& curFaces = cells[slaveFc[facei]];
+        const labelList& curFaces = cells[slaveFci];
 
-        forAll(curFaces, facei)
+        for (const label facei : curFaces)
         {
             // Check if the face belongs to the slave face zone; and
             // if it has been removed; if not add it
             if
             (
-                mesh.faceZones().whichZone(curFaces[facei])
+                faceZones.whichZone(facei)
              != slaveFaceZoneID_.index()
-             && !ref.faceRemoved(curFaces[facei])
-
+             && !ref.faceRemoved(facei)
             )
             {
-                slaveLayerCellFaceMap.insert(curFaces[facei]);
+                slaveLayerCellFaceMap.insert(facei);
             }
         }
     }
@@ -292,12 +280,10 @@ void Foam::slidingInterface::decoupleInterface
     // Grab master point mapping
     const Map<label>& masterPm = masterPatch.meshPointMap();
 
-    forAll(slaveStickOuts, facei)
+    for (const label curFaceID : slaveStickOuts)
     {
         // Renumber the face and remove additional points
 
-        const label curFaceID = slaveStickOuts[facei];
-
         const face& oldFace = faces[curFaceID];
 
         DynamicList<label> newFaceLabels(oldFace.size());
@@ -307,16 +293,19 @@ void Foam::slidingInterface::decoupleInterface
         forAll(oldFace, pointi)
         {
             // Check if the point is removed or retired
-            if (rpm.found(oldFace[pointi]))
+
+            const label retiredPointi = rpm.lookup(oldFace[pointi], -1);
+
+            if (retiredPointi != -1)
             {
                 // Master of retired point; grab its original
                 changed = true;
 
                 // Pout<< "Reinstating retired point: " << oldFace[pointi]
-                //     << " with old: " << rpm.find(oldFace[pointi])()
+                //     << " with old: " << retiredPointi
                 //     << endl;
 
-                newFaceLabels.append(rpm.find(oldFace[pointi])());
+                newFaceLabels.append(retiredPointi);
             }
             else if (ref.pointRemoved(oldFace[pointi]))
             {
@@ -346,17 +335,19 @@ void Foam::slidingInterface::decoupleInterface
             }
 
             // Get face zone and its flip
-            label modifiedFaceZone = mesh.faceZones().whichZone(curFaceID);
-            bool modifiedFaceZoneFlip = false;
+            const label modifiedFaceZone =
+                faceZones.whichZone(curFaceID);
 
-            if (modifiedFaceZone >= 0)
-            {
-                modifiedFaceZoneFlip =
-                    mesh.faceZones()[modifiedFaceZone].flipMap()
-                    [
-                        mesh.faceZones()[modifiedFaceZone].whichFace(curFaceID)
-                    ];
-            }
+            const bool modifiedFaceZoneFlip =
+            (
+                modifiedFaceZone >= 0
+              ?
+                faceZones[modifiedFaceZone].flipMap()
+                [
+                    faceZones[modifiedFaceZone].whichFace(curFaceID)
+                ]
+              : false
+            );
 
             face newFace;
             newFace.transfer(newFaceLabels);
@@ -386,22 +377,20 @@ void Foam::slidingInterface::decoupleInterface
     }
 
     // Bring all slave patch points back to life
-    const pointField& points = mesh.points();
-
     const labelList& slaveMeshPoints =
-        mesh.faceZones()[slaveFaceZoneID_.index()]().meshPoints();
+        faceZones[slaveFaceZoneID_.index()]().meshPoints();
 
-    forAll(slaveMeshPoints, pointi)
+    for (const label slavePointi : slaveMeshPoints)
     {
         ref.setAction
         (
             polyModifyPoint
             (
-                slaveMeshPoints[pointi],             // point ID
-                points[slaveMeshPoints[pointi]],     // point
-                false,                               // remove from zone
-                mesh.pointZones().whichZone(slaveMeshPoints[pointi]), // zone
-                true                                // in a cell
+                slavePointi,            // point ID
+                points[slavePointi],    // point
+                false,                  // remove from zone
+                mesh.pointZones().whichZone(slavePointi), // zone
+                true                    // in a cell
             )
         );
     }
@@ -414,9 +403,8 @@ void Foam::slidingInterface::decoupleInterface
 
     if (debug)
     {
-        Pout<< "void slidingInterface::coupleInterface("
-            << "polyTopoChange& ref) const : "
-            << "Finished decoupling sliding interface " << name() << endl;
+        Pout<< FUNCTION_NAME << nl
+            << ": Finished decoupling sliding interface " << name() << endl;
     }
 }
 
diff --git a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatch.C b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatch.C
index e0ba1fd8e2..e4df4bc133 100644
--- a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatch.C
+++ b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatch.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  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,7 +47,7 @@ void Foam::enrichedPatch::calcMeshPoints() const
             << abort(FatalError);
     }
 
-    meshPointsPtr_ = new labelList(pointMap().toc());
+    meshPointsPtr_ = new labelList(std::move(pointMap().toc()));
     labelList& mp = *meshPointsPtr_;
 
     sort(mp);
@@ -88,7 +88,7 @@ void Foam::enrichedPatch::calcLocalFaces() const
 
         forAll(f, pointi)
         {
-            curlf[pointi] = mpLookup.find(f[pointi])();
+            curlf[pointi] = mpLookup.cfind(f[pointi])();
         }
     }
 }
@@ -110,7 +110,7 @@ void Foam::enrichedPatch::calcLocalPoints() const
 
     forAll(lp, i)
     {
-        lp[i] = pointMap().find(mp[i])();
+        lp[i] = pointMap().cfind(mp[i])();
     }
 }
 
@@ -131,14 +131,13 @@ void Foam::enrichedPatch::clearOut()
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::enrichedPatch::enrichedPatch
 (
     const primitiveFacePatch& masterPatch,
     const primitiveFacePatch& slavePatch,
-    const labelList& slavePointPointHits,
-    const labelList& slavePointEdgeHits,
-    const List<objectHit>& slavePointFaceHits
+    const labelUList& slavePointPointHits,
+    const labelUList& slavePointEdgeHits,
+    const UList<objectHit>& slavePointFaceHits
 )
 :
     masterPatch_(masterPatch),
@@ -252,36 +251,20 @@ void Foam::enrichedPatch::writeOBJ(const fileName& fName) const
 {
     OFstream str(fName);
 
-    const pointField& lp = localPoints();
-
-    forAll(lp, pointi)
-    {
-        meshTools::writeOBJ(str, lp[pointi]);
-    }
+    meshTools::writeOBJ(str, localPoints());
 
     const faceList& faces = localFaces();
 
-    forAll(faces, facei)
+    for (const face& f : faces)
     {
-        const face& f = faces[facei];
-
         str << 'f';
-        forAll(f, fp)
+        for (const label fp : f)
         {
-            str << ' ' << f[fp]+1;
+            str << ' ' << fp+1;
         }
         str << nl;
     }
 }
 
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //
diff --git a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatch.H b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatch.H
index f47d46940a..6d04a55a33 100644
--- a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatch.H
+++ b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatch.H
@@ -36,8 +36,7 @@ Description
 
     Note:
     If new points are created during master-slave edge cutting, they
-    should be registred with the pointMap.
-
+    should be registered with the pointMap.
 
 SourceFiles
     enrichedPatch.C
@@ -134,10 +133,10 @@ class enrichedPatch
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
-        enrichedPatch(const enrichedPatch&);
+        enrichedPatch(const enrichedPatch&) = delete;
 
         //- Disallow default bitwise assignment
-        void operator=(const enrichedPatch&);
+        void operator=(const enrichedPatch&) = delete;
 
         // Creation of demand-driven private data
 
@@ -196,11 +195,11 @@ public:
         (
             const primitiveFacePatch& masterPatch,
             const primitiveFacePatch& slavePatch,
-            const labelList& slavePointPointHits,
+            const labelUList& slavePointPointHits,
                                             // -1 or common point snapped to
-            const labelList& slavePointEdgeHits,
+            const labelUList& slavePointEdgeHits,
                                             // -1 or common edge snapped to
-            const List<objectHit>& slavePointFaceHits
+            const UList<objectHit>& slavePointFaceHits
                                             // master face snapped to
         );
 
@@ -213,23 +212,17 @@ public:
 
         // Access
 
-            //- Return non-const access to point map to add points
-            Map<point>& pointMap();
-
             //- Return map of points
-            const Map<point>& pointMap() const;
+            inline const Map<point>& pointMap() const;
+
+            //- Return non-const access to point map to add points
+            inline Map<point>& pointMap();
 
             //- Return map of point merges
-            Map<label>& pointMergeMap()
-            {
-                return pointMergeMap_;
-            }
+            inline const Map<label>& pointMergeMap() const;
 
             //- Return map of point merges
-            const Map<label>& pointMergeMap() const
-            {
-                return pointMergeMap_;
-            }
+            inline Map<label>& pointMergeMap();
 
 
         // Topological data
@@ -278,7 +271,7 @@ public:
 
 
         //- Debugging: dump graphical representation to obj format file
-        void writeOBJ(const fileName&) const;
+        void writeOBJ(const fileName& fName) const;
 };
 
 
@@ -288,6 +281,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "enrichedPatchI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //
diff --git a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchCutFaces.C b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchCutFaces.C
index 50c1af29f3..4c1fd31096 100644
--- a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchCutFaces.C
+++ b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchCutFaces.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  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,7 +32,7 @@ Description
 #include "DynamicList.H"
 #include "labelPair.H"
 #include "primitiveMesh.H"
-#include "HashSet.H"
+#include "edgeHashes.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -94,9 +94,8 @@ void Foam::enrichedPatch::calcCutFaces() const
     //    the points projected onto the face.
 
     // Create a set of edge usage parameters
-    HashSet<edge, Hash<edge>> edgesUsedOnce(pp.size());
-    HashSet<edge, Hash<edge>> edgesUsedTwice
-        (pp.size()*primitiveMesh::edgesPerPoint_);
+    edgeHashSet edgesUsedOnce(pp.size());
+    edgeHashSet edgesUsedTwice(pp.size()*primitiveMesh::edgesPerPoint_);
 
 
     forAll(lf, facei)
@@ -145,9 +144,9 @@ void Foam::enrichedPatch::calcCutFaces() const
 
         // Insert the edges of current face into the seed list.
         edgeList cfe = curLocalFace.edges();
-        forAll(curLocalFace, edgeI)
+        for (const edge& e : cfe)
         {
-            edgeSeeds.append(cfe[edgeI]);
+            edgeSeeds.append(e);
         }
 
         // Grab face normal
@@ -156,9 +155,7 @@ void Foam::enrichedPatch::calcCutFaces() const
 
         while (edgeSeeds.size())
         {
-            // Pout<< "edgeSeeds.size(): "
-            //     << edgeSeeds.size()
-            //     << endl;
+            // Pout<< "edgeSeeds.size(): " << edgeSeeds.size() << endl;
 
             const edge curEdge = edgeSeeds.removeHead();
 
@@ -198,6 +195,7 @@ void Foam::enrichedPatch::calcCutFaces() const
             cutFaceGlobalPoints.append(mp[prevPointLabel]);
             cutFaceLocalPoints.append(prevPointLabel);
             // Pout<< "prevPointLabel: " << mp[prevPointLabel] << endl;
+
             // Grab current point and append it to the list
             label curPointLabel = curEdge.end();
             point curPoint = lp[curPointLabel];
@@ -266,7 +264,7 @@ void Foam::enrichedPatch::calcCutFaces() const
 
                         newDir /= magNewDir;
 
-                        scalar curAtanTurn =
+                        const scalar curAtanTurn =
                             atan2(newDir & right, newDir & ahead);
 
                         // Pout<< " atan: " << curAtanTurn << endl;
@@ -352,19 +350,18 @@ void Foam::enrichedPatch::calcCutFaces() const
                     // If the edge corresponds to a starting face edge,
                     // mark the starting face edge as true
 
-                    forAll(cutFaceLocal, cutI)
+                    forAll(cutFaceLocal, cuti)
                     {
                         const edge curCutFaceEdge
                         (
-                            cutFaceLocal[cutI],
-                            cutFaceLocal.nextLabel(cutI)
+                            cutFaceLocal[cuti],
+                            cutFaceLocal.nextLabel(cuti)
                         );
 
                         // Increment the usage count using two hash sets
-                        HashSet<edge, Hash<edge>>::iterator euoIter =
-                            edgesUsedOnce.find(curCutFaceEdge);
+                        auto euoIter = edgesUsedOnce.find(curCutFaceEdge);
 
-                        if (euoIter == edgesUsedOnce.end())
+                        if (!euoIter.found())
                         {
                             // Pout<< "Found edge not used before: "
                             //     << curCutFaceEdge
@@ -435,12 +432,12 @@ void Foam::enrichedPatch::calcCutFaces() const
 
                     if (facei < slavePatch_.size())
                     {
-                        Map<labelList>::const_iterator mpfAddrIter =
-                            masterPointFaceAddr.find(cutFaceGlobal[0]);
+                        const auto mpfAddrIter =
+                            masterPointFaceAddr.cfind(cutFaceGlobal[0]);
 
                         bool otherSideFound = false;
 
-                        if (mpfAddrIter != masterPointFaceAddr.end())
+                        if (mpfAddrIter.found())
                         {
                             bool miss = false;
 
@@ -456,18 +453,13 @@ void Foam::enrichedPatch::calcCutFaces() const
                                 pointi++
                             )
                             {
-                                Map<labelList>::const_iterator
-                                    mpfAddrPointIter =
-                                        masterPointFaceAddr.find
-                                        (
-                                            cutFaceGlobal[pointi]
-                                        );
+                                const auto mpfAddrPointIter =
+                                    masterPointFaceAddr.cfind
+                                    (
+                                        cutFaceGlobal[pointi]
+                                    );
 
-                                if
-                                (
-                                    mpfAddrPointIter
-                                 == masterPointFaceAddr.end()
-                                )
+                                if (!mpfAddrPointIter.found())
                                 {
                                     // Point is off the master patch. Skip
                                     miss = true;
@@ -479,15 +471,11 @@ void Foam::enrichedPatch::calcCutFaces() const
 
                                 // For every current face, try to find it in the
                                 // zero-list
-                                forAll(curMasterFaces, i)
+                                for (const label facei : curMasterFaces)
                                 {
                                     forAll(masterFacesOfPZero, j)
                                     {
-                                        if
-                                        (
-                                            curMasterFaces[i]
-                                         == masterFacesOfPZero[j]
-                                        )
+                                        if (facei == masterFacesOfPZero[j])
                                         {
                                             hits[j]++;
                                             break;
@@ -644,7 +632,7 @@ void Foam::enrichedPatch::calcCutFaces() const
 
     } // end of local faces
 
-    // Re-pack the list into compact storage
+    // Re-pack lists into compact storage
     cutFacesPtr_ = new faceList();
     cutFacesPtr_->transfer(cf);
 
diff --git a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchFaces.C b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchFaces.C
index f6eb9de42e..af4947e912 100644
--- a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchFaces.C
+++ b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchFaces.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  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,11 +26,11 @@ License
 #include "enrichedPatch.H"
 #include "DynamicList.H"
 
-
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 const Foam::label Foam::enrichedPatch::enrichedFaceRatio_ = 3;
 
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::enrichedPatch::calcEnrichedFaces
@@ -79,68 +79,57 @@ void Foam::enrichedPatch::calcEnrichedFaces
 
     forAll(slavePatch_, facei)
     {
-        const face oldFace = slavePatch_[facei];
-        const face oldLocalFace = slaveLocalFaces[facei];
-//         Info<< "old slave face " << facei << ": " << oldFace << endl;
+        const face& oldFace = slavePatch_[facei];
+        const face& oldLocalFace = slaveLocalFaces[facei];
         const labelList& curEdges = slaveFaceEdges[facei];
 
+        // Info<< "old slave face[" << facei << "] " << oldFace << endl;
+
         DynamicList<label> newFace(oldFace.size()*enrichedFaceRatio_);
 
         // Note: The number of points and edges in a face is always identical
         // so both can be done is the same loop
         forAll(oldFace, i)
         {
-            // Add the point
-            Map<label>::const_iterator mpIter =
-                pmm.find(oldFace[i]);
+            // Add the point.
+            // Using the mapped point id if possible
 
-            if (mpIter == pmm.end())
-            {
-                // Point not mapped
-                newFace.append(oldFace[i]);
+            const label mappedPointi = pmm.lookup(oldFace[i], oldFace[i]);
 
-                // Add the projected point into the patch support
-                pointMap().insert
-                (
-                    oldFace[i],    // Global label of point
-                    projectedSlavePoints[oldLocalFace[i]] // Projected position
-                );
-            }
-            else
-            {
-                // Point mapped
-                newFace.append(mpIter());
+            newFace.append(mappedPointi);
 
-                // Add the projected point into the patch support
-                pointMap().insert
-                (
-                    mpIter(),    // Merged global label of point
-                    projectedSlavePoints[oldLocalFace[i]] // Projected position
-                );
-            }
+            // Add the projected point into the patch support
+            pointMap().insert
+            (
+                mappedPointi,   // Global label of point
+                projectedSlavePoints[oldLocalFace[i]] // Projected position
+            );
 
             // Grab the edge points
 
-            const labelList& slavePointsOnEdge =
+            const labelList& pointsOnEdge =
                 pointsIntoSlaveEdges[curEdges[i]];
 
-            // Info<< "slavePointsOnEdge for "
-            //     << curEdges[i] << ": " << slavePointsOnEdge
+            // Info<< "slave pointsOnEdge for "
+            //     << curEdges[i] << ": " << pointsOnEdge
             //     << endl;
 
             // If there are no points on the edge, skip everything
             // If there is only one point, no need for sorting
-            if (slavePointsOnEdge.size())
+            if (pointsOnEdge.size())
             {
                 // Sort edge points in order
-                scalarField edgePointWeights(slavePointsOnEdge.size());
-                const point& startPoint = projectedSlavePoints[oldLocalFace[i]];
+                scalarField weightOnEdge(pointsOnEdge.size());
 
-                vector e =
-                    projectedSlavePoints[oldLocalFace.nextLabel(i)]
-                  - startPoint;
+                const point& startPoint =
+                    projectedSlavePoints[oldLocalFace[i]];
 
-                scalar magSqrE = magSqr(e);
+                const point& endPoint =
+                    projectedSlavePoints[oldLocalFace.nextLabel(i)];
+
+                vector e = (endPoint - startPoint);
+
+                const scalar magSqrE = magSqr(e);
 
                 if (magSqrE > SMALL)
                 {
@@ -154,27 +143,34 @@ void Foam::enrichedPatch::calcEnrichedFaces
                         << abort(FatalError);
                 }
 
-                pointField slavePosOnEdge(slavePointsOnEdge.size());
+                pointField positionOnEdge(pointsOnEdge.size());
 
-                forAll(slavePointsOnEdge, edgePointi)
+                forAll(pointsOnEdge, edgePointi)
                 {
-                    slavePosOnEdge[edgePointi] =
-                        pointMap().find(slavePointsOnEdge[edgePointi])();
+                    positionOnEdge[edgePointi] =
+                        pointMap()[pointsOnEdge[edgePointi]];
 
-                    edgePointWeights[edgePointi] =
-                        (e & (slavePosOnEdge[edgePointi] - startPoint));
+                    weightOnEdge[edgePointi] =
+                    (
+                        e
+                      &
+                        (
+                            positionOnEdge[edgePointi]
+                          - startPoint
+                        )
+                    );
                 }
 
                 if (debug)
                 {
                     // Check weights: all new points should be on the edge
-                    if (min(edgePointWeights) < 0 || max(edgePointWeights) > 1)
+                    if (min(weightOnEdge) < 0 || max(weightOnEdge) > 1)
                     {
                         FatalErrorInFunction
                             << " not on the edge for edge " << curEdges[i]
                             << " of face " << facei << " in slave patch." << nl
-                            << "Min weight: " << min(edgePointWeights)
-                            << " Max weight: " << max(edgePointWeights)
+                            << "Min weight: " << min(weightOnEdge)
+                            << " Max weight: " << max(weightOnEdge)
                             << abort(FatalError);
                     }
                 }
@@ -182,101 +178,98 @@ void Foam::enrichedPatch::calcEnrichedFaces
                 // Go through the points and collect them based on
                 // weights from lower to higher.  This gives the
                 // correct order of points along the edge.
-                forAll(edgePointWeights, passI)
+                forAll(weightOnEdge, passI)
                 {
                     // Max weight can only be one, so the sorting is
                     // done by elimination.
                     label nextPoint = -1;
                     scalar dist = 2;
 
-                    forAll(edgePointWeights, wI)
+                    forAll(weightOnEdge, wI)
                     {
-                        if (edgePointWeights[wI] < dist)
+                        if (weightOnEdge[wI] < dist)
                         {
-                            dist = edgePointWeights[wI];
+                            dist = weightOnEdge[wI];
                             nextPoint = wI;
                         }
                     }
 
                     // Insert the next point and reset its weight to exclude it
                     // from future picks
-                    newFace.append(slavePointsOnEdge[nextPoint]);
-                    edgePointWeights[nextPoint] = GREAT;
+                    newFace.append(pointsOnEdge[nextPoint]);
+                    weightOnEdge[nextPoint] = GREAT;
 
                     // Add the point into patch support
                     pointMap().insert
                     (
-                        slavePointsOnEdge[nextPoint],
-                        slavePosOnEdge[nextPoint]
+                        pointsOnEdge[nextPoint],
+                        positionOnEdge[nextPoint]
                     );
                 }
             }
         }
-        // Info<< "New slave face " << facei << ": " << newFace << endl;
+
+        // Info<< "New slave face[" << facei << "] "
+        //     << flatOutput(newFace) << " was " << flatOutput(oldFace)
+        //     << endl;
 
         // Add the new face to the list
         enrichedFaces[nEnrichedFaces].transfer(newFace);
         nEnrichedFaces++;
     }
 
+
     // Add master faces into the enriched faces list
 
     forAll(masterPatch_, facei)
     {
         const face& oldFace = masterPatch_[facei];
         const face& oldLocalFace = masterLocalFaces[facei];
-//         Info<< "old master face: " << oldFace << endl;
         const labelList& curEdges = masterFaceEdges[facei];
 
+        // Info<< "old master face[" << facei << "] " << oldFace << endl;
+
         DynamicList<label> newFace(oldFace.size()*enrichedFaceRatio_);
 
         // Note: The number of points and edges in a face is always identical
         // so both can be done is the same loop
         forAll(oldFace, i)
         {
-            // Add the point
-            Map<label>::const_iterator mpIter =
-                pmm.find(oldFace[i]);
+            // Add the point.
+            // Using the mapped point id if possible
 
-            if (mpIter == pmm.end())
-            {
-                // Point not mapped
-                newFace.append(oldFace[i]);
+            const label mappedPointi = pmm.lookup(oldFace[i], oldFace[i]);
 
-                // Add the point into patch support
-                pointMap().insert
-                (
-                    oldFace[i],
-                    masterLocalPoints[oldLocalFace[i]]
-                );
-            }
-            else
-            {
-                // Point mapped
-                newFace.append(mpIter());
+            newFace.append(mappedPointi);
 
-                // Add the point into support
-                pointMap().insert(mpIter(), masterLocalPoints[oldLocalFace[i]]);
-            }
+            // Add the point into patch support
+            pointMap().insert
+            (
+                mappedPointi,   // Global label of point
+                masterLocalPoints[oldLocalFace[i]]
+            );
 
             // Grab the edge points
 
-            const labelList& masterPointsOnEdge =
+            const labelList& pointsOnEdge =
                 pointsIntoMasterEdges[curEdges[i]];
 
             // If there are no points on the edge, skip everything
             // If there is only one point, no need for sorting
-            if (masterPointsOnEdge.size())
+            if (pointsOnEdge.size())
             {
                 // Sort edge points in order
-                scalarField edgePointWeights(masterPointsOnEdge.size());
-                const point& startPoint = masterLocalPoints[oldLocalFace[i]];
+                scalarField weightOnEdge(pointsOnEdge.size());
 
-                vector e =
-                    masterLocalPoints[oldLocalFace.nextLabel(i)]
-                  - startPoint;
+                const point& startPoint =
+                    masterLocalPoints[oldLocalFace[i]];
 
-                scalar magSqrE = magSqr(e);
+                const point& endPoint =
+                    masterLocalPoints[oldLocalFace.nextLabel(i)];
+
+                vector e = (endPoint - startPoint);
+
+                const scalar magSqrE = magSqr(e);
 
                 if (magSqrE > SMALL)
                 {
@@ -290,27 +283,33 @@ void Foam::enrichedPatch::calcEnrichedFaces
                         << abort(FatalError);
                 }
 
-                pointField masterPosOnEdge(masterPointsOnEdge.size());
+                pointField positionOnEdge(pointsOnEdge.size());
 
-                forAll(masterPointsOnEdge, edgePointi)
+                forAll(pointsOnEdge, edgePointi)
                 {
-                    masterPosOnEdge[edgePointi] =
-                        pointMap().find(masterPointsOnEdge[edgePointi])();
+                    positionOnEdge[edgePointi] =
+                        pointMap()[pointsOnEdge[edgePointi]];
 
-                    edgePointWeights[edgePointi] =
-                        (e & (masterPosOnEdge[edgePointi] - startPoint));
+                    weightOnEdge[edgePointi] =
+                    (
+                        e
+                      &
+                        (
+                            positionOnEdge[edgePointi] - startPoint
+                        )
+                    );
                 }
 
                 if (debug)
                 {
                     // Check weights: all new points should be on the edge
-                    if (min(edgePointWeights) < 0 || max(edgePointWeights) > 1)
+                    if (min(weightOnEdge) < 0 || max(weightOnEdge) > 1)
                     {
                         FatalErrorInFunction
                             << " not on the edge for edge " << curEdges[i]
                             << " of face " << facei << " in master patch." << nl
-                            << "Min weight: " << min(edgePointWeights)
-                            << " Max weight: " << max(edgePointWeights)
+                            << "Min weight: " << min(weightOnEdge)
+                            << " Max weight: " << max(weightOnEdge)
                             << abort(FatalError);
                     }
                 }
@@ -318,37 +317,40 @@ void Foam::enrichedPatch::calcEnrichedFaces
                 // Go through the points and collect them based on
                 // weights from lower to higher.  This gives the
                 // correct order of points along the edge.
-                forAll(edgePointWeights, passI)
+                forAll(weightOnEdge, passI)
                 {
                     // Max weight can only be one, so the sorting is
                     // done by elimination.
                     label nextPoint = -1;
                     scalar dist = 2;
 
-                    forAll(edgePointWeights, wI)
+                    forAll(weightOnEdge, wI)
                     {
-                        if (edgePointWeights[wI] < dist)
+                        if (weightOnEdge[wI] < dist)
                         {
-                            dist = edgePointWeights[wI];
+                            dist = weightOnEdge[wI];
                             nextPoint = wI;
                         }
                     }
 
                     // Insert the next point and reset its weight to exclude it
                     // from future picks
-                    newFace.append(masterPointsOnEdge[nextPoint]);
-                    edgePointWeights[nextPoint] = GREAT;
+                    newFace.append(pointsOnEdge[nextPoint]);
+                    weightOnEdge[nextPoint] = GREAT;
 
                     // Add the point into patch support
                     pointMap().insert
                     (
-                        masterPointsOnEdge[nextPoint],
-                        masterPosOnEdge[nextPoint]
+                        pointsOnEdge[nextPoint],
+                        positionOnEdge[nextPoint]
                     );
                 }
             }
         }
-        // Info<< "New master face: " << newFace << endl;
+
+        // Info<< "New master face[" << facei << "] "
+        //     << flatOutput(newFace) << " was " << flatOutput(oldFace)
+        //     << endl;
 
         // Add the new face to the list
         enrichedFaces[nEnrichedFaces].transfer(newFace);
diff --git a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchI.H b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchI.H
new file mode 100644
index 0000000000..0a03eb2bd2
--- /dev/null
+++ b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchI.H
@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::Map<Foam::point>& Foam::enrichedPatch::pointMap() const
+{
+    if (!pointMapComplete_)
+    {
+        completePointMap();
+    }
+
+    return pointMap_;
+}
+
+
+inline Foam::Map<Foam::point>& Foam::enrichedPatch::pointMap()
+{
+    if (!pointMapComplete_)
+    {
+        completePointMap();
+    }
+
+    return pointMap_;
+}
+
+
+inline const Foam::Map<Foam::label>& Foam::enrichedPatch::pointMergeMap() const
+{
+    return pointMergeMap_;
+}
+
+
+inline Foam::Map<Foam::label>& Foam::enrichedPatch::pointMergeMap()
+{
+    return pointMergeMap_;
+}
+
+
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchMasterPoints.C b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchMasterPoints.C
index 280f4324d2..18edef858e 100644
--- a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchMasterPoints.C
+++ b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchMasterPoints.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  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,6 +32,7 @@ License
 
 const Foam::label Foam::enrichedPatch::nFaceHits_ = 4;
 
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::enrichedPatch::calcMasterPointFaces() const
@@ -51,7 +52,7 @@ void Foam::enrichedPatch::calcMasterPointFaces() const
     // Master face points lists the points of the enriched master face plus
     // points projected into the master face
 
-    Map<DynamicList<label>> mpf(meshPoints().size());
+    Map<DynamicList<label>> mpf(2*meshPoints().size());
 
     const faceList& ef = enrichedFaces();
 
@@ -59,28 +60,14 @@ void Foam::enrichedPatch::calcMasterPointFaces() const
     forAll(masterPatch_, facei)
     {
         const face& curFace = ef[facei + slavePatch_.size()];
-//         Pout<< "Cur face in pfAddr: " << curFace << endl;
-        forAll(curFace, pointi)
+
+        for (const label pointi : curFace)
         {
-            Map<DynamicList<label>>::iterator mpfIter =
-                mpf.find(curFace[pointi]);
+            DynamicList<label>& dynLst = mpf(pointi); // Get or create
 
-            if (mpfIter == mpf.end())
-            {
-                // Not found, add new dynamic list
-                mpf.insert
-                (
-                    curFace[pointi],
-                    DynamicList<label>(primitiveMesh::facesPerPoint_)
-                );
+            dynLst.reserve(primitiveMesh::facesPerPoint_); // Min size for list
 
-                // Iterator is invalidated - have to find again
-                mpf.find(curFace[pointi])().append(facei);
-            }
-            else
-            {
-                mpfIter().append(facei);
-            }
+            dynLst.append(facei);
         }
     }
 
@@ -96,50 +83,27 @@ void Foam::enrichedPatch::calcMasterPointFaces() const
          && slavePointFaceHits_[pointi].hit()
         )
         {
-            // Get the index of projected point corresponding to this slave
-            // point
-            const label mergedSmp =
-                pointMergeMap().find(slaveMeshPoints[pointi])();
+            // Index of projected point corresponding to this slave point
+            const label mergedPointi = pointMergeMap()[slaveMeshPoints[pointi]];
 
-            Map<DynamicList<label>>::iterator mpfIter =
-                mpf.find(mergedSmp);
+            DynamicList<label>& dynLst = mpf(mergedPointi); // Get or create
 
-            if (mpfIter == mpf.end())
-            {
-                // Not found, add new dynamic list
-                mpf.insert
-                (
-                    mergedSmp,
-                    DynamicList<label>(primitiveMesh::facesPerPoint_)
-                );
+            dynLst.reserve(primitiveMesh::facesPerPoint_); // Min size for list
 
-                // Iterator is invalidated - have to find again
-                mpf.find(mergedSmp)().append
-                (
-                    slavePointFaceHits_[pointi].hitObject()
-                );
-            }
-            else
-            {
-                mpfIter().append(slavePointFaceHits_[pointi].hitObject());
-            }
+            dynLst.append(slavePointFaceHits_[pointi].hitObject());
         }
     }
 
     // Re-pack dynamic lists into normal lists
-    const labelList mpfToc = mpf.toc();
 
-    masterPointFacesPtr_ = new Map<labelList>(2*mpfToc.size());
-    Map<labelList>& masterPointFaceAddr = *masterPointFacesPtr_;
+    masterPointFacesPtr_ = new Map<labelList>(2*mpf.size());
+    Map<labelList>& masterPointFaceMap = *masterPointFacesPtr_;
 
-    forAll(mpfToc, mpfTocI)
+    forAllIters(mpf, mpfIter)
     {
-        labelList l;
-        l.transfer(mpf.find(mpfToc[mpfTocI])());
-
-        masterPointFaceAddr.insert(mpfToc[mpfTocI], l);
+        masterPointFaceMap(mpfIter.key()).transfer(mpfIter.object());
     }
-    // Pout<< "masterPointFaceAddr: " << masterPointFaceAddr << endl;
+    // Pout<< "masterPointFaceMap: " << masterPointFaceMap << endl;
 }
 
 
diff --git a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchPointMap.C b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchPointMap.C
index bab79f6926..632578586e 100644
--- a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchPointMap.C
+++ b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchPointMap.C
@@ -40,8 +40,8 @@ void Foam::enrichedPatch::completePointMap() const
 
     const Map<label>& pmm = pointMergeMap();
 
-    // Get the mesh points for both patches.  If the point has not been
-    // merged away, add it to the map
+    // Get the mesh points for both patches.
+    // If the point has not been merged away, add it to the map
 
     // Do master patch
     const labelList& masterMeshPoints = masterPatch_.meshPoints();
@@ -77,29 +77,4 @@ void Foam::enrichedPatch::completePointMap() const
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
-Foam::Map<Foam::point>& Foam::enrichedPatch::pointMap()
-{
-    if (!pointMapComplete_)
-    {
-        completePointMap();
-    }
-
-    return pointMap_;
-}
-
-
-const Foam::Map<Foam::point>& Foam::enrichedPatch::pointMap() const
-{
-    if (!pointMapComplete_)
-    {
-        completePointMap();
-    }
-
-    return pointMap_;
-}
-
-
 // ************************************************************************* //
diff --git a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchPointPoints.C b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchPointPoints.C
index 7a04712c02..0bd562692a 100644
--- a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchPointPoints.C
+++ b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchPointPoints.C
@@ -27,9 +27,6 @@ License
 #include "primitiveMesh.H"
 #include "DynamicList.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::enrichedPatch::calcPointPoints() const
@@ -52,16 +49,14 @@ void Foam::enrichedPatch::calcPointPoints() const
 
     bool found = false;
 
-    forAll(lf, facei)
+    for (const face& curFace : lf)
     {
-        const face& curFace = lf[facei];
-
         forAll(curFace, pointi)
         {
             DynamicList<label>& curPp = pp[curFace[pointi]];
 
             // Do next label
-            label next = curFace.nextLabel(pointi);
+            const label next = curFace.nextLabel(pointi);
 
             found = false;
 
@@ -80,7 +75,7 @@ void Foam::enrichedPatch::calcPointPoints() const
             }
 
             // Do previous label
-            label prev = curFace.prevLabel(pointi);
+            const label prev = curFace.prevLabel(pointi);
             found = false;
 
             forAll(curPp, i)
diff --git a/src/dynamicMesh/slidingInterface/slidingInterfaceAttachedAddressing.C b/src/dynamicMesh/slidingInterface/slidingInterfaceAttachedAddressing.C
index 34504028ce..305db3f376 100644
--- a/src/dynamicMesh/slidingInterface/slidingInterfaceAttachedAddressing.C
+++ b/src/dynamicMesh/slidingInterface/slidingInterfaceAttachedAddressing.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  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,7 +34,7 @@ void Foam::slidingInterface::calcAttachedAddressing() const
 {
     if (debug)
     {
-        Pout<< "void Foam::slidingInterface::calcAttachedAddressing() const "
+        Pout<< FUNCTION_NAME
             << " for object " << name() << " : "
             << "Calculating zone face-cell addressing."
             << endl;
@@ -140,64 +140,85 @@ void Foam::slidingInterface::calcAttachedAddressing() const
         // Calculate stick-out faces
         const labelListList& pointFaces = mesh.pointFaces();
 
-        // Master side
-        labelHashSet masterStickOutFaceMap
+        labelHashSet stickOutFaceMap
         (
-            primitiveMesh::facesPerCell_*(masterPatch.size())
+            primitiveMesh::facesPerCell_
+          * max(masterPatch.size(), slavePatch.size())
         );
 
+        // Master side
         const labelList& masterMeshPoints = masterPatch.meshPoints();
 
-        forAll(masterMeshPoints, pointi)
-        {
-            const labelList& curFaces = pointFaces[masterMeshPoints[pointi]];
+        stickOutFaceMap.clear();
 
-            forAll(curFaces, facei)
+        for (const label pointi : masterMeshPoints)
+        {
+            for (const label facei : pointFaces[pointi])
             {
-                // Check if the face belongs to the master face zone;
-                // if not add it
+                const label zoneIdx = faceZones.whichZone(facei);
+
+                // Add if face not already part of master or slave face zone
+                // This handles partially attached faces.
                 if
                 (
-                    faceZones.whichZone(curFaces[facei])
-                 != masterFaceZoneID_.index()
+                    zoneIdx != masterFaceZoneID_.index()
+                 && zoneIdx != slaveFaceZoneID_.index()
                 )
                 {
-                    masterStickOutFaceMap.insert(curFaces[facei]);
+                    stickOutFaceMap.insert(facei);
                 }
             }
         }
 
-        masterStickOutFacesPtr_ = new labelList(masterStickOutFaceMap.toc());
+        // Sort in debug mode for easier diagnostics
+        if (debug)
+        {
+            masterStickOutFacesPtr_ =
+                new labelList(std::move(stickOutFaceMap.sortedToc()));
+        }
+        else
+        {
+            masterStickOutFacesPtr_ =
+                new labelList(std::move(stickOutFaceMap.toc()));
+        }
 
         // Slave side
-        labelHashSet slaveStickOutFaceMap
-        (
-            primitiveMesh::facesPerCell_*(slavePatch.size())
-        );
-
         const labelList& slaveMeshPoints = slavePatch.meshPoints();
 
-        forAll(slaveMeshPoints, pointi)
-        {
-            const labelList& curFaces = pointFaces[slaveMeshPoints[pointi]];
+        stickOutFaceMap.clear();
 
-            forAll(curFaces, facei)
+        for (const label pointi : slaveMeshPoints)
+        {
+            for (const label facei : pointFaces[pointi])
             {
-                // Check if the face belongs to the slave face zone;
-                // if not add it
+                const label zoneIdx = faceZones.whichZone(facei);
+
+                // Add if face not already part of master or slave face zone
+                // This handles partially attached faces.
                 if
                 (
-                    faceZones.whichZone(curFaces[facei])
-                 != slaveFaceZoneID_.index()
+                    zoneIdx != masterFaceZoneID_.index()
+                 && zoneIdx != slaveFaceZoneID_.index()
                 )
                 {
-                    slaveStickOutFaceMap.insert(curFaces[facei]);
+                    stickOutFaceMap.insert(facei);
                 }
             }
         }
 
-        slaveStickOutFacesPtr_ = new labelList(slaveStickOutFaceMap.toc());
+        // Sort in debug mode for easier diagnostics
+        if (debug)
+        {
+            slaveStickOutFacesPtr_ =
+                new labelList(std::move(stickOutFaceMap.sortedToc()));
+        }
+        else
+        {
+            slaveStickOutFacesPtr_ =
+                new labelList(std::move(stickOutFaceMap.toc()));
+        }
 
+        stickOutFaceMap.clear();
 
         // Retired point addressing does not exist at this stage.
         // It will be filled when the interface is coupled.
@@ -208,7 +229,6 @@ void Foam::slidingInterface::calcAttachedAddressing() const
             );
 
         // Ditto for cut point edge map.  This is a rough guess of its size
-        //
         cutPointEdgePairMapPtr_ =
             new Map<Pair<edge>>
             (
@@ -224,7 +244,7 @@ void Foam::slidingInterface::calcAttachedAddressing() const
 
     if (debug)
     {
-        Pout<< "void Foam::slidingInterface::calcAttachedAddressing() const "
+        Pout<< FUNCTION_NAME
             << " for object " << name() << " : "
             << "Finished calculating zone face-cell addressing."
             << endl;
diff --git a/src/dynamicMesh/slidingInterface/slidingInterfaceProjectPoints.C b/src/dynamicMesh/slidingInterface/slidingInterfaceProjectPoints.C
index 7a0cfe1a08..9e52c5f103 100644
--- a/src/dynamicMesh/slidingInterface/slidingInterfaceProjectPoints.C
+++ b/src/dynamicMesh/slidingInterface/slidingInterfaceProjectPoints.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  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -58,8 +58,8 @@ bool Foam::slidingInterface::projectPoints() const
 {
     if (debug)
     {
-        Pout<< "bool slidingInterface::projectPoints() : "
-            << " for object " << name() << " : "
+        Pout<< FUNCTION_NAME << nl
+            << ": for object " << name() << " : "
             << "Projecting slave points onto master surface." << endl;
     }
 
@@ -126,12 +126,11 @@ bool Foam::slidingInterface::projectPoints() const
     scalarField minMasterPointLength(masterLocalPoints.size(), GREAT);
     scalarField minMasterFaceLength(masterPatch.size(), GREAT);
 
-    forAll(masterEdges, edgeI)
+    forAll(masterEdges, edgei)
     {
-        const edge& curEdge = masterEdges[edgeI];
+        const edge& curEdge = masterEdges[edgei];
 
-        const scalar curLength =
-            masterEdges[edgeI].mag(masterLocalPoints);
+        const scalar curLength = curEdge.mag(masterLocalPoints);
 
         // Do points
         minMasterPointLength[curEdge.start()] =
@@ -149,14 +148,14 @@ bool Foam::slidingInterface::projectPoints() const
             );
 
         // Do faces
-        const labelList& curFaces = masterEdgeFaces[edgeI];
+        const labelList& curFaces = masterEdgeFaces[edgei];
 
-        forAll(curFaces, facei)
+        for (const label facei : curFaces)
         {
-            minMasterFaceLength[curFaces[facei]] =
+            minMasterFaceLength[facei] =
                 min
                 (
-                    minMasterFaceLength[curFaces[facei]],
+                    minMasterFaceLength[facei],
                     curLength
                 );
         }
@@ -169,12 +168,11 @@ bool Foam::slidingInterface::projectPoints() const
     scalarField minSlavePointLength(slaveLocalPoints.size(), GREAT);
     scalarField minSlaveFaceLength(slavePatch.size(), GREAT);
 
-    forAll(slaveEdges, edgeI)
+    forAll(slaveEdges, edgei)
     {
-        const edge& curEdge = slaveEdges[edgeI];
+        const edge& curEdge = slaveEdges[edgei];
 
-        const scalar curLength =
-            slaveEdges[edgeI].mag(slaveLocalPoints);
+        const scalar curLength = curEdge.mag(slaveLocalPoints);
 
         // Do points
         minSlavePointLength[curEdge.start()] =
@@ -192,14 +190,14 @@ bool Foam::slidingInterface::projectPoints() const
             );
 
         // Do faces
-        const labelList& curFaces = slaveEdgeFaces[edgeI];
+        const labelList& curFaces = slaveEdgeFaces[edgei];
 
-        forAll(curFaces, facei)
+        for (const label facei : curFaces)
         {
-            minSlaveFaceLength[curFaces[facei]] =
+            minSlaveFaceLength[facei] =
                 min
                 (
-                    minSlaveFaceLength[curFaces[facei]],
+                    minSlaveFaceLength[facei],
                     curLength
                 );
         }
@@ -219,25 +217,15 @@ bool Foam::slidingInterface::projectPoints() const
             projectionAlgo_
         );
 
-//     Pout<< "USING N-SQUARED!!!" << endl;
-//     List<objectHit> slavePointFaceHits =
-//         projectPointsNSquared<face, List, const pointField&>
-//         (
-//             slavePatch,
-//             masterPatch,
-//             slavePointNormals,
-//             projectionAlgo_
-//         );
-
     if (debug)
     {
         label nHits = 0;
 
-        forAll(slavePointFaceHits, pointi)
+        for (const auto& hit : slavePointFaceHits)
         {
-            if (slavePointFaceHits[pointi].hit())
+            if (hit.hit())
             {
-                nHits++;
+                ++nHits;
             }
         }
 
@@ -247,8 +235,7 @@ bool Foam::slidingInterface::projectPoints() const
     }
 
     // Projected slave points are stored for mesh motion correction
-    if (projectedSlavePointsPtr_) delete projectedSlavePointsPtr_;
-
+    deleteDemandDrivenData(projectedSlavePointsPtr_);
     projectedSlavePointsPtr_ =
         new pointField(slavePointFaceHits.size(), Zero);
     pointField& projectedSlavePoints = *projectedSlavePointsPtr_;
@@ -380,29 +367,27 @@ bool Foam::slidingInterface::projectPoints() const
 
         // Check for zero-length edges in slave projection
         scalar minEdgeLength = GREAT;
-        scalar el = 0;
         label nShortEdges = 0;
 
-        forAll(slaveEdges, edgeI)
+        for (const edge& e : slaveEdges)
         {
-            el = slaveEdges[edgeI].mag(projectedSlavePoints);
+            const scalar len = e.mag(projectedSlavePoints);
+            minEdgeLength = min(minEdgeLength, len);
 
-            if (el < SMALL)
+            if (len < SMALL)
             {
                 Pout<< "Point projection problems for edge: "
-                    << slaveEdges[edgeI] << ". Length = " << el
+                    << e << ". Length = " << len
                     << endl;
 
                 nShortEdges++;
             }
-
-            minEdgeLength = min(minEdgeLength, el);
         }
 
         if (nShortEdges > 0)
         {
             FatalErrorInFunction
-                << " short projected edges "
+                << nShortEdges << " short projected edges "
                 << "after adjustment for object " << name()
                 << abort(FatalError);
         }
@@ -439,10 +424,10 @@ bool Foam::slidingInterface::projectPoints() const
     {
         if (slavePointFaceHits[pointi].hit())
         {
-            // Taking a non-const reference so the point can be adjusted
+            // Use non-const reference so the point can be adjusted
             point& curPoint = projectedSlavePoints[pointi];
 
-            // Get the hit face
+            // Get the hit face (on master)
             const face& hitFace =
                 masterLocalFaces[slavePointFaceHits[pointi].hitObject()];
 
@@ -450,10 +435,10 @@ bool Foam::slidingInterface::projectPoints() const
             scalar mergeDist = GREAT;
 
             // Try all point before deciding on best fit.
-            forAll(hitFace, hitPointi)
+            for (const label hitPointi : hitFace)
             {
-                scalar dist =
-                    mag(masterLocalPoints[hitFace[hitPointi]] - curPoint);
+                const scalar dist =
+                    mag(masterLocalPoints[hitPointi] - curPoint);
 
                 // Calculate the tolerance
                 const scalar mergeTol =
@@ -461,12 +446,12 @@ bool Foam::slidingInterface::projectPoints() const
                     min
                     (
                         minSlavePointLength[pointi],
-                        minMasterPointLength[hitFace[hitPointi]]
+                        minMasterPointLength[hitPointi]
                     );
 
                 if (dist < mergeTol && dist < mergeDist)
                 {
-                    mergePoint = hitFace[hitPointi];
+                    mergePoint = hitPointi;
                     mergeDist = dist;
 
 //                     Pout<< "Merging slave point "
@@ -481,12 +466,17 @@ bool Foam::slidingInterface::projectPoints() const
 
             if (mergePoint > -1)
             {
-                // Point is to be merged with master point
+                // Slave point is to be merged with master point.
+                // This may also include a false positive when the two points
+                // already point to the same global point, but this will need
+                // to be addressed by the caller.
                 nMergedPoints++;
 
                 slavePointPointHits[pointi] = mergePoint;
-                curPoint = masterLocalPoints[mergePoint];
                 masterPointPointHits[mergePoint] = pointi;
+
+                // Adjust (snap) slave point
+                curPoint = masterLocalPoints[mergePoint];
             }
         }
     }
@@ -498,20 +488,18 @@ bool Foam::slidingInterface::projectPoints() const
     {
         // Check for zero-length edges in slave projection
         scalar minEdgeLength = GREAT;
-        scalar el = 0;
 
-        forAll(slaveEdges, edgeI)
+        for (const edge& e : slaveEdges)
         {
-            el = slaveEdges[edgeI].mag(projectedSlavePoints);
+            const scalar len = e.mag(projectedSlavePoints);
+            minEdgeLength = min(minEdgeLength, len);
 
-            if (el < SMALL)
+            if (len < SMALL)
             {
                 Pout<< "Point projection problems for edge: "
-                    << slaveEdges[edgeI] << ". Length = " << el
+                    << e << ". Length = " << len
                     << endl;
             }
-
-            minEdgeLength = min(minEdgeLength, el);
         }
 
         if (minEdgeLength < SMALL)
@@ -556,16 +544,16 @@ bool Foam::slidingInterface::projectPoints() const
 
             scalar minDistance = GREAT;
 
-            forAll(hitFaceEdges, edgeI)
+            for (const label edgei : hitFaceEdges)
             {
-                const edge& curEdge = masterEdges[hitFaceEdges[edgeI]];
+                const edge& curEdge = masterEdges[edgei];
 
                 pointHit edgeHit =
                     curEdge.line(masterLocalPoints).nearestDist(curPoint);
 
                 if (edgeHit.hit())
                 {
-                    scalar dist =
+                    const scalar dist =
                         mag(edgeHit.hitPoint() - projectedSlavePoints[pointi]);
 
                     if (dist < mergeTol && dist < minDistance)
@@ -573,7 +561,7 @@ bool Foam::slidingInterface::projectPoints() const
                         // Point is to be moved onto master edge
                         nMovedPoints++;
 
-                        slavePointEdgeHits[pointi] = hitFaceEdges[edgeI];
+                        slavePointEdgeHits[pointi] = edgei;
                         projectedSlavePoints[pointi] = edgeHit.hitPoint();
 
                         minDistance = dist;
@@ -582,7 +570,7 @@ bool Foam::slidingInterface::projectPoints() const
 //                             << slavePatch.meshPoints()[pointi]
 //                             << " (" << pointi
 //                             << ") at " << slaveLocalPoints[pointi]
-//                             << " onto master edge " << hitFaceEdges[edgeI]
+//                             << " onto master edge " << edgei
 //                             << " or ("
 //                             << masterLocalPoints[curEdge.start()]
 //                             << masterLocalPoints[curEdge.end()]
@@ -607,7 +595,7 @@ bool Foam::slidingInterface::projectPoints() const
 
                 forAll(hitMasterEdge, hmeI)
                 {
-                    scalar hmeDist =
+                    const scalar hmeDist =
                         mag(masterLocalPoints[hitMasterEdge[hmeI]] - curPoint);
 
                     // Calculate the tolerance
@@ -660,27 +648,25 @@ bool Foam::slidingInterface::projectPoints() const
 
         // Check for zero-length edges in slave projection
         scalar minEdgeLength = GREAT;
-        scalar el = 0;
 
-        forAll(slaveEdges, edgeI)
+        for (const edge& e : slaveEdges)
         {
-            el = slaveEdges[edgeI].mag(projectedSlavePoints);
+            const scalar len = e.mag(projectedSlavePoints);
+            minEdgeLength = min(minEdgeLength, len);
 
-            if (el < SMALL)
+            if (len < SMALL)
             {
                 Pout<< "Point projection problems for edge: "
-                    << slaveEdges[edgeI] << ". Length = " << el
+                    << e << ". Length = " << len
                     << endl;
             }
-
-            minEdgeLength = min(minEdgeLength, el);
         }
 
         if (minEdgeLength < SMALL)
         {
             FatalErrorInFunction
-            << " after correction for object " << name()
-            << abort(FatalError);
+                << " after correction for object " << name()
+                << abort(FatalError);
         }
     }
 
@@ -747,25 +733,19 @@ bool Foam::slidingInterface::projectPoints() const
 
     labelHashSet addedFaces(2*primitiveMesh::edgesPerFace_);
 
-    forAll(slaveEdges, edgeI)
+    forAll(slaveEdges, edgei)
     {
-        const edge& curEdge = slaveEdges[edgeI];
+        const edge& curEdge = slaveEdges[edgei];
 
         {
-            // Clear the maps
-            curFaceMap.clear();
-            addedFaces.clear();
-
             // Grab the faces for start and end points
             const label startFace =
                 slavePointFaceHits[curEdge.start()].hitObject();
-            const label endFace = slavePointFaceHits[curEdge.end()].hitObject();
 
-            // Insert the start face into the list
-            curFaceMap.insert(startFace);
-            addedFaces.insert(startFace);
+            const label endFace =
+                slavePointFaceHits[curEdge.end()].hitObject();
 
-            // Pout<< "Doing edge " << edgeI
+            // Pout<< "Doing edge " << edgei
             //     << " or " << curEdge
             //     << " start: "
             //     << slavePointFaceHits[curEdge.start()].hitObject()
@@ -774,46 +754,59 @@ bool Foam::slidingInterface::projectPoints() const
             //     << endl;
 
             // If the end face is on the list, the face collection is finished
-            label nSweeps = 0;
             bool completed = false;
 
-            while (nSweeps < edgeFaceEscapeLimit_)
+            if (!completed)
             {
-                nSweeps++;
+                // Forward sweep
 
-                if (addedFaces.found(endFace))
-                {
-                    completed = true;
-                }
-
-                // Add all face neighbours of face in the map
-                const labelList cf = addedFaces.toc();
+                // Clear the maps
+                curFaceMap.clear();
                 addedFaces.clear();
 
-                forAll(cf, cfI)
-                {
-                    const labelList& curNbrs = masterFaceFaces[cf[cfI]];
+                // Insert the start face into the list
+                curFaceMap.insert(startFace);
+                addedFaces.insert(startFace);
 
-                    forAll(curNbrs, nbrI)
+                for
+                (
+                    label nSweeps = 0;
+                    nSweeps < edgeFaceEscapeLimit_;
+                    ++nSweeps
+                )
+                {
+                    completed = addedFaces.found(endFace);
+
+                    // Add all face neighbours of face in the map
+                    const labelList cf(std::move(addedFaces.toc()));
+                    addedFaces.clear();
+
+                    for (const label cfi : cf)
                     {
-                        if (!curFaceMap.found(curNbrs[nbrI]))
+                        const labelList& curNbrs = masterFaceFaces[cfi];
+
+                        for (const label nbri : curNbrs)
                         {
-                            curFaceMap.insert(curNbrs[nbrI]);
-                            addedFaces.insert(curNbrs[nbrI]);
+                            if (curFaceMap.insert(nbri))
+                            {
+                                addedFaces.insert(nbri);
+                            }
                         }
                     }
-                }
 
-                if (completed) break;
+                    if (completed) break;
 
-                if (debug)
-                {
-                    Pout<< ".";
+                    if (debug)
+                    {
+                        Pout<< ".";
+                    }
                 }
             }
 
             if (!completed)
             {
+                // Reverse sweep
+
                 if (debug)
                 {
                     Pout<< "x";
@@ -822,35 +815,32 @@ bool Foam::slidingInterface::projectPoints() const
                 // It is impossible to reach the end from the start, probably
                 // due to disconnected domain.  Do search in opposite direction
 
-                label nReverseSweeps = 0;
-
                 addedFaces.clear();
                 curFaceMap.insert(endFace);
                 addedFaces.insert(endFace);
 
-                while (nReverseSweeps < edgeFaceEscapeLimit_)
+                for
+                (
+                    label nSweeps = 0;
+                    nSweeps < edgeFaceEscapeLimit_;
+                    ++nSweeps
+                )
                 {
-                    nReverseSweeps++;
-
-                    if (addedFaces.found(startFace))
-                    {
-                        completed = true;
-                    }
+                    completed = addedFaces.found(startFace);
 
                     // Add all face neighbours of face in the map
-                    const labelList cf = addedFaces.toc();
+                    const labelList cf(std::move(addedFaces.toc()));
                     addedFaces.clear();
 
-                    forAll(cf, cfI)
+                    for (const label cfi : cf)
                     {
-                        const labelList& curNbrs = masterFaceFaces[cf[cfI]];
+                        const labelList& curNbrs = masterFaceFaces[cfi];
 
-                        forAll(curNbrs, nbrI)
+                        for (const label nbri : curNbrs)
                         {
-                            if (!curFaceMap.found(curNbrs[nbrI]))
+                            if (curFaceMap.insert(nbri))
                             {
-                                curFaceMap.insert(curNbrs[nbrI]);
-                                addedFaces.insert(curNbrs[nbrI]);
+                                addedFaces.insert(nbri);
                             }
                         }
                     }
@@ -887,16 +877,11 @@ bool Foam::slidingInterface::projectPoints() const
                 nFacesPerSlaveEdge_*primitiveMesh::pointsPerFace_
             );
 
-            const labelList curFaces = curFaceMap.toc();
-//             Pout<< "curFaces: " << curFaces << endl;
-            forAll(curFaces, facei)
+            for (const label facei : curFaceMap)
             {
-                const face& f = masterLocalFaces[curFaces[facei]];
+                const face& f = masterLocalFaces[facei];
 
-                forAll(f, pointi)
-                {
-                    curPointMap.insert(f[pointi]);
-                }
+                curPointMap.insert(f);  // Insert all face points
             }
 
             const labelList curMasterPoints = curPointMap.toc();
@@ -912,7 +897,7 @@ bool Foam::slidingInterface::projectPoints() const
             // is used to reject master points out of reach.
             // Calculated as a combination of travel distance in projection and
             // edge length
-            scalar slaveCatchDist =
+            const scalar slaveCatchDist =
                 edgeMasterCatchFraction_*edgeMag
               + 0.5*
                 (
@@ -945,10 +930,8 @@ bool Foam::slidingInterface::projectPoints() const
 
             edgeNormalInPlane /= mag(edgeNormalInPlane);
 
-            forAll(curMasterPoints, pointi)
+            for (const label cmp : curMasterPoints)
             {
-                const label cmp = curMasterPoints[pointi];
-
                 // Skip the current point if the edge start or end has
                 // been adjusted onto in
                 if
@@ -974,11 +957,11 @@ bool Foam::slidingInterface::projectPoints() const
 
                     // Strict checking of slave cut to avoid capturing
                     // end points.
-                    scalar cutOnSlave =
+                    const scalar cutOnSlave =
                         ((edgeLineHit.hitPoint() - edgeLine.start()) & edgeVec)
                         /sqr(edgeMag);
 
-                    scalar distInEdgePlane =
+                    const scalar distInEdgePlane =
                         min
                         (
                             edgeLineHit.distance(),
@@ -1032,14 +1015,14 @@ bool Foam::slidingInterface::projectPoints() const
                         }
 
                         // Snap to point onto edge
-                        masterPointEdgeHits[cmp] = edgeI;
+                        masterPointEdgeHits[cmp] = edgei;
                         masterPointEdgeDist[cmp] = edgeLineHit.distance();
 
 //                         Pout<< "Inserting master point "
 //                             << masterPatch.meshPoints()[cmp]
 //                             << " (" << cmp
 //                             << ") at " << masterLocalPoints[cmp]
-//                             << " into slave edge " << edgeI
+//                             << " into slave edge " << edgei
 //                             << " " << curEdge
 //                             << " cutOnSlave: " << cutOnSlave
 //                             << " distInEdgePlane: " << distInEdgePlane
@@ -1096,16 +1079,20 @@ bool Foam::slidingInterface::projectPoints() const
 
         // Store the addressing arrays and projected points
         deleteDemandDrivenData(slavePointPointHitsPtr_);
-        slavePointPointHitsPtr_ = new labelList(slavePointPointHits);
+        slavePointPointHitsPtr_ = new labelList();
+        slavePointPointHitsPtr_->transfer(slavePointPointHits);
 
         deleteDemandDrivenData(slavePointEdgeHitsPtr_);
-        slavePointEdgeHitsPtr_ = new labelList(slavePointEdgeHits);
+        slavePointEdgeHitsPtr_ = new labelList();
+        slavePointEdgeHitsPtr_->transfer(slavePointEdgeHits);
 
         deleteDemandDrivenData(slavePointFaceHitsPtr_);
-        slavePointFaceHitsPtr_ = new List<objectHit>(slavePointFaceHits);
+        slavePointFaceHitsPtr_ = new List<objectHit>();
+        slavePointFaceHitsPtr_->transfer(slavePointFaceHits);
 
         deleteDemandDrivenData(masterPointEdgeHitsPtr_);
-        masterPointEdgeHitsPtr_ = new labelList(masterPointEdgeHits);
+        masterPointEdgeHitsPtr_ = new labelList();
+        masterPointEdgeHitsPtr_->transfer(masterPointEdgeHits);
 
         if (debug)
         {
@@ -1128,16 +1115,20 @@ bool Foam::slidingInterface::projectPoints() const
         {
             // Must be restart.  Force topological change.
             deleteDemandDrivenData(slavePointPointHitsPtr_);
-            slavePointPointHitsPtr_ = new labelList(slavePointPointHits);
+            slavePointPointHitsPtr_ = new labelList();
+            slavePointPointHitsPtr_->transfer(slavePointPointHits);
 
             deleteDemandDrivenData(slavePointEdgeHitsPtr_);
-            slavePointEdgeHitsPtr_ = new labelList(slavePointEdgeHits);
+            slavePointEdgeHitsPtr_ = new labelList();
+            slavePointEdgeHitsPtr_->transfer(slavePointEdgeHits);
 
             deleteDemandDrivenData(slavePointFaceHitsPtr_);
-            slavePointFaceHitsPtr_ = new List<objectHit>(slavePointFaceHits);
+            slavePointFaceHitsPtr_ = new List<objectHit>();
+            slavePointFaceHitsPtr_->transfer(slavePointFaceHits);
 
             deleteDemandDrivenData(masterPointEdgeHitsPtr_);
-            masterPointEdgeHitsPtr_ = new labelList(masterPointEdgeHits);
+            masterPointEdgeHitsPtr_ = new labelList();
+            masterPointEdgeHitsPtr_->transfer(masterPointEdgeHits);
 
             if (debug)
             {
@@ -1218,16 +1209,20 @@ bool Foam::slidingInterface::projectPoints() const
         {
             // Grab new data
             deleteDemandDrivenData(slavePointPointHitsPtr_);
-            slavePointPointHitsPtr_ = new labelList(slavePointPointHits);
+            slavePointPointHitsPtr_ = new labelList();
+            slavePointPointHitsPtr_->transfer(slavePointPointHits);
 
             deleteDemandDrivenData(slavePointEdgeHitsPtr_);
-            slavePointEdgeHitsPtr_ = new labelList(slavePointEdgeHits);
+            slavePointEdgeHitsPtr_ = new labelList();
+            slavePointEdgeHitsPtr_->transfer(slavePointEdgeHits);
 
             deleteDemandDrivenData(slavePointFaceHitsPtr_);
-            slavePointFaceHitsPtr_ = new List<objectHit>(slavePointFaceHits);
+            slavePointFaceHitsPtr_ = new List<objectHit>();
+            slavePointFaceHitsPtr_->transfer(slavePointFaceHits);
 
             deleteDemandDrivenData(masterPointEdgeHitsPtr_);
-            masterPointEdgeHitsPtr_ = new labelList(masterPointEdgeHits);
+            masterPointEdgeHitsPtr_ = new labelList();
+            masterPointEdgeHitsPtr_->transfer(masterPointEdgeHits);
 
             if (debug)
             {
diff --git a/tutorials/mesh/stitchMesh/simple-cube1/Allclean b/tutorials/mesh/stitchMesh/simple-cube1/Allclean
new file mode 100755
index 0000000000..d6d329b5ff
--- /dev/null
+++ b/tutorials/mesh/stitchMesh/simple-cube1/Allclean
@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1                        # Run from this directory
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions  # Tutorial clean functions
+
+cleanCase0
+rmdir constant 2> /dev/null
+
+# -----------------------------------------------------------------------------
diff --git a/tutorials/mesh/stitchMesh/simple-cube1/Allmesh b/tutorials/mesh/stitchMesh/simple-cube1/Allmesh
new file mode 100755
index 0000000000..a3e4e0da15
--- /dev/null
+++ b/tutorials/mesh/stitchMesh/simple-cube1/Allmesh
@@ -0,0 +1,19 @@
+#!/bin/sh
+cd ${0%/*} || exit 1                        # Run from this directory
+. $WM_PROJECT_DIR/bin/tools/RunFunctions    # Tutorial run functions
+
+runApplication blockMesh
+runApplication topoSet
+
+runApplication subsetMesh final -patch outer -resultTime 0
+
+# Avoid any inconsistencies
+rm -f constant/polyMesh/{cellZones,faceZones,pointZones}
+
+runApplication -s patches topoSet -dict system/topoSetDict.patches
+
+runApplication createPatch -overwrite
+
+# rm -f 0/polyMesh/{cellZones,faceZones,pointZones}
+
+# -----------------------------------------------------------------------------
diff --git a/tutorials/mesh/stitchMesh/simple-cube1/Allrun b/tutorials/mesh/stitchMesh/simple-cube1/Allrun
new file mode 100755
index 0000000000..1f37298922
--- /dev/null
+++ b/tutorials/mesh/stitchMesh/simple-cube1/Allrun
@@ -0,0 +1,12 @@
+#!/bin/sh
+cd ${0%/*} || exit 1                        # Run from this directory
+. $WM_PROJECT_DIR/bin/tools/RunFunctions    # Tutorial run functions
+
+runApplication ./Allmesh
+
+for dir in x y z
+do
+    runApplication -s dir-$dir stitchMesh -partial outer$dir inner$dir
+done
+
+# -----------------------------------------------------------------------------
diff --git a/tutorials/mesh/stitchMesh/simple-cube1/scripts/stitch b/tutorials/mesh/stitchMesh/simple-cube1/scripts/stitch
new file mode 100755
index 0000000000..fd63d6d682
--- /dev/null
+++ b/tutorials/mesh/stitchMesh/simple-cube1/scripts/stitch
@@ -0,0 +1,50 @@
+#!/bin/sh
+cd ${0%/*} || exit 1                        # Run from this directory
+. $WM_PROJECT_DIR/bin/tools/RunFunctions    # Tutorial run functions
+
+[ "$#" -gt 0 ] || {
+    echo "provide x, y, z or -all directions"
+    exit 1
+}
+
+unset optAll
+for dir
+do
+    case "$dir" in
+    -x | x) rm -rf 1 2 3 ;;
+    -y | y) rm -rf 2 3 ;;
+    -z | z) rm -rf 3 ;;
+    -all | all)
+        optAll=true
+        rm -rf 1 2 3
+        ;;
+    *)
+        echo "provide x, y, z or -all directions"
+        exit 1
+        ;;
+    esac
+done
+
+if [ "$optAll" = true ]
+then
+    set -- x y z
+fi
+
+stitch()
+{
+    local dir=$1
+    stitchMesh -partial outer$dir inner$dir | tee log.stitch-$dir
+}
+
+
+for dir
+do
+    case "$dir" in
+    -x | x) [ -d 1 ] || stitch x;;
+    -y | y) [ -d 2 ] || stitch y;;
+    -z | z) [ -d 3 ] || stitch z;;
+    esac
+done
+
+
+# -----------------------------------------------------------------------------
diff --git a/tutorials/mesh/stitchMesh/simple-cube1/system/blockMeshDict b/tutorials/mesh/stitchMesh/simple-cube1/system/blockMeshDict
new file mode 100644
index 0000000000..5c02992369
--- /dev/null
+++ b/tutorials/mesh/stitchMesh/simple-cube1/system/blockMeshDict
@@ -0,0 +1,132 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   1;
+
+// Dimensions for outer box
+lo      0;      // outer lower corner
+uo      1;      // outer upper corner
+
+// Dimensions for inner box
+li      -0.1;   // inner lower corner
+
+ui      0.497;  // inner upper corner
+//      ^^^^^ very poorly handled gaps
+
+ui      0.498;  // inner upper corner
+//      ^^^^^ slave patch probably does not project onto master
+
+ui      0.499;  // inner upper corner
+//      ^^^^^ OK
+
+
+vertices
+(
+    ( $lo $lo $lo )
+    ( $uo $lo $lo )
+    ( $uo $uo $lo )
+    ( $lo $uo $lo )
+    ( $lo $lo $uo )
+    ( $uo $lo $uo )
+    ( $uo $uo $uo )
+    ( $lo $uo $uo )
+
+    ( $li $li $li )
+    ( $ui $li $li )
+    ( $ui $ui $li )
+    ( $li $ui $li )
+    ( $li $li $ui )
+    ( $ui $li $ui )
+    ( $ui $ui $ui )
+    ( $li $ui $ui )
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7)       (6 6 6)    simpleGrading (1 1 1)
+    hex (8 9 10 11 12 13 14 15) (10 10 10) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    wallso
+    {
+        type patch;
+        faces
+        (
+            // Walls of the outside block
+            (0 4 7 3)
+            (1 2 6 5)
+            (0 1 5 4)
+            (3 7 6 2)
+            (0 3 2 1)
+            (4 5 6 7)
+        );
+    }
+    wallsi
+    {
+        type patch;
+        faces
+        (
+            // Walls of the inside block
+            (8 12 15 11)
+            (8 9 13 12)
+            (8 11 10 9)
+        );
+    }
+
+    // Outer faces populated by subsetMesh
+    outer
+    {
+        type patch;
+        faces ();
+    }
+
+    innerx
+    {
+        type patch;
+        inGroups (inner);
+        faces
+        (
+            (9 10 14 13)
+        );
+    }
+    innery
+    {
+        type patch;
+        inGroups (inner);
+        faces
+        (
+            (11 15 14 10)
+        );
+    }
+    innerz
+    {
+        type patch;
+        inGroups (inner);
+        faces
+        (
+            (12 13 14 15)
+        );
+    }
+
+);
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/stitchMesh/simple-cube1/system/controlDict b/tutorials/mesh/stitchMesh/simple-cube1/system/controlDict
new file mode 100644
index 0000000000..b0f6027ead
--- /dev/null
+++ b/tutorials/mesh/stitchMesh/simple-cube1/system/controlDict
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     simpleFoam;
+
+startFrom       latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         10;
+
+deltaT          1;
+
+writeControl    runTime;
+
+writeInterval   1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/stitchMesh/simple-cube1/system/createPatchDict b/tutorials/mesh/stitchMesh/simple-cube1/system/createPatchDict
new file mode 100644
index 0000000000..deae7b1a60
--- /dev/null
+++ b/tutorials/mesh/stitchMesh/simple-cube1/system/createPatchDict
@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      createPatchDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+pointSync   false;
+
+outerSet
+{
+    patchInfo
+    {
+        type patch;
+        inGroups (outer);
+    }
+
+    constructFrom set;
+}
+
+
+// Patches to create.
+patches
+(
+    {
+        name outerx;
+        set  outerx;
+        $outerSet;
+    }
+
+    {
+        name outery;
+        set  outery;
+        $outerSet;
+    }
+
+    {
+        name outerz;
+        set  outerz;
+        $outerSet;
+    }
+);
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/stitchMesh/simple-cube1/system/fvSchemes b/tutorials/mesh/stitchMesh/simple-cube1/system/fvSchemes
new file mode 100644
index 0000000000..9c127712e7
--- /dev/null
+++ b/tutorials/mesh/stitchMesh/simple-cube1/system/fvSchemes
@@ -0,0 +1,55 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      Gauss cubic;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/stitchMesh/simple-cube1/system/fvSolution b/tutorials/mesh/stitchMesh/simple-cube1/system/fvSolution
new file mode 100644
index 0000000000..fce1b29ce9
--- /dev/null
+++ b/tutorials/mesh/stitchMesh/simple-cube1/system/fvSolution
@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-06;
+        relTol          0;
+    }
+
+    U
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-05;
+        relTol          0;
+    }
+}
+
+PISO
+{
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/stitchMesh/simple-cube1/system/topoSetDict b/tutorials/mesh/stitchMesh/simple-cube1/system/topoSetDict
new file mode 100644
index 0000000000..9602c4e147
--- /dev/null
+++ b/tutorials/mesh/stitchMesh/simple-cube1/system/topoSetDict
@@ -0,0 +1,106 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      topoSetDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+    // Outer cube
+    {
+        name    outer;
+        type    cellSet;
+        action  new;
+        source  regionToCell;
+        sourceInfo
+        {
+            insidePoints
+            (
+                (0.99 0.99 0.99)
+            );
+        }
+    }
+
+    {
+        name    corner;
+        type    cellSet;
+        action  new;
+        source  boxToCell;
+        sourceInfo
+        {
+            box (-1 -1 -1)(0.5 0.5 0.5);
+        }
+    }
+
+    {
+        name    inner;
+        type    cellSet;
+        action  new;
+        source  cellToCell;
+        sourceInfo
+        {
+            set corner;
+        }
+    }
+
+    {
+        name    inner;
+        type    cellSet;
+        action  delete;
+        source  cellToCell;
+        sourceInfo
+        {
+            set outer;
+        }
+    }
+
+    {
+        name    outer;
+        type    cellSet;
+        action  delete;
+        source  cellToCell;
+        sourceInfo
+        {
+            set corner;
+        }
+    }
+
+    {
+        name    corner;
+        type    cellSet;
+        action  clear;
+    }
+
+    {
+        name    final;
+        type    cellSet;
+        action  new;
+        source  cellToCell;
+        sourceInfo
+        {
+            set outer;
+        }
+    }
+    {
+        name    final;
+        type    cellSet;
+        action  add;
+        source  cellToCell;
+        sourceInfo
+        {
+            set inner;
+        }
+    }
+);
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/stitchMesh/simple-cube1/system/topoSetDict.patches b/tutorials/mesh/stitchMesh/simple-cube1/system/topoSetDict.patches
new file mode 100644
index 0000000000..36a415b66b
--- /dev/null
+++ b/tutorials/mesh/stitchMesh/simple-cube1/system/topoSetDict.patches
@@ -0,0 +1,88 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      topoSetDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+newFromPatch
+{
+    type    faceSet;
+    action  new;
+    source  patchToFace;
+    sourceInfo
+    {
+        name outer;
+    }
+}
+
+subsetNormal
+{
+    type    faceSet;
+    action  subset;
+    source  normalToFace;
+    sourceInfo
+    {
+        cos     0.1;
+    }
+}
+
+actions
+(
+    // outerx
+    {
+        name  outerx;
+        $newFromPatch
+    }
+
+    {
+        name  outerx;
+        $subsetNormal;
+        sourceInfo
+        {
+            normal (-1 0 0);
+        }
+    }
+
+    // outery
+    {
+        name  outery;
+        $newFromPatch
+    }
+
+    {
+        name  outery;
+        $subsetNormal;
+        sourceInfo
+        {
+            normal (0 -1 0);
+        }
+    }
+
+    // outerz
+    {
+        name  outerz;
+        $newFromPatch
+    }
+
+    {
+        name  outerz;
+        $subsetNormal;
+        sourceInfo
+        {
+            normal (0 0 -1);
+        }
+    }
+
+);
+
+// ************************************************************************* //