diff --git a/src/meshTools/regionSplit/regionSplit.C b/src/meshTools/regionSplit/regionSplit.C
index a901e8d763..eb399eae96 100644
--- a/src/meshTools/regionSplit/regionSplit.C
+++ b/src/meshTools/regionSplit/regionSplit.C
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,12 +24,11 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "regionSplit.H"
-#include "FaceCellWave.H"
 #include "cyclicPolyPatch.H"
 #include "processorPolyPatch.H"
 #include "globalIndex.H"
 #include "syncTools.H"
-#include "minData.H"
+#include "clockValue.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -38,176 +37,361 @@ namespace Foam
     defineTypeNameAndDebug(regionSplit, 0);
 }
 
+static constexpr Foam::label UNASSIGNED = -1;
+static constexpr Foam::label BLOCKED = -2;
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+//- Use key/value from mapper to modify the values for map
+
+template<class Container>
+static void inplaceMapValue(const Map<label>& mapper, Container& input)
+{
+    if (mapper.empty())
+    {
+        return;
+    }
+
+    for (auto iter = input.begin(); iter != input.end(); ++iter)
+    {
+        label& value = *iter;
+
+        auto mapIter = mapper.find(value);
+        if (mapIter.found())
+        {
+            value = *mapIter;
+        }
+    }
+}
+
+
+//- The sizes of a List of containers (eg, labelHashSet)
+template<class Container>
+static labelList containerSizes(const UList<Container>& input)
+{
+    const label len = input.size();
+
+    labelList output(len);
+
+    for (label i = 0; i < len; ++i)
+    {
+        output[i] = input[i].size();
+    }
+
+    return output;
+}
+
+} // end of namespace Foam
+
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::regionSplit::calcNonCompactRegionSplit
+void Foam::regionSplit::updateFacePair
+(
+    const label face0,
+    const label face1,
+
+    labelList& faceRegion,
+    DynamicList<label>& facesChanged
+) const
+{
+    if (faceRegion[face0] == UNASSIGNED)
+    {
+        if (faceRegion[face1] >= 0)
+        {
+            faceRegion[face0] = faceRegion[face1];
+            facesChanged.append(face0);
+        }
+        else if (faceRegion[face1] == BLOCKED)
+        {
+            // face1 blocked but not face0.
+            // - illegal for coupled faces, OK for explicit connections.
+        }
+    }
+    else if (faceRegion[face0] >= 0)
+    {
+        if (faceRegion[face1] == UNASSIGNED)
+        {
+            faceRegion[face1] = faceRegion[face0];
+            facesChanged.append(face1);
+        }
+        else if (faceRegion[face1] == BLOCKED)
+        {
+            // face1 blocked but not face0.
+            // - illegal for coupled faces, OK for explicit connections.
+        }
+        else if (faceRegion[face1] != faceRegion[face0])
+        {
+            FatalErrorInFunction
+                << "Problem : coupled face " << face0
+                << " on patch " << mesh().boundaryMesh().whichPatch(face0)
+                << " has region " << faceRegion[face0]
+                << " but coupled face " << face1
+                << " has region " << faceRegion[face1] << nl
+                << "Is your blocked faces specification"
+                << " synchronized across coupled boundaries?" << endl
+                << abort(FatalError);
+        }
+    }
+}
+
+
+void Foam::regionSplit::fillSeedMask
+(
+    const List<labelPair>& explicitConnections,
+    const label seedCellId,
+    const label markValue,
+    labelList& cellRegion,
+    labelList& faceRegion
+) const
+{
+    // Seed cell
+    cellRegion[seedCellId] = markValue;
+
+    // Faces on seed cell
+    changedFaces_.clear();
+    for (const label facei : mesh().cells()[seedCellId])
+    {
+        if (faceRegion[facei] == UNASSIGNED)
+        {
+            faceRegion[facei] = markValue;
+            changedFaces_.append(facei);
+        }
+    }
+
+    const polyBoundaryMesh& patches = mesh().boundaryMesh();
+
+    // Loop over changed faces. FaceCellWave in small.
+
+    while (changedFaces_.size())
+    {
+        changedCells_.clear();
+
+        for (const label facei : changedFaces_)
+        {
+            const label own = mesh().faceOwner()[facei];
+
+            if (cellRegion[own] == UNASSIGNED)
+            {
+                cellRegion[own] = markValue;
+                changedCells_.append(own);
+            }
+
+            if (mesh().isInternalFace(facei))
+            {
+                const label nei = mesh().faceNeighbour()[facei];
+
+                if (cellRegion[nei] == UNASSIGNED)
+                {
+                    cellRegion[nei] = markValue;
+                    changedCells_.append(nei);
+                }
+            }
+        }
+
+        if (debug & 2)
+        {
+            Pout<< " Changed cells / faces : "
+                << changedCells_.size() << " / " << changedFaces_.size()
+                << " before sync" << endl;
+        }
+
+        // Loop over changedCells and collect faces
+        changedFaces_.clear();
+        for (const label celli : changedCells_)
+        {
+            for (const label facei : mesh().cells()[celli])
+            {
+                if (faceRegion[facei] == UNASSIGNED)
+                {
+                    faceRegion[facei] = markValue;
+                    changedFaces_.append(facei);
+                }
+            }
+        }
+
+
+        // Update locally coupled faces
+        // Global connections are done later.
+
+        for (const polyPatch& pp : patches)
+        {
+            if
+            (
+                isA<cyclicPolyPatch>(pp)
+             && refCast<const cyclicPolyPatch>(pp).owner()
+            )
+            {
+                // Transfer from neighbourPatch to here or vice versa.
+
+                const cyclicPolyPatch& cycPatch =
+                    refCast<const cyclicPolyPatch>(pp);
+
+                label face0 = cycPatch.start();
+
+                forAll(cycPatch, i)
+                {
+                    const label face1 = cycPatch.transformGlobalFace(face0);
+
+                    updateFacePair
+                    (
+                        face0,
+                        face1,
+                        faceRegion,
+                        changedFaces_
+                    );
+
+                    ++face0;
+                }
+            }
+        }
+
+        for (const labelPair& pr : explicitConnections)
+        {
+            updateFacePair
+            (
+                pr.first(),
+                pr.second(),
+                faceRegion,
+                changedFaces_
+            );
+        }
+
+        if (debug & 2)
+        {
+            Pout<< " Changed faces : "
+                << changedFaces_.size()
+                << " after sync" << endl;
+        }
+    }
+}
+
+
+Foam::label Foam::regionSplit::calcLocalRegionSplit
 (
-    const globalIndex& globalFaces,
     const boolList& blockedFace,
     const List<labelPair>& explicitConnections,
 
     labelList& cellRegion
 ) const
 {
-    // Field on cells and faces.
-    List<minData> cellData(mesh().nCells());
-    List<minData> faceData(mesh().nFaces());
+    clockValue timing(debug);
 
-    // Take over blockedFaces by seeding a negative number
-    // (so is always less than the decomposition)
-
-    forAll(blockedFace, facei)
+    if (debug)
     {
-        if (blockedFace[facei])
+        if (blockedFace.size())
         {
-            faceData[facei] = minData(-2);
+            // Check that blockedFace is synced.
+            boolList syncBlockedFace(blockedFace);
+            syncTools::swapFaceList(mesh(), syncBlockedFace);
+
+            forAll(syncBlockedFace, facei)
+            {
+                if (syncBlockedFace[facei] != blockedFace[facei])
+                {
+                    FatalErrorInFunction
+                        << "Face " << facei << " not synchronised. My value:"
+                        << blockedFace[facei] << "  coupled value:"
+                        << syncBlockedFace[facei]
+                        << abort(FatalError);
+                }
+            }
         }
     }
 
-    // Seed all faces on (real) boundaries and cell faces next to blockFace,
-    // since regions can only occur because of boundaries (or blocked faces)
+    changedCells_.reserve(mesh_.nCells());
+    changedFaces_.reserve(mesh_.nFaces());
 
-    bitSet isSeed(mesh().nFaces());
-
-    // Get internal or coupled faces
-    bitSet isConnection(syncTools::getInternalOrCoupledFaces(mesh()));
-
-    // 1. Seed (real) boundaries
-    for
-    (
-        label facei = mesh().nInternalFaces();
-        facei < mesh().nFaces();
-        ++facei
-    )
-    {
-        if (!isConnection[facei])
-        {
-            isSeed.set(facei);
-        }
-    }
-
-    // 2. Seed (internal) faces on cells next to blocked (internal) faces
-    //    (since we can't seed the blocked faces themselves)
+    // Region per face.
+    // -1 = unassigned
+    // -2 = blocked
+    labelList faceRegion(mesh().nFaces(), UNASSIGNED);
 
     if (blockedFace.size())
     {
-        label nBlockedCells = 0;
-
-        label celli = 0;
-
-        for (const cell& cFaces : mesh().cells())
+        forAll(blockedFace, facei)
         {
-            bool blockedCell = false;
-            label connectedFacei = -1;
-
-            for (const label facei : cFaces)
+            if (blockedFace[facei])
             {
-                if (blockedFace[facei])
-                {
-                    blockedCell = true;
-                }
-                else if (isConnection[facei])
-                {
-                    connectedFacei = facei;
-                }
+                faceRegion[facei] = BLOCKED;
             }
-
-            if (blockedCell)
-            {
-                if (connectedFacei < 0)
-                {
-                    ++nBlockedCells;
-
-                    if (debug && nBlockedCells <= 10)
-                    {
-                        // Report a few, but not all
-                        WarningInFunction
-                            << "Cell " << celli << " has all faces blocked."
-                            << endl;
-                    }
-                }
-                else
-                {
-                    isSeed.set(connectedFacei);
-                }
-            }
-
-            ++celli;
-        }
-
-        reduce(nBlockedCells, sumOp<label>());
-
-        if (nBlockedCells)
-        {
-            Info<< "Detected " << nBlockedCells << " from "
-                << returnReduce(mesh().nCells(), sumOp<label>())
-                << " cells with all faces blocked." << nl;
         }
     }
 
-    List<label> seedFaces(isSeed.toc());
-    List<minData> seedData(seedFaces.size());
 
-    // Seed face with globally unique number
-    forAll(seedFaces, i)
+    // Assign local regions
+    // ~~~~~~~~~~~~~~~~~~~~
+
+    // Start with region 0
+    label nLocalRegions = 0;
+
+    for (label seedCellId = 0; seedCellId < cellRegion.size(); ++seedCellId)
     {
-        const label facei = seedFaces[i];
-        seedData[i] = minData(globalFaces.toGlobal(facei));
+        // Find next unset cell - use as seed
+
+        for (; seedCellId < cellRegion.size(); ++seedCellId)
+        {
+            if (cellRegion[seedCellId] == UNASSIGNED)
+            {
+                break;
+            }
+        }
+
+        if (seedCellId >= cellRegion.size())
+        {
+            break;
+        }
+
+        fillSeedMask
+        (
+            explicitConnections,
+            seedCellId,
+            nLocalRegions,
+            cellRegion,
+            faceRegion
+        );
+
+        ++nLocalRegions; // Next region
+    }
+
+    // Discard temporary working data
+    changedCells_.clearStorage();
+    changedFaces_.clearStorage();
+
+    if (debug)
+    {
+        forAll(cellRegion, celli)
+        {
+            if (cellRegion[celli] < 0)
+            {
+                FatalErrorInFunction
+                    << "cell:" << celli << " region:" << cellRegion[celli]
+                    << abort(FatalError);
+            }
+        }
+
+        forAll(faceRegion, facei)
+        {
+            if (faceRegion[facei] == UNASSIGNED)
+            {
+                FatalErrorInFunction
+                    << "face:" << facei << " region:" << faceRegion[facei]
+                    << abort(FatalError);
+            }
+        }
     }
 
     if (debug)
     {
-        Pout<< "Seeded " << seedFaces.size()
-            << " boundary and internal faces out of"
-            << " total:" << mesh().nInternalFaces()
-            << " boundary:" << mesh().nFaces()-mesh().nInternalFaces()
-            << endl;
+        Info<<"regionSplit = " << double(timing.elapsed()) << "s\n";
     }
 
-    // Propagate information inwards
-    FaceCellWave<minData> deltaCalc
-    (
-        mesh(),
-        explicitConnections,
-        false,  // disable walking through cyclicAMI for backwards compatibility
-        seedFaces,
-        seedData,
-        faceData,
-        cellData,
-        mesh().globalData().nTotalCells()+1
-    );
-
-
-    // And extract
-    cellRegion.setSize(mesh().nCells());
-    forAll(cellRegion, celli)
-    {
-        if (cellData[celli].valid(deltaCalc.data()))
-        {
-            cellRegion[celli] = cellData[celli].data();
-        }
-        else
-        {
-            // Unvisited cell -> only possible if surrounded by blocked faces
-            // since:
-            // - all walls become seed faces
-            // - all faces on cells with blocked faces become seed faces
-            // If so make up region from any of the faces
-            const cell& cFaces = mesh().cells()[celli];
-            label facei = cFaces[0];
-
-            if (blockedFace.size() && !blockedFace[facei])
-            {
-                FatalErrorInFunction
-                    << "Problem: unblocked face " << facei
-                    << " at " << mesh().faceCentres()[facei]
-                    << " on unassigned cell " << celli
-                    << mesh().cellCentres()[celli]
-                    << exit(FatalError);
-            }
-            cellRegion[celli] = globalFaces.toGlobal(facei);
-        }
-    }
+    return nLocalRegions;
 }
 
 
@@ -220,261 +404,29 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
     labelList& cellRegion
 ) const
 {
-    // See header in regionSplit.H
-
-
-    if (!doGlobalRegions)
-    {
-        // Block all parallel faces to avoid comms across
-        boolList coupledOrBlockedFace(blockedFace);
-        const polyBoundaryMesh& pbm = mesh().boundaryMesh();
-
-        if (coupledOrBlockedFace.size())
-        {
-            forAll(pbm, patchi)
-            {
-                const polyPatch& pp = pbm[patchi];
-                if (isA<processorPolyPatch>(pp))
-                {
-                    label facei = pp.start();
-                    forAll(pp, i)
-                    {
-                        coupledOrBlockedFace[facei++] = true;
-                    }
-                }
-            }
-        }
-
-        // Create dummy (local only) globalIndex
-        labelList offsets(Pstream::nProcs()+1, 0);
-        for (label i = Pstream::myProcNo()+1; i < offsets.size(); i++)
-        {
-            offsets[i] = mesh().nFaces();
-        }
-        const globalIndex globalRegions(std::move(offsets));
-
-        // Minimise regions across connected cells
-        // Note: still uses global decisions so all processors are running
-        //       in lock-step, i.e. slowest determines overall time.
-        //       To avoid this we could switch off Pstream::parRun.
-        calcNonCompactRegionSplit
-        (
-            globalRegions,
-            coupledOrBlockedFace,
-            explicitConnections,
-            cellRegion
-        );
-
-        // Compact
-        Map<label> globalToCompact(mesh().nCells()/8);
-        forAll(cellRegion, celli)
-        {
-            label region = cellRegion[celli];
-
-            label globalRegion;
-
-            const auto fnd = globalToCompact.cfind(region);
-            if (fnd.found())
-            {
-                globalRegion = *fnd;
-            }
-            else
-            {
-                globalRegion = globalRegions.toGlobal(globalToCompact.size());
-                globalToCompact.insert(region, globalRegion);
-            }
-            cellRegion[celli] = globalRegion;
-        }
-
-
-        // Return globalIndex with size = localSize and all regions local
-        labelList compactOffsets(Pstream::nProcs()+1, 0);
-        for (label i = Pstream::myProcNo()+1; i < compactOffsets.size(); i++)
-        {
-            compactOffsets[i] = globalToCompact.size();
-        }
-
-        return autoPtr<globalIndex>::New(std::move(compactOffsets));
-    }
-
-
-
-    // Initial global region numbers
-    const globalIndex globalRegions(mesh().nFaces());
-
-    // Minimise regions across connected cells (including parallel)
-    calcNonCompactRegionSplit
+    const label nLocalRegions = calcLocalRegionSplit
     (
-        globalRegions,
         blockedFace,
         explicitConnections,
         cellRegion
     );
 
-
-    // Now our cellRegion will have
-    // - non-local regions (i.e. originating from other processors)
-    // - non-compact locally originating regions
-    // so we'll need to compact
-
-    // 4a: count per originating processor the number of regions
-    labelList nOriginating(Pstream::nProcs(), 0);
+    if (!doGlobalRegions)
     {
-        labelHashSet haveRegion(mesh().nCells()/8);
-
-        forAll(cellRegion, celli)
-        {
-            label region = cellRegion[celli];
-
-            // Count originating processor. Use isLocal as efficiency since
-            // most cells are locally originating.
-            if (globalRegions.isLocal(region))
-            {
-                if (haveRegion.insert(region))
-                {
-                    nOriginating[Pstream::myProcNo()]++;
-                }
-            }
-            else
-            {
-                label proci = globalRegions.whichProcID(region);
-                if (haveRegion.insert(region))
-                {
-                    nOriginating[proci]++;
-                }
-            }
-        }
+        return autoPtr<globalIndex>::New(nLocalRegions);
     }
 
-    if (debug)
-    {
-        Pout<< "Counted " << nOriginating[Pstream::myProcNo()]
-            << " local regions." << endl;
-    }
-
-
-    // Global numbering for compacted local regions
-    autoPtr<globalIndex> globalCompactPtr
-    (
-        new globalIndex(nOriginating[Pstream::myProcNo()])
-    );
-    const globalIndex& globalCompact = globalCompactPtr();
-
-
-    // 4b: renumber
-    // Renumber into compact indices. Note that since we've already made
-    // all regions global we now need a Map to store the compacting information
-    // instead of a labelList - otherwise we could have used a straight
-    // labelList.
-
-    // Local compaction map
-    Map<label> globalToCompact(2*nOriginating[Pstream::myProcNo()]);
-    // Remote regions we want the compact number for
-    List<labelHashSet> nonLocal(Pstream::nProcs());
-    forAll(nonLocal, proci)
-    {
-        if (proci != Pstream::myProcNo())
-        {
-            nonLocal[proci].resize(2*nOriginating[proci]);
-        }
-    }
-
-    forAll(cellRegion, celli)
-    {
-        label region = cellRegion[celli];
-        if (globalRegions.isLocal(region))
-        {
-            // Insert new compact region (if not yet present)
-            globalToCompact.insert
-            (
-                region,
-                globalCompact.toGlobal(globalToCompact.size())
-            );
-        }
-        else
-        {
-            nonLocal[globalRegions.whichProcID(region)].insert(region);
-        }
-    }
-
-
-    // Now we have all the local regions compacted. Now we need to get the
-    // non-local ones from the processors to whom they are local.
-    // Convert the nonLocal (labelHashSets) to labelLists.
-
-    labelListList sendNonLocal(Pstream::nProcs());
-    forAll(sendNonLocal, proci)
-    {
-        sendNonLocal[proci] = nonLocal[proci].toc();
-    }
-
-    if (debug)
-    {
-        forAll(sendNonLocal, proci)
-        {
-            Pout<< "    from processor " << proci
-                << " want " << sendNonLocal[proci].size()
-                << " region numbers."
-                << endl;
-        }
-        Pout<< endl;
-    }
-
-
-    // Get the wanted region labels into recvNonLocal
-    labelListList recvNonLocal;
-    Pstream::exchange<labelList, label>(sendNonLocal, recvNonLocal);
-
-    // Now we have the wanted compact region labels that proci wants in
-    // recvNonLocal[proci]. Construct corresponding list of compact
-    // region labels to send back.
-
-    labelListList sendWantedLocal(Pstream::nProcs());
-    forAll(recvNonLocal, proci)
-    {
-        const labelList& nonLocal = recvNonLocal[proci];
-        sendWantedLocal[proci].setSize(nonLocal.size());
-
-        forAll(nonLocal, i)
-        {
-            sendWantedLocal[proci][i] = globalToCompact[nonLocal[i]];
-        }
-    }
-
-
-    // Send back (into recvNonLocal)
-    recvNonLocal.clear();
-    Pstream::exchange<labelList, label>(sendWantedLocal, recvNonLocal);
-    sendWantedLocal.clear();
-
-    // Now recvNonLocal contains for every element in setNonLocal the
-    // corresponding compact number. Insert these into the local compaction
-    // map.
-
-    forAll(recvNonLocal, proci)
-    {
-        const labelList& wantedRegions = sendNonLocal[proci];
-        const labelList& compactRegions = recvNonLocal[proci];
-
-        forAll(wantedRegions, i)
-        {
-            globalToCompact.insert(wantedRegions[i], compactRegions[i]);
-        }
-    }
-
-    // Finally renumber the regions
-    forAll(cellRegion, celli)
-    {
-        cellRegion[celli] = globalToCompact[cellRegion[celli]];
-    }
-
-    return globalCompactPtr;
+    return reduceRegions(nLocalRegions, blockedFace, cellRegion);
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::regionSplit::regionSplit(const polyMesh& mesh, const bool doGlobalRegions)
+Foam::regionSplit::regionSplit
+(
+    const polyMesh& mesh,
+    const bool doGlobalRegions
+)
 :
     MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
     labelList(mesh.nCells(), -1)
@@ -530,4 +482,354 @@ Foam::regionSplit::regionSplit
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::reduceRegions
+(
+    const label numLocalRegions,
+    const boolList& blockedFace,
+
+    labelList& cellRegion
+) const
+{
+    clockValue timing(debug);
+
+    if (cellRegion.size() != mesh().nCells())
+    {
+        FatalErrorInFunction
+            << "The cellRegion size " << cellRegion.size()
+            << " is not equal to the of number of cells "
+            << mesh().nCells() << endl
+            << abort(FatalError);
+    }
+
+
+    // (numLocalRegions < 0) to signal that region information should be
+    // determined ourselves. This is not really efficient, but can be useful
+
+    const label nLocalRegions =
+    (
+        numLocalRegions < 0
+      ? labelHashSet(cellRegion).size()
+      : numLocalRegions
+    );
+
+
+    // Preliminary global region numbers
+    const globalIndex globalRegions(nLocalRegions);
+
+
+    // Lookup table of local region to global region.
+    // Initially an identity mapping of the uncombined global values.
+
+    Map<label> localToGlobal(2*nLocalRegions);
+    for (const label regioni : cellRegion)
+    {
+        localToGlobal.insert(regioni, globalRegions.toGlobal(regioni));
+    }
+
+    // To update the localToGlobal mapping during traversal of the boundaries
+    // and later when finalizing things.
+    Map<label> updateLookup(2*nLocalRegions);
+
+
+    // Note that we use two separate maps during the process.
+    // The localToGlobal is used to map the local to global regions.
+    // Merging across processors will normally make this a many->few mapping.
+    // However, we may need to walk up and down processor boundaries several
+    // times before all the information propagates through.
+    // During these traversals, it will normally be more efficient to just
+    // update the mapping without updating the cellRegion immediately.
+    // Only after everything is finalized do we renumber all of the cell
+    // regions.
+
+
+    // Merge global regions
+    // ~~~~~~~~~~~~~~~~~~~~
+    // Regions across non-blocked proc patches get merged.
+    // This will set merged global regions to be the min of both.
+    // (this will create gaps in the global region list so they will get
+    // merged later on)
+
+    const polyBoundaryMesh& patches = mesh().boundaryMesh();
+
+    // Buffer for swapping boundary information
+    labelList nbrRegion(mesh().nFaces()-mesh().nInternalFaces());
+
+    do
+    {
+        if (debug)
+        {
+            Pout<< nl << "-- Starting Iteration --" << endl;
+        }
+
+        updateLookup.clear();
+        nbrRegion = UNASSIGNED;
+
+        // Region information to send
+        for (const polyPatch& pp : patches)
+        {
+            if (pp.coupled())
+            {
+                const labelUList& faceCells = pp.faceCells();
+                SubList<label> patchNbrRegion
+                (
+                    nbrRegion,
+                    pp.size(),
+                    pp.start()-mesh().nInternalFaces()
+                );
+
+                forAll(faceCells, patchFacei)
+                {
+                    const label meshFacei = pp.start()+patchFacei;
+
+                    if (!blockedFace[meshFacei])
+                    {
+                        // Send the most currently updated region Id
+                        const label orig = cellRegion[faceCells[patchFacei]];
+
+                        patchNbrRegion[patchFacei] = localToGlobal[orig];
+                    }
+                }
+            }
+        }
+        syncTools::swapBoundaryFaceList(mesh(), nbrRegion);
+
+        // Receive and reduce region information
+        for (const polyPatch& pp : patches)
+        {
+            if (pp.coupled())
+            {
+                const labelUList& faceCells = pp.faceCells();
+                SubList<label> patchNbrRegion
+                (
+                    nbrRegion,
+                    pp.size(),
+                    pp.start()-mesh().nInternalFaces()
+                );
+
+                forAll(faceCells, patchFacei)
+                {
+                    const label meshFacei = pp.start()+patchFacei;
+
+                    if (!blockedFace[meshFacei])
+                    {
+                        // Reduction by retaining the min region id.
+
+                        const label orig = cellRegion[faceCells[patchFacei]];
+
+                        const label sent = localToGlobal[orig];
+                        const label recv = patchNbrRegion[patchFacei];
+
+                        // Record the minimum value seen
+                        if (recv < sent)
+                        {
+                            auto fnd = updateLookup.find(sent);
+                            if (!fnd.found())
+                            {
+                                updateLookup.insert(sent, recv);
+                            }
+                            else if (recv < *fnd)
+                            {
+                                *fnd = recv;
+                            }
+                        }
+                    }
+                }
+            }
+        }
+
+
+        // Note: by always using the minimum region number across the
+        // processor faces, we effect a consolidation of connected regions
+        // and converge to a unique number for each distinct region.
+
+
+        // Update localToGlobal according to newly exchanged information
+
+        inplaceMapValue(updateLookup, localToGlobal);
+
+        if (debug & 2)
+        {
+            labelList keys(localToGlobal.sortedToc());
+            labelList vals(keys.size());
+            forAll(keys, i)
+            {
+                vals[i] = localToGlobal[keys[i]];
+            }
+
+            Pout<< "Updated local regions:" << nl
+                << "old: " << flatOutput(keys) << nl
+                << "new: " << flatOutput(vals) << endl;
+        }
+        else if (debug)
+        {
+            Pout<< "Updated " << localToGlobal.size()
+                << " local regions" << endl;
+        }
+
+        // Continue until there are no further changes
+    }
+    while (returnReduce(!updateLookup.empty(), orOp<bool>()));
+
+
+    //
+    // We will need compact numbers locally and non-locally
+    //
+
+    // Determine the local compact numbering
+    label nCompact = 0;
+    {
+        labelHashSet localRegion(2*localToGlobal.size());
+
+        forAllConstIters(localToGlobal, iter)
+        {
+            const label regioni = iter.object();
+
+            if (globalRegions.isLocal(regioni))
+            {
+                localRegion.insert(regioni);
+            }
+        }
+
+        nCompact = localRegion.size();
+    }
+
+
+    // The new global numbering using compacted local regions
+    auto globalCompactPtr = autoPtr<globalIndex>::New(nCompact);
+    const auto& globalCompact = *globalCompactPtr;
+
+
+    // Determine the following:
+    // - the local compact regions (store as updateLookup)
+    // - the non-local regions, ordered according to the processor on which
+    //   they are local.
+
+
+    // The local compaction map (updated local to compact local numbering)
+    updateLookup.clear();
+
+    labelListList sendNonLocal(Pstream::nProcs());
+
+    {
+        List<labelHashSet> nonLocal(Pstream::nProcs(), labelHashSet(0));
+
+        // Use estimate of sizing for non-local regions
+        forAll(nonLocal, proci)
+        {
+            if (proci != Pstream::myProcNo())
+            {
+                nonLocal[proci].resize
+                (
+                    2*((nLocalRegions-nCompact)/Pstream::nProcs())
+                );
+            }
+        }
+
+
+        forAllConstIters(localToGlobal, iter)
+        {
+            const label regioni = iter.object();
+
+            if (globalRegions.isLocal(regioni))
+            {
+                updateLookup.insert
+                (
+                    regioni,
+                    globalCompact.toGlobal(updateLookup.size())
+                );
+            }
+            else
+            {
+                nonLocal[globalRegions.whichProcID(regioni)].insert(regioni);
+            }
+        }
+
+        if (debug)
+        {
+            Pout<< " per processor nonLocal regions: "
+                << flatOutput(containerSizes(nonLocal)) << endl;
+        }
+
+
+        // Convert to label list
+        forAll(sendNonLocal, proci)
+        {
+            sendNonLocal[proci] = nonLocal[proci].toc();
+        }
+    }
+
+
+    // Get the wanted region labels into recvNonLocal
+    labelListList recvNonLocal;
+    Pstream::exchange<labelList, label>
+    (
+        sendNonLocal,
+        recvNonLocal
+    );
+
+
+    // The recvNonLocal[proci] region labels are what proci requires.
+    // Transcribe into their compacted number.
+
+    {
+        labelListList sendLocal(std::move(recvNonLocal));
+
+        for (labelList& send : sendLocal)
+        {
+            for (label& regioni : send)
+            {
+                regioni = updateLookup[regioni];
+            }
+        }
+
+        // Send back (into recvNonLocal)
+        Pstream::exchange<labelList, label>
+        (
+            sendLocal,
+            recvNonLocal
+        );
+    }
+
+
+    // Now recvNonLocal and sendNonLocal contain matched pairs with
+    // sendNonLocal being the non-compact region and recvNonLocal being
+    // the compact region.
+    //
+    // Insert these into the local compaction map.
+
+    forAll(recvNonLocal, proci)
+    {
+        const labelList& send = sendNonLocal[proci];
+        const labelList& recv = recvNonLocal[proci];
+
+        forAll(send, i)
+        {
+            updateLookup.insert(send[i], recv[i]);
+        }
+    }
+
+
+    // Now renumber the localToGlobal to use the final compact global values
+    inplaceMapValue(updateLookup, localToGlobal);
+
+
+    // Can now finally use localToGlobal to renumber cellRegion
+
+    forAll(cellRegion, celli)
+    {
+        cellRegion[celli] = localToGlobal[cellRegion[celli]];
+    }
+
+    if (debug)
+    {
+        Info<<"regionSplit::reduceRegions = "
+            << double(timing.elapsed()) << "s\n";
+    }
+
+    return globalCompactPtr;
+}
+
+
 // ************************************************************************* //
diff --git a/src/meshTools/regionSplit/regionSplit.H b/src/meshTools/regionSplit/regionSplit.H
index 1051491f19..494e69cbf5 100644
--- a/src/meshTools/regionSplit/regionSplit.H
+++ b/src/meshTools/regionSplit/regionSplit.H
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -103,12 +103,11 @@ Description
 
     \endverbatim
 
-
     Can optionally keep all regions local to the processor.
 
-    Note: does not walk across cyclicAMI/cyclicACMI - since not 'coupled()'
-          at the patch level.
-
+Note
+    does not walk across cyclicAMI/cyclicACMI - since these are not
+    \c coupled() at the patch level.
 
 SourceFiles
     regionSplit.C
@@ -128,6 +127,8 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declarations
+
 class polyMesh;
 
 /*---------------------------------------------------------------------------*\
@@ -141,21 +142,52 @@ class regionSplit
 {
     // Private data
 
-        mutable autoPtr<globalIndex> globalNumberingPtr_;
+        autoPtr<globalIndex> globalNumberingPtr_;
+
+        //- Temporary list of cells that have changed
+        mutable DynamicList<label> changedCells_;
+
+        //- Temporary list of faces that have changed
+        mutable DynamicList<label> changedFaces_;
 
 
     // Private Member Functions
 
-        //- Calculate region split in non-compact (global) numbering.
-        void calcNonCompactRegionSplit
+        //- Update faceRegion data between (non-processor) coupled faces.
+        void updateFacePair
+        (
+            const label face0,
+            const label face1,
+            labelList& faceRegion,
+            DynamicList<label>& facesChanged
+        ) const;
+
+        //- Given a seed cell label, fill cellRegion/faceRegion with markValue
+        //- for contiguous region around it
+        void fillSeedMask
+        (
+            const List<labelPair>& explicitConnections,
+            const label seedCellID,
+            const label markValue,
+            labelList& cellRegion,
+            labelList& faceRegion
+        ) const;
+
+
+        //- Calculate the local region split.
+        //  \return number of processor-local regions,
+        //      without consolidation between procesors
+        label calcLocalRegionSplit
         (
-            const globalIndex& globalFaces,
             const boolList& blockedFace,
             const List<labelPair>& explicitConnections,
             labelList& cellRegion
         ) const;
 
-        //- Calculate global region split. Return globalIndex.
+
+        //- Calculate the local region split.
+        //  \return number of processor-local regions,
+        //      without consolidation between procesors
         autoPtr<globalIndex> calcRegionSplit
         (
             const bool doGlobalRegions,
@@ -176,27 +208,27 @@ public:
         //- Construct from mesh
         regionSplit
         (
-            const polyMesh&,
+            const polyMesh& mesh,
             const bool doGlobalRegions = Pstream::parRun()
         );
 
         //- Construct from mesh and whether face is blocked
-        //  \note blockedFace has to be consistent across coupled faces!
+        //  \note blockedFace must be consistent across coupled faces!
         regionSplit
         (
-            const polyMesh&,
+            const polyMesh& mesh,
             const boolList& blockedFace,
             const bool doGlobalRegions = Pstream::parRun()
         );
 
-        //- Construct from mesh and whether face is blocked.
-        //  Additional explicit connections between normal boundary faces.
-        //  \note blockedFace has to be consistent across coupled faces!
+        //- Construct from mesh and whether face is blocked, with additional
+        //- explicit connections between normal boundary faces.
+        //  \note blockedFace must be consistent across coupled faces!
         regionSplit
         (
-            const polyMesh&,
+            const polyMesh& mesh,
             const boolList& blockedFace,
-            const List<labelPair>&,
+            const List<labelPair>& explicitConnections,
             const bool doGlobalRegions = Pstream::parRun()
         );
 
@@ -220,6 +252,19 @@ public:
         {
             return globalNumbering().size();
         }
+
+
+        //- Manually consolidate the regions globally by swapping information
+        //  between processor domains and reducing the regions accordingly.
+        //
+        //  \return number of local regions after reduction.
+        autoPtr<globalIndex> reduceRegions
+        (
+            const label numLocalRegions,
+            const boolList& blockedFace,
+            labelList& cellRegion
+        ) const;
+
 };