diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index a3629c8e13..2dbac062ae 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -623,336 +623,349 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
 }
 
 
-// Determine for multi-processor regions the lowest numbered cell on the lowest
-// numbered processor.
-void Foam::meshRefinement::getCoupledRegionMaster
-(
-    const globalIndex& globalCells,
-    const boolList& blockedFace,
-    const regionSplit& globalRegion,
-    Map<label>& regionToMaster
-) const
-{
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (isA<processorPolyPatch>(pp))
-        {
-            forAll(pp, i)
-            {
-                label faceI = pp.start()+i;
-
-                if (!blockedFace[faceI])
-                {
-                    // Only if there is a connection to the neighbour
-                    // will there be a multi-domain region; if not through
-                    // this face then through another.
-
-                    label cellI = mesh_.faceOwner()[faceI];
-                    label globalCellI = globalCells.toGlobal(cellI);
-
-                    Map<label>::iterator iter =
-                        regionToMaster.find(globalRegion[cellI]);
-
-                    if (iter != regionToMaster.end())
-                    {
-                        label master = iter();
-                        iter() = min(master, globalCellI);
-                    }
-                    else
-                    {
-                        regionToMaster.insert
-                        (
-                            globalRegion[cellI],
-                            globalCellI
-                        );
-                    }
-                }
-            }
-        }
-    }
-
-    // Do reduction
-    Pstream::mapCombineGather(regionToMaster, minEqOp<label>());
-    Pstream::mapCombineScatter(regionToMaster);
-}
-
-
-void Foam::meshRefinement::calcLocalRegions
-(
-    const globalIndex& globalCells,
-    const labelList& globalRegion,
-    const Map<label>& coupledRegionToMaster,
-    const scalarField& cellWeights,
-
-    Map<label>& globalToLocalRegion,
-    pointField& localPoints,
-    scalarField& localWeights
-) const
-{
-    globalToLocalRegion.resize(globalRegion.size());
-    DynamicList<point> localCc(globalRegion.size()/2);
-    DynamicList<scalar> localWts(globalRegion.size()/2);
-
-    forAll(globalRegion, cellI)
-    {
-        Map<label>::const_iterator fndMaster =
-            coupledRegionToMaster.find(globalRegion[cellI]);
-
-        if (fndMaster != coupledRegionToMaster.end())
-        {
-            // Multi-processor region.
-            if (globalCells.toGlobal(cellI) == fndMaster())
-            {
-                // I am master. Allocate region for me.
-                globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
-                localCc.append(mesh_.cellCentres()[cellI]);
-                localWts.append(cellWeights[cellI]);
-            }
-        }
-        else
-        {
-            // Single processor region.
-            if (globalToLocalRegion.insert(globalRegion[cellI], localCc.size()))
-            {
-                localCc.append(mesh_.cellCentres()[cellI]);
-                localWts.append(cellWeights[cellI]);
-            }
-        }
-    }
-
-    localPoints.transfer(localCc);
-    localWeights.transfer(localWts);
-
-    if (localPoints.size() != globalToLocalRegion.size())
-    {
-        FatalErrorIn("calcLocalRegions(..)")
-            << "localPoints:" << localPoints.size()
-            << " globalToLocalRegion:" << globalToLocalRegion.size()
-            << exit(FatalError);
-    }
-}
-
-
-Foam::label Foam::meshRefinement::getShiftedRegion
-(
-    const globalIndex& indexer,
-    const Map<label>& globalToLocalRegion,
-    const Map<label>& coupledRegionToShifted,
-    const label globalRegion
-)
-{
-    Map<label>::const_iterator iter =
-        globalToLocalRegion.find(globalRegion);
-
-    if (iter != globalToLocalRegion.end())
-    {
-        // Region is 'owned' locally. Convert local region index into global.
-        return indexer.toGlobal(iter());
-    }
-    else
-    {
-        return coupledRegionToShifted[globalRegion];
-    }
-}
-
-
-// Add if not yet present
-void Foam::meshRefinement::addUnique(const label elem, labelList& lst)
-{
-    if (findIndex(lst, elem) == -1)
-    {
-        label sz = lst.size();
-        lst.setSize(sz+1);
-        lst[sz] = elem;
-    }
-}
-
-
-void Foam::meshRefinement::calcRegionRegions
-(
-    const labelList& globalRegion,
-    const Map<label>& globalToLocalRegion,
-    const Map<label>& coupledRegionToMaster,
-    labelListList& regionRegions
-) const
-{
-    // Global region indexing since we now know the shifted regions.
-    globalIndex shiftIndexer(globalToLocalRegion.size());
-
-    // Redo the coupledRegionToMaster to be in shifted region indexing.
-    Map<label> coupledRegionToShifted(coupledRegionToMaster.size());
-    forAllConstIter(Map<label>, coupledRegionToMaster, iter)
-    {
-        label region = iter.key();
-
-        Map<label>::const_iterator fndRegion = globalToLocalRegion.find(region);
-
-        if (fndRegion != globalToLocalRegion.end())
-        {
-            // A local cell is master of this region. Get its shifted region.
-            coupledRegionToShifted.insert
-            (
-                region,
-                shiftIndexer.toGlobal(fndRegion())
-            );
-        }
-        Pstream::mapCombineGather(coupledRegionToShifted, minEqOp<label>());
-        Pstream::mapCombineScatter(coupledRegionToShifted);
-    }
-
-
-    // Determine region-region connectivity.
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // This is for all my regions (so my local ones or the ones I am master
-    // of) the neighbouring region indices.
-
-
-    // Transfer lists.
-    PtrList<HashSet<edge, Hash<edge> > > regionConnectivity(Pstream::nProcs());
-    forAll(regionConnectivity, procI)
-    {
-        if (procI != Pstream::myProcNo())
-        {
-            regionConnectivity.set
-            (
-                procI,
-                new HashSet<edge, Hash<edge> >
-                (
-                    coupledRegionToShifted.size()
-                  / Pstream::nProcs()
-                )
-            );
-        }
-    }
-
-
-    // Connectivity. For all my local regions the connected regions.
-    regionRegions.setSize(globalToLocalRegion.size());
-
-    // Add all local connectivity to regionRegions; add all non-local
-    // to the transferlists.
-    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
-    {
-        label ownRegion = globalRegion[mesh_.faceOwner()[faceI]];
-        label neiRegion = globalRegion[mesh_.faceNeighbour()[faceI]];
-
-        if (ownRegion != neiRegion)
-        {
-            label shiftOwnRegion = getShiftedRegion
-            (
-                shiftIndexer,
-                globalToLocalRegion,
-                coupledRegionToShifted,
-                ownRegion
-            );
-            label shiftNeiRegion = getShiftedRegion
-            (
-                shiftIndexer,
-                globalToLocalRegion,
-                coupledRegionToShifted,
-                neiRegion
-            );
-
-
-            // Connection between two regions. Send to owner of region.
-            // - is ownRegion 'owned' by me
-            // - is neiRegion 'owned' by me
-
-            if (shiftIndexer.isLocal(shiftOwnRegion))
-            {
-                label localI = shiftIndexer.toLocal(shiftOwnRegion);
-                addUnique(shiftNeiRegion, regionRegions[localI]);
-            }
-            else
-            {
-                label masterProc = shiftIndexer.whichProcID(shiftOwnRegion);
-                regionConnectivity[masterProc].insert
-                (
-                    edge(shiftOwnRegion, shiftNeiRegion)
-                );
-            }
-
-            if (shiftIndexer.isLocal(shiftNeiRegion))
-            {
-                label localI = shiftIndexer.toLocal(shiftNeiRegion);
-                addUnique(shiftOwnRegion, regionRegions[localI]);
-            }
-            else
-            {
-                label masterProc = shiftIndexer.whichProcID(shiftNeiRegion);
-                regionConnectivity[masterProc].insert
-                (
-                    edge(shiftOwnRegion, shiftNeiRegion)
-                );
-            }
-        }
-    }
-
-
-    // Send
-    forAll(regionConnectivity, procI)
-    {
-        if (procI != Pstream::myProcNo())
-        {
-            OPstream str(Pstream::blocking, procI);
-            str << regionConnectivity[procI];
-        }
-    }
-    // Receive
-    forAll(regionConnectivity, procI)
-    {
-        if (procI != Pstream::myProcNo())
-        {
-            IPstream str(Pstream::blocking, procI);
-            str >> regionConnectivity[procI];
-        }
-    }
-
-    // Add to addressing.
-    forAll(regionConnectivity, procI)
-    {
-        if (procI != Pstream::myProcNo())
-        {
-            for
-            (
-                HashSet<edge, Hash<edge> >::const_iterator iter =
-                    regionConnectivity[procI].begin();
-                iter != regionConnectivity[procI].end();
-                ++iter
-            )
-            {
-                const edge& e = iter.key();
-
-                bool someLocal = false;
-                if (shiftIndexer.isLocal(e[0]))
-                {
-                    label localI = shiftIndexer.toLocal(e[0]);
-                    addUnique(e[1], regionRegions[localI]);
-                    someLocal = true;
-                }
-                if (shiftIndexer.isLocal(e[1]))
-                {
-                    label localI = shiftIndexer.toLocal(e[1]);
-                    addUnique(e[0], regionRegions[localI]);
-                    someLocal = true;
-                }
-
-                if (!someLocal)
-                {
-                    FatalErrorIn("calcRegionRegions(..)")
-                        << "Received from processor " << procI
-                        << " connection " << e
-                        << " where none of the elements is local to me."
-                        << abort(FatalError);
-                }
-            }
-        }
-    }
-}
+//// Determine for multi-processor regions the lowest numbered cell on the
+//// lowest numbered processor.
+//void Foam::meshRefinement::getCoupledRegionMaster
+//(
+//    const globalIndex& globalCells,
+//    const boolList& blockedFace,
+//    const regionSplit& globalRegion,
+//    Map<label>& regionToMaster
+//) const
+//{
+//    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+//
+//    forAll(patches, patchI)
+//    {
+//        const polyPatch& pp = patches[patchI];
+//
+//        if (isA<processorPolyPatch>(pp))
+//        {
+//            forAll(pp, i)
+//            {
+//                label faceI = pp.start()+i;
+//
+//                if (!blockedFace[faceI])
+//                {
+//                    // Only if there is a connection to the neighbour
+//                    // will there be a multi-domain region; if not through
+//                    // this face then through another.
+//
+//                    label cellI = mesh_.faceOwner()[faceI];
+//                    label globalCellI = globalCells.toGlobal(cellI);
+//
+//                    Map<label>::iterator iter =
+//                        regionToMaster.find(globalRegion[cellI]);
+//
+//                    if (iter != regionToMaster.end())
+//                    {
+//                        label master = iter();
+//                        iter() = min(master, globalCellI);
+//                    }
+//                    else
+//                    {
+//                        regionToMaster.insert
+//                        (
+//                            globalRegion[cellI],
+//                            globalCellI
+//                        );
+//                    }
+//                }
+//            }
+//        }
+//    }
+//
+//    // Do reduction
+//    Pstream::mapCombineGather(regionToMaster, minEqOp<label>());
+//    Pstream::mapCombineScatter(regionToMaster);
+//}
+//
+//
+//void Foam::meshRefinement::calcLocalRegions
+//(
+//    const globalIndex& globalCells,
+//    const labelList& globalRegion,
+//    const Map<label>& coupledRegionToMaster,
+//    const scalarField& cellWeights,
+//
+//    Map<label>& globalToLocalRegion,
+//    pointField& localPoints,
+//    scalarField& localWeights
+//) const
+//{
+//    globalToLocalRegion.resize(globalRegion.size());
+//    DynamicList<point> localCc(globalRegion.size()/2);
+//    DynamicList<scalar> localWts(globalRegion.size()/2);
+//
+//    forAll(globalRegion, cellI)
+//    {
+//        Map<label>::const_iterator fndMaster =
+//            coupledRegionToMaster.find(globalRegion[cellI]);
+//
+//        if (fndMaster != coupledRegionToMaster.end())
+//        {
+//            // Multi-processor region.
+//            if (globalCells.toGlobal(cellI) == fndMaster())
+//            {
+//                // I am master. Allocate region for me.
+//                globalToLocalRegion.insert
+//                (
+//                    globalRegion[cellI],
+//                    localCc.size()
+//                );
+//                localCc.append(mesh_.cellCentres()[cellI]);
+//                localWts.append(cellWeights[cellI]);
+//            }
+//        }
+//        else
+//        {
+//            // Single processor region.
+//            if
+//            (
+//                globalToLocalRegion.insert
+//                (
+//                    globalRegion[cellI],
+//                    localCc.size()
+//                )
+//            )
+//            {
+//                localCc.append(mesh_.cellCentres()[cellI]);
+//                localWts.append(cellWeights[cellI]);
+//            }
+//        }
+//    }
+//
+//    localPoints.transfer(localCc);
+//    localWeights.transfer(localWts);
+//
+//    if (localPoints.size() != globalToLocalRegion.size())
+//    {
+//        FatalErrorIn("calcLocalRegions(..)")
+//            << "localPoints:" << localPoints.size()
+//            << " globalToLocalRegion:" << globalToLocalRegion.size()
+//            << exit(FatalError);
+//    }
+//}
+//
+//
+//Foam::label Foam::meshRefinement::getShiftedRegion
+//(
+//    const globalIndex& indexer,
+//    const Map<label>& globalToLocalRegion,
+//    const Map<label>& coupledRegionToShifted,
+//    const label globalRegion
+//)
+//{
+//    Map<label>::const_iterator iter =
+//        globalToLocalRegion.find(globalRegion);
+//
+//    if (iter != globalToLocalRegion.end())
+//    {
+//        // Region is 'owned' locally. Convert local region index into global.
+//        return indexer.toGlobal(iter());
+//    }
+//    else
+//    {
+//        return coupledRegionToShifted[globalRegion];
+//    }
+//}
+//
+//
+//// Add if not yet present
+//void Foam::meshRefinement::addUnique(const label elem, labelList& lst)
+//{
+//    if (findIndex(lst, elem) == -1)
+//    {
+//        label sz = lst.size();
+//        lst.setSize(sz+1);
+//        lst[sz] = elem;
+//    }
+//}
+//
+//
+//void Foam::meshRefinement::calcRegionRegions
+//(
+//    const labelList& globalRegion,
+//    const Map<label>& globalToLocalRegion,
+//    const Map<label>& coupledRegionToMaster,
+//    labelListList& regionRegions
+//) const
+//{
+//    // Global region indexing since we now know the shifted regions.
+//    globalIndex shiftIndexer(globalToLocalRegion.size());
+//
+//    // Redo the coupledRegionToMaster to be in shifted region indexing.
+//    Map<label> coupledRegionToShifted(coupledRegionToMaster.size());
+//    forAllConstIter(Map<label>, coupledRegionToMaster, iter)
+//    {
+//        label region = iter.key();
+//
+//        Map<label>::const_iterator fndRegion = globalToLocalRegion.find
+//        (region);
+//
+//        if (fndRegion != globalToLocalRegion.end())
+//        {
+//            // A local cell is master of this region. Get its shifted region.
+//            coupledRegionToShifted.insert
+//            (
+//                region,
+//                shiftIndexer.toGlobal(fndRegion())
+//            );
+//        }
+//        Pstream::mapCombineGather(coupledRegionToShifted, minEqOp<label>());
+//        Pstream::mapCombineScatter(coupledRegionToShifted);
+//    }
+//
+//
+//    // Determine region-region connectivity.
+//    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+//    // This is for all my regions (so my local ones or the ones I am master
+//    // of) the neighbouring region indices.
+//
+//
+//    // Transfer lists.
+//    PtrList<HashSet<edge, Hash<edge> > > regionConnectivity
+//    (Pstream::nProcs());
+//    forAll(regionConnectivity, procI)
+//    {
+//        if (procI != Pstream::myProcNo())
+//        {
+//            regionConnectivity.set
+//            (
+//                procI,
+//                new HashSet<edge, Hash<edge> >
+//                (
+//                    coupledRegionToShifted.size()
+//                  / Pstream::nProcs()
+//                )
+//            );
+//        }
+//    }
+//
+//
+//    // Connectivity. For all my local regions the connected regions.
+//    regionRegions.setSize(globalToLocalRegion.size());
+//
+//    // Add all local connectivity to regionRegions; add all non-local
+//    // to the transferlists.
+//    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+//    {
+//        label ownRegion = globalRegion[mesh_.faceOwner()[faceI]];
+//        label neiRegion = globalRegion[mesh_.faceNeighbour()[faceI]];
+//
+//        if (ownRegion != neiRegion)
+//        {
+//            label shiftOwnRegion = getShiftedRegion
+//            (
+//                shiftIndexer,
+//                globalToLocalRegion,
+//                coupledRegionToShifted,
+//                ownRegion
+//            );
+//            label shiftNeiRegion = getShiftedRegion
+//            (
+//                shiftIndexer,
+//                globalToLocalRegion,
+//                coupledRegionToShifted,
+//                neiRegion
+//            );
+//
+//
+//            // Connection between two regions. Send to owner of region.
+//            // - is ownRegion 'owned' by me
+//            // - is neiRegion 'owned' by me
+//
+//            if (shiftIndexer.isLocal(shiftOwnRegion))
+//            {
+//                label localI = shiftIndexer.toLocal(shiftOwnRegion);
+//                addUnique(shiftNeiRegion, regionRegions[localI]);
+//            }
+//            else
+//            {
+//                label masterProc = shiftIndexer.whichProcID(shiftOwnRegion);
+//                regionConnectivity[masterProc].insert
+//                (
+//                    edge(shiftOwnRegion, shiftNeiRegion)
+//                );
+//            }
+//
+//            if (shiftIndexer.isLocal(shiftNeiRegion))
+//            {
+//                label localI = shiftIndexer.toLocal(shiftNeiRegion);
+//                addUnique(shiftOwnRegion, regionRegions[localI]);
+//            }
+//            else
+//            {
+//                label masterProc = shiftIndexer.whichProcID(shiftNeiRegion);
+//                regionConnectivity[masterProc].insert
+//                (
+//                    edge(shiftOwnRegion, shiftNeiRegion)
+//                );
+//            }
+//        }
+//    }
+//
+//
+//    // Send
+//    forAll(regionConnectivity, procI)
+//    {
+//        if (procI != Pstream::myProcNo())
+//        {
+//            OPstream str(Pstream::blocking, procI);
+//            str << regionConnectivity[procI];
+//        }
+//    }
+//    // Receive
+//    forAll(regionConnectivity, procI)
+//    {
+//        if (procI != Pstream::myProcNo())
+//        {
+//            IPstream str(Pstream::blocking, procI);
+//            str >> regionConnectivity[procI];
+//        }
+//    }
+//
+//    // Add to addressing.
+//    forAll(regionConnectivity, procI)
+//    {
+//        if (procI != Pstream::myProcNo())
+//        {
+//            for
+//            (
+//                HashSet<edge, Hash<edge> >::const_iterator iter =
+//                    regionConnectivity[procI].begin();
+//                iter != regionConnectivity[procI].end();
+//                ++iter
+//            )
+//            {
+//                const edge& e = iter.key();
+//
+//                bool someLocal = false;
+//                if (shiftIndexer.isLocal(e[0]))
+//                {
+//                    label localI = shiftIndexer.toLocal(e[0]);
+//                    addUnique(e[1], regionRegions[localI]);
+//                    someLocal = true;
+//                }
+//                if (shiftIndexer.isLocal(e[1]))
+//                {
+//                    label localI = shiftIndexer.toLocal(e[1]);
+//                    addUnique(e[0], regionRegions[localI]);
+//                    someLocal = true;
+//                }
+//
+//                if (!someLocal)
+//                {
+//                    FatalErrorIn("calcRegionRegions(..)")
+//                        << "Received from processor " << procI
+//                        << " connection " << e
+//                        << " where none of the elements is local to me."
+//                        << abort(FatalError);
+//                }
+//            }
+//        }
+//    }
+//}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -1021,146 +1034,148 @@ Foam::label Foam::meshRefinement::countHits() const
 }
 
 
-// Determine distribution to move connected regions onto one processor.
-Foam::labelList Foam::meshRefinement::decomposeCombineRegions
-(
-    const scalarField& cellWeights,
-    const boolList& blockedFace,
-    const List<labelPair>& explicitConnections,
-    decompositionMethod& decomposer
-) const
-{
-    // Determine global regions, separated by blockedFaces
-    regionSplit globalRegion(mesh_, blockedFace, explicitConnections);
-
-    // Now globalRegion has global region per cell. Problem is that
-    // the region might span multiple domains so we want to get
-    // a "region master" per domain. Note that multi-processor
-    // regions can only occur on cells on coupled patches.
-    // Note: since the number of regions does not change by this the
-    // process can be seen as just a shift of a region onto a single
-    // processor.
-
-
-    // Global cell numbering engine
-    globalIndex globalCells(mesh_.nCells());
-
-
-    // Determine per coupled region the master cell (lowest numbered cell
-    // on lowest numbered processor)
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // (does not determine master for non-coupled=fully-local regions)
-
-    Map<label> coupledRegionToMaster(mesh_.nFaces()-mesh_.nInternalFaces());
-    getCoupledRegionMaster
-    (
-        globalCells,
-        blockedFace,
-        globalRegion,
-        coupledRegionToMaster
-    );
-
-    // Determine my regions
-    // ~~~~~~~~~~~~~~~~~~~~
-    // These are the fully local ones or the coupled ones of which I am master.
-
-    Map<label> globalToLocalRegion;
-    pointField localPoints;
-    scalarField localWeights;
-    calcLocalRegions
-    (
-        globalCells,
-        globalRegion,
-        coupledRegionToMaster,
-        cellWeights,
-
-        globalToLocalRegion,
-        localPoints,
-        localWeights
-    );
-
-
-
-    // Find distribution for regions
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    labelList regionDistribution;
-
-    if (isA<geomDecomp>(decomposer))
-    {
-        regionDistribution = decomposer.decompose(localPoints, localWeights);
-    }
-    else
-    {
-        labelListList regionRegions;
-        calcRegionRegions
-        (
-            globalRegion,
-            globalToLocalRegion,
-            coupledRegionToMaster,
-
-            regionRegions
-        );
-
-        regionDistribution = decomposer.decompose
-        (
-            regionRegions,
-            localPoints,
-            localWeights
-        );
-    }
-
-
-
-    // Convert region-based decomposition back to cell-based one
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Transfer destination processor back to all. Use global reduce for now.
-    Map<label> regionToDist(coupledRegionToMaster.size());
-    forAllConstIter(Map<label>, coupledRegionToMaster, iter)
-    {
-        label region = iter.key();
-
-        Map<label>::const_iterator regionFnd = globalToLocalRegion.find(region);
-
-        if (regionFnd != globalToLocalRegion.end())
-        {
-            // Master cell is local. Store my distribution.
-            regionToDist.insert(iter.key(), regionDistribution[regionFnd()]);
-        }
-        else
-        {
-            // Master cell is not on this processor. Make sure it is overridden.
-            regionToDist.insert(iter.key(), labelMax);
-        }
-    }
-    Pstream::mapCombineGather(regionToDist, minEqOp<label>());
-    Pstream::mapCombineScatter(regionToDist);
-
-
-    // Determine destination for all cells
-    labelList distribution(mesh_.nCells());
-
-    forAll(globalRegion, cellI)
-    {
-        Map<label>::const_iterator fndRegion =
-            regionToDist.find(globalRegion[cellI]);
-
-        if (fndRegion != regionToDist.end())
-        {
-            distribution[cellI] = fndRegion();
-        }
-        else
-        {
-            // region is local to the processor.
-            label localRegionI = globalToLocalRegion[globalRegion[cellI]];
-
-            distribution[cellI] = regionDistribution[localRegionI];
-        }
-    }
-
-    return distribution;
-}
+//// Determine distribution to move connected regions onto one processor.
+//Foam::labelList Foam::meshRefinement::decomposeCombineRegions
+//(
+//    const scalarField& cellWeights,
+//    const boolList& blockedFace,
+//    const List<labelPair>& explicitConnections,
+//    decompositionMethod& decomposer
+//) const
+//{
+//    // Determine global regions, separated by blockedFaces
+//    regionSplit globalRegion(mesh_, blockedFace, explicitConnections);
+//
+//    // Now globalRegion has global region per cell. Problem is that
+//    // the region might span multiple domains so we want to get
+//    // a "region master" per domain. Note that multi-processor
+//    // regions can only occur on cells on coupled patches.
+//    // Note: since the number of regions does not change by this the
+//    // process can be seen as just a shift of a region onto a single
+//    // processor.
+//
+//
+//    // Global cell numbering engine
+//    globalIndex globalCells(mesh_.nCells());
+//
+//
+//    // Determine per coupled region the master cell (lowest numbered cell
+//    // on lowest numbered processor)
+//    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+//    // (does not determine master for non-coupled=fully-local regions)
+//
+//    Map<label> coupledRegionToMaster(mesh_.nFaces()-mesh_.nInternalFaces());
+//    getCoupledRegionMaster
+//    (
+//        globalCells,
+//        blockedFace,
+//        globalRegion,
+//        coupledRegionToMaster
+//    );
+//
+//    // Determine my regions
+//    // ~~~~~~~~~~~~~~~~~~~~
+//    // These are the fully local ones or the coupled ones of which I am master
+//
+//    Map<label> globalToLocalRegion;
+//    pointField localPoints;
+//    scalarField localWeights;
+//    calcLocalRegions
+//    (
+//        globalCells,
+//        globalRegion,
+//        coupledRegionToMaster,
+//        cellWeights,
+//
+//        globalToLocalRegion,
+//        localPoints,
+//        localWeights
+//    );
+//
+//
+//
+//    // Find distribution for regions
+//    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+//
+//    labelList regionDistribution;
+//
+//    if (isA<geomDecomp>(decomposer))
+//    {
+//        regionDistribution = decomposer.decompose(localPoints, localWeights);
+//    }
+//    else
+//    {
+//        labelListList regionRegions;
+//        calcRegionRegions
+//        (
+//            globalRegion,
+//            globalToLocalRegion,
+//            coupledRegionToMaster,
+//
+//            regionRegions
+//        );
+//
+//        regionDistribution = decomposer.decompose
+//        (
+//            regionRegions,
+//            localPoints,
+//            localWeights
+//        );
+//    }
+//
+//
+//
+//    // Convert region-based decomposition back to cell-based one
+//    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+//
+//    // Transfer destination processor back to all. Use global reduce for now.
+//    Map<label> regionToDist(coupledRegionToMaster.size());
+//    forAllConstIter(Map<label>, coupledRegionToMaster, iter)
+//    {
+//        label region = iter.key();
+//
+//        Map<label>::const_iterator regionFnd = globalToLocalRegion.find
+//        (region);
+//
+//        if (regionFnd != globalToLocalRegion.end())
+//        {
+//            // Master cell is local. Store my distribution.
+//            regionToDist.insert(iter.key(), regionDistribution[regionFnd()]);
+//        }
+//        else
+//        {
+//            // Master cell is not on this processor. Make sure it is
+//            // overridden.
+//            regionToDist.insert(iter.key(), labelMax);
+//        }
+//    }
+//    Pstream::mapCombineGather(regionToDist, minEqOp<label>());
+//    Pstream::mapCombineScatter(regionToDist);
+//
+//
+//    // Determine destination for all cells
+//    labelList distribution(mesh_.nCells());
+//
+//    forAll(globalRegion, cellI)
+//    {
+//        Map<label>::const_iterator fndRegion =
+//            regionToDist.find(globalRegion[cellI]);
+//
+//        if (fndRegion != regionToDist.end())
+//        {
+//            distribution[cellI] = fndRegion();
+//        }
+//        else
+//        {
+//            // region is local to the processor.
+//            label localRegionI = globalToLocalRegion[globalRegion[cellI]];
+//
+//            distribution[cellI] = regionDistribution[localRegionI];
+//        }
+//    }
+//
+//    return distribution;
+//}
 
 
 Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
@@ -1184,22 +1199,38 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
         // Wanted distribution
         labelList distribution;
 
+
+        // Faces where owner and neighbour are not 'connected' so can
+        // go to different processors.
+        boolList blockedFace;
+        label nUnblocked = 0;
+
+        // Faces that move as block onto single processor
+        PtrList<labelList> specifiedProcessorFaces;
+        labelList specifiedProcessor;
+
+        // Constraints from decomposeParDict
+        decomposer.setConstraints
+        (
+            mesh_,
+            blockedFace,
+            specifiedProcessorFaces,
+            specifiedProcessor
+        );
+
+        // Pairs of baffles
+        List<labelPair> couples;
+
+
         if (keepZoneFaces || keepBaffles)
         {
-            // Faces where owner and neighbour are not 'connected' so can
-            // go to different processors.
-            boolList blockedFace(mesh_.nFaces(), true);
-            label nUnblocked = 0;
-            // Pairs of baffles
-            List<labelPair> couples;
-
             if (keepZoneFaces)
             {
                 // Determine decomposition to keep/move surface zones
                 // on one processor. The reason is that snapping will make these
                 // into baffles, move and convert them back so if they were
                 // proc boundaries after baffling&moving the points might be no
-                // longer snychronised so recoupling will fail. To prevent this
+                // longer synchronised so recoupling will fail. To prevent this
                 // keep owner&neighbour of such a surface zone on the same
                 // processor.
 
@@ -1255,6 +1286,32 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
                     << " zoned faces to keep together." << endl;
             }
 
+
+            // Extend un-blockedFaces with any cyclics
+            {
+                boolList separatedCoupledFace(mesh_.nFaces(), false);
+                selectSeparatedCoupledFaces(separatedCoupledFace);
+
+                label nSeparated = 0;
+                forAll(separatedCoupledFace, faceI)
+                {
+                    if (separatedCoupledFace[faceI])
+                    {
+                        if (blockedFace[faceI])
+                        {
+                            blockedFace[faceI] = false;
+                            nSeparated++;
+                        }
+                    }
+                }
+                reduce(nSeparated, sumOp<label>());
+                Info<< "Found " << nSeparated
+                    << " separated coupled faces to keep together." << endl;
+
+                nUnblocked += nSeparated;
+            }
+
+
             if (keepBaffles)
             {
                 // Get boundary baffles that need to stay together.
@@ -1272,51 +1329,61 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
                     << endl;
             }
 
-            if (nUnblocked > 0 || nCouples > 0)
-            {
-                Info<< "Applying special decomposition to keep baffles"
-                    << " and zoned faces together." << endl;
-
-                distribution = decomposeCombineRegions
-                (
-                    cellWeights,
-                    blockedFace,
-                    couples,
-                    decomposer
-                );
-
-                labelList nProcCells(distributor.countCells(distribution));
-                Pstream::listCombineGather(nProcCells, plusEqOp<label>());
-                Pstream::listCombineScatter(nProcCells);
-
-                Info<< "Calculated decomposition:" << endl;
-                forAll(nProcCells, procI)
-                {
-                    Info<< "    " << procI << '\t' << nProcCells[procI] << endl;
-                }
-                Info<< endl;
-            }
-            else
-            {
-                // Normal decomposition
-                distribution = decomposer.decompose
-                (
-                    mesh_,
-                    mesh_.cellCentres(),
-                    cellWeights
-                );
-            }
-        }
-        else
-        {
-            // Normal decomposition
-            distribution = decomposer.decompose
-            (
-                mesh_,
-                mesh_.cellCentres(),
-                cellWeights
-            );
+            //if (nUnblocked > 0 || nCouples > 0)
+            //{
+            //    Info<< "Applying special decomposition to keep baffles"
+            //        << " and zoned faces together." << endl;
+            //
+            //    distribution = decomposeCombineRegions
+            //    (
+            //        cellWeights,
+            //        blockedFace,
+            //        couples,
+            //        decomposer
+            //    );
+            //
+            //    labelList nProcCells(distributor.countCells(distribution));
+            //    Pstream::listCombineGather(nProcCells, plusEqOp<label>());
+            //    Pstream::listCombineScatter(nProcCells);
+            //
+            //    Info<< "Calculated decomposition:" << endl;
+            //    forAll(nProcCells, procI)
+            //    {
+            //        Info<< "    " << procI << '\t' << nProcCells[procI]
+            //            << endl;
+            //    }
+            //    Info<< endl;
+            //}
+            //else
+            //{
+            //    // Normal decomposition
+            //    distribution = decomposer.decompose
+            //    (
+            //        mesh_,
+            //        mesh_.cellCentres(),
+            //        cellWeights
+            //    );
+            //}
         }
+        //else
+        //{
+        //    // Normal decomposition
+        //    distribution = decomposer.decompose
+        //    (
+        //        mesh_,
+        //        cellWeights
+        //    );
+        //}
+
+        distribution = decomposer.decompose
+        (
+            mesh_,
+            cellWeights,
+            blockedFace,
+            specifiedProcessorFaces,
+            specifiedProcessor,
+            couples                 // explicit connections
+        );
 
         if (debug)
         {
@@ -1343,6 +1410,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
         mesh_.setInstance(timeName());
         setInstance(mesh_.facesInstance());
     }
+
     return map;
 }
 
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 1b5662db32..8d68415e22 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -193,54 +193,6 @@ private:
             const label exposedPatchI
         );
 
-        // For decomposeCombineRegions
-
-            //- Used in decomposeCombineRegions. Given global region per cell
-            //  determines master processor/cell for regions straddling
-            //  procboundaries.
-            void getCoupledRegionMaster
-            (
-                const globalIndex& globalCells,
-                const boolList& blockedFace,
-                const regionSplit& globalRegion,
-                Map<label>& regionToMaster
-            ) const;
-
-            //- Determine regions that are local to me or coupled ones that
-            //  are owned by me. Determine representative location.
-            void calcLocalRegions
-            (
-                const globalIndex& globalCells,
-                const labelList& globalRegion,
-                const Map<label>& coupledRegionToMaster,
-                const scalarField& cellWeights,
-
-                Map<label>& globalToLocalRegion,
-                pointField& localPoints,
-                scalarField& localWeights
-            ) const;
-
-            //- Convert region into global index.
-            static label getShiftedRegion
-            (
-                const globalIndex& indexer,
-                const Map<label>& globalToLocalRegion,
-                const Map<label>& coupledRegionToShifted,
-                const label globalRegion
-            );
-
-            //- helper: add element if not in list. Linear search.
-            static void addUnique(const label, labelList&);
-
-            //- Calculate region connectivity. Major communication.
-            void calcRegionRegions
-            (
-                const labelList& globalRegion,
-                const Map<label>& globalToLocalRegion,
-                const Map<label>& coupledRegionToMaster,
-                labelListList& regionRegions
-            ) const;
-
 
         // Refinement candidate selection
 
@@ -668,19 +620,6 @@ public:
             //- Count number of intersections (local)
             label countHits() const;
 
-            //- Helper function to get decomposition such that all connected
-            //  regions get moved onto one processor. Used to prevent baffles
-            //  straddling processor boundaries. explicitConnections is to
-            //  keep pairs of non-coupled boundary faces together
-            //  (e.g. to keep baffles together)
-            labelList decomposeCombineRegions
-            (
-                const scalarField& cellWeights,
-                const boolList& blockedFace,
-                const List<labelPair>& explicitConnections,
-                decompositionMethod&
-            ) const;
-
             //- Redecompose according to cell count
             //  keepZoneFaces : find all faceZones from zoned surfaces and keep
             //                  owner and neighbour together
diff --git a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
index 9ae4704239..a04edc3149 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
+++ b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
@@ -745,133 +745,64 @@ void Foam::decompositionMethod::calcCellCells
 Foam::labelList Foam::decompositionMethod::decompose
 (
     const polyMesh& mesh,
-    const scalarField& cellWeights
+    const scalarField& cellWeights,
+
+    //- Whether owner and neighbour should be on same processor
+    //  (takes priority over explicitConnections)
+    const boolList& blockedFace,
+
+    //- Whether whole sets of faces (and point neighbours) need to be kept
+    //  on single processor
+    const PtrList<labelList>& specifiedProcessorFaces,
+    const labelList& specifiedProcessor,
+
+    //- Additional connections between boundary faces
+    const List<labelPair>& explicitConnections
 )
 {
-    labelHashSet sameProcFaces;
+    // Any weights specified?
+    label nWeights = returnReduce(cellWeights.size(), sumOp<label>());
 
-    if (decompositionDict_.found("preservePatches"))
+    // Any processor sets?
+    label nProcSets = 0;
+    forAll(specifiedProcessorFaces, setI)
     {
-        wordList pNames(decompositionDict_.lookup("preservePatches"));
-
-        Info<< nl
-            << "Keeping owner of faces in patches " << pNames
-            << " on same processor. This only makes sense for cyclics." << endl;
-
-        const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-        forAll(pNames, i)
-        {
-            const label patchI = patches.findPatchID(pNames[i]);
-
-            if (patchI == -1)
-            {
-                FatalErrorIn("decompositionMethod::decompose(const polyMesh&)")
-                    << "Unknown preservePatch " << pNames[i]
-                    << endl << "Valid patches are " << patches.names()
-                    << exit(FatalError);
-            }
-
-            const polyPatch& pp = patches[patchI];
-
-            forAll(pp, i)
-            {
-                sameProcFaces.insert(pp.start() + i);
-            }
-        }
+        nProcSets += specifiedProcessorFaces[setI].size();
     }
-    if (decompositionDict_.found("preserveFaceZones"))
-    {
-        wordList zNames(decompositionDict_.lookup("preserveFaceZones"));
+    reduce(nProcSets, sumOp<label>());
 
-        Info<< nl
-            << "Keeping owner and neighbour of faces in zones " << zNames
-            << " on same processor" << endl;
-
-        const faceZoneMesh& fZones = mesh.faceZones();
-
-        forAll(zNames, i)
-        {
-            label zoneI = fZones.findZoneID(zNames[i]);
-
-            if (zoneI == -1)
-            {
-                FatalErrorIn("decompositionMethod::decompose(const polyMesh&)")
-                    << "Unknown preserveFaceZone " << zNames[i]
-                    << endl << "Valid faceZones are " << fZones.names()
-                    << exit(FatalError);
-            }
-
-            const faceZone& fz = fZones[zoneI];
-
-            forAll(fz, i)
-            {
-                sameProcFaces.insert(fz[i]);
-            }
-        }
-    }
-
-
-    // Specified processor for group of cells connected to faces
-
-    //- Sets of faces to move together
-    PtrList<labelList> specifiedProcessorFaces;
-    //- Destination processor
-    labelList specifiedProcessor;
-
-    if (decompositionDict_.found("singleProcessorFaceSets"))
-    {
-        List<Tuple2<word, label> > zNameAndProcs
-        (
-            decompositionDict_.lookup("singleProcessorFaceSets")
-        );
-
-        specifiedProcessorFaces.setSize(zNameAndProcs.size());
-        specifiedProcessor.setSize(zNameAndProcs.size());
-
-        forAll(zNameAndProcs, setI)
-        {
-            Info<< "Keeping all cells connected to faceSet "
-                << zNameAndProcs[setI].first()
-                << " on processor " << zNameAndProcs[setI].second() << endl;
-
-            // Read faceSet
-            faceSet fz(mesh, zNameAndProcs[setI].first());
-
-            specifiedProcessorFaces.set(setI, new labelList(fz.sortedToc()));
-            specifiedProcessor[setI] = zNameAndProcs[setI].second();
-        }
-    }
-
-
-    // Construct decomposition method and either do decomposition on
-    // cell centres or on agglomeration
-
-
-    autoPtr<decompositionMethod> decomposePtr = decompositionMethod::New
+    // Any non-mesh connections?
+    label nConnections = returnReduce
     (
-        decompositionDict_
+        explicitConnections.size(),
+        sumOp<label>()
     );
 
+    // Any faces not blocked?
+    label nUnblocked = 0;
+    forAll(blockedFace, faceI)
+    {
+        if (!blockedFace[faceI])
+        {
+            nUnblocked++;
+        }
+    }
+    reduce(nUnblocked, sumOp<label>());
+
+
+
+    // Either do decomposition on cell centres or on agglomeration
 
     labelList finalDecomp;
 
 
-    label nConstraints = returnReduce
-    (
-        sameProcFaces.size()
-      + specifiedProcessorFaces.size(),
-        sumOp<label>()
-    );
-
-
-    label nWeights = returnReduce(cellWeights.size(), sumOp<label>());
-
-    if (nConstraints == 0)
+    if (nProcSets+nConnections+nUnblocked == 0)
     {
+        // No constraints, possibly weights
+
         if (nWeights > 0)
         {
-            finalDecomp = decomposePtr().decompose
+            finalDecomp = decompose
             (
                 mesh,
                 mesh.cellCentres(),
@@ -880,65 +811,24 @@ Foam::labelList Foam::decompositionMethod::decompose
         }
         else
         {
-            finalDecomp = decomposePtr().decompose(mesh, mesh.cellCentres());
+            finalDecomp = decompose(mesh, mesh.cellCentres());
         }
-
     }
     else
     {
-        Info<< "Constrained decomposition:" << endl
-            << "    faces with same processor owner and neighbour : "
-            << sameProcFaces.size() << endl
-            << "    faces all on same processor                   : "
-            << specifiedProcessorFaces.size() << endl << endl;
-
-        // Faces where owner and neighbour are not 'connected' (= all except
-        // sameProcFaces)
-        boolList blockedFace(mesh.nFaces(), true);
-
-        forAllConstIter(labelHashSet, sameProcFaces, iter)
+        if (debug)
         {
-            blockedFace[iter.key()] = false;
-        }
-
-
-        // For specifiedProcessorFaces add all point connected faces
-        forAll(specifiedProcessorFaces, setI)
-        {
-            const labelList& set = specifiedProcessorFaces[setI];
-            forAll(set, fI)
-            {
-                const face& f = mesh.faces()[set[fI]];
-                forAll(f, fp)
-                {
-                    const labelList& pFaces = mesh.pointFaces()[f[fp]];
-                    forAll(pFaces, i)
-                    {
-                        blockedFace[pFaces[i]] = false;
-                    }
-                }
-            }
-        }
-
-
-        // Connect coupled boundary faces
-        const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-        forAll(patches, patchI)
-        {
-            const polyPatch& pp = patches[patchI];
-
-            if (pp.coupled())
-            {
-                forAll(pp, i)
-                {
-                    blockedFace[pp.start()+i] = false;
-                }
-            }
+            Info<< "Constrained decomposition:" << endl
+                << "    faces with same owner and neighbour processor : "
+                << nUnblocked << endl
+                << "    baffle faces with same owner processor        : "
+                << nConnections << endl
+                << "    faces all on same processor                   : "
+                << nProcSets << endl << endl;
         }
 
         // Determine global regions, separated by blockedFaces
-        regionSplit globalRegion(mesh, blockedFace);
+        regionSplit globalRegion(mesh, blockedFace, explicitConnections);
 
 
         // Determine region cell centres
@@ -1044,4 +934,354 @@ Foam::labelList Foam::decompositionMethod::decompose
 }
 
 
+void Foam::decompositionMethod::setConstraints
+(
+    const polyMesh& mesh,
+    boolList& blockedFace,
+    PtrList<labelList>& specifiedProcessorFaces,
+    labelList& specifiedProcessor
+)
+{
+    blockedFace.setSize(mesh.nFaces());
+    blockedFace = true;
+    label nUnblocked = 0;
+
+    specifiedProcessorFaces.clear();
+    specifiedProcessor.setSize(0);
+
+
+    if (decompositionDict_.found("preservePatches"))
+    {
+        wordList pNames(decompositionDict_.lookup("preservePatches"));
+
+        Info<< nl
+            << "Keeping owner of faces in patches " << pNames
+            << " on same processor. This only makes sense for cyclics." << endl;
+
+        const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+        forAll(pNames, i)
+        {
+            const label patchI = patches.findPatchID(pNames[i]);
+
+            if (patchI == -1)
+            {
+                FatalErrorIn("decompositionMethod::decompose(const polyMesh&)")
+                    << "Unknown preservePatch " << pNames[i]
+                    << endl << "Valid patches are " << patches.names()
+                    << exit(FatalError);
+            }
+
+            const polyPatch& pp = patches[patchI];
+
+            forAll(pp, i)
+            {
+                if (blockedFace[pp.start() + i])
+                {
+                    blockedFace[pp.start() + i] = false;
+                    nUnblocked++;
+                }
+            }
+        }
+    }
+    if (decompositionDict_.found("preserveFaceZones"))
+    {
+        wordList zNames(decompositionDict_.lookup("preserveFaceZones"));
+
+        Info<< nl
+            << "Keeping owner and neighbour of faces in zones " << zNames
+            << " on same processor" << endl;
+
+        const faceZoneMesh& fZones = mesh.faceZones();
+
+        forAll(zNames, i)
+        {
+            label zoneI = fZones.findZoneID(zNames[i]);
+
+            if (zoneI == -1)
+            {
+                FatalErrorIn("decompositionMethod::decompose(const polyMesh&)")
+                    << "Unknown preserveFaceZone " << zNames[i]
+                    << endl << "Valid faceZones are " << fZones.names()
+                    << exit(FatalError);
+            }
+
+            const faceZone& fz = fZones[zoneI];
+
+            forAll(fz, i)
+            {
+                if (blockedFace[fz[i]])
+                {
+                    blockedFace[fz[i]] = false;
+                    nUnblocked++;
+                }
+            }
+        }
+    }
+
+    if
+    (
+        decompositionDict_.found("preservePatches")
+     || decompositionDict_.found("preserveFaceZones")
+    )
+    {
+        syncTools::syncFaceList(mesh, blockedFace, andEqOp<bool>());
+        reduce(nUnblocked, sumOp<label>());
+    }
+
+
+
+    // Specified processor for group of cells connected to faces
+
+    label nProcSets = 0;
+    if (decompositionDict_.found("singleProcessorFaceSets"))
+    {
+        List<Tuple2<word, label> > zNameAndProcs
+        (
+            decompositionDict_.lookup("singleProcessorFaceSets")
+        );
+
+        specifiedProcessorFaces.setSize(zNameAndProcs.size());
+        specifiedProcessor.setSize(zNameAndProcs.size());
+
+        forAll(zNameAndProcs, setI)
+        {
+            Info<< "Keeping all cells connected to faceSet "
+                << zNameAndProcs[setI].first()
+                << " on processor " << zNameAndProcs[setI].second() << endl;
+
+            // Read faceSet
+            faceSet fz(mesh, zNameAndProcs[setI].first());
+
+            specifiedProcessorFaces.set(setI, new labelList(fz.sortedToc()));
+            specifiedProcessor[setI] = zNameAndProcs[setI].second();
+            nProcSets += fz.size();
+        }
+        reduce(nProcSets, sumOp<label>());
+
+
+        // Unblock all point connected faces
+        // 1. Mark all points on specifiedProcessorFaces
+        boolList procFacePoint(mesh.nPoints(), false);
+        forAll(specifiedProcessorFaces, setI)
+        {
+            const labelList& set = specifiedProcessorFaces[setI];
+            forAll(set, fI)
+            {
+                const face& f = mesh.faces()[set[fI]];
+                forAll(f, fp)
+                {
+                    procFacePoint[f[fp]] = true;
+                }
+            }
+        }
+        syncTools::syncPointList(mesh, procFacePoint, orEqOp<bool>(), false);
+
+        // 2. Unblock all faces on procFacePoint
+        forAll(procFacePoint, pointI)
+        {
+            if (procFacePoint[pointI])
+            {
+                const labelList& pFaces = mesh.pointFaces()[pointI];
+                forAll(pFaces, i)
+                {
+                    blockedFace[pFaces[i]] = false;
+                }
+            }
+        }
+        syncTools::syncFaceList(mesh, blockedFace, andEqOp<bool>());
+    }
+}
+
+
+Foam::labelList Foam::decompositionMethod::decompose
+(
+    const polyMesh& mesh,
+    const scalarField& cellWeights
+)
+{
+    //labelHashSet sameProcFaces;
+    //
+    //if (decompositionDict_.found("preservePatches"))
+    //{
+    //    wordList pNames(decompositionDict_.lookup("preservePatches"));
+    //
+    //    Info<< nl
+    //        << "Keeping owner of faces in patches " << pNames
+    //        << " on same processor. This only makes sense for cyclics."
+    //        << endl;
+    //
+    //    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+    //
+    //    forAll(pNames, i)
+    //    {
+    //        const label patchI = patches.findPatchID(pNames[i]);
+    //
+    //        if (patchI == -1)
+    //        {
+    //            FatalErrorIn
+    //            (
+    //                "decompositionMethod::decompose(const polyMesh&)")
+    //                << "Unknown preservePatch " << pNames[i]
+    //                << endl << "Valid patches are " << patches.names()
+    //                << exit(FatalError);
+    //        }
+    //
+    //        const polyPatch& pp = patches[patchI];
+    //
+    //        forAll(pp, i)
+    //        {
+    //            sameProcFaces.insert(pp.start() + i);
+    //        }
+    //    }
+    //}
+    //if (decompositionDict_.found("preserveFaceZones"))
+    //{
+    //    wordList zNames(decompositionDict_.lookup("preserveFaceZones"));
+    //
+    //    Info<< nl
+    //        << "Keeping owner and neighbour of faces in zones " << zNames
+    //        << " on same processor" << endl;
+    //
+    //    const faceZoneMesh& fZones = mesh.faceZones();
+    //
+    //    forAll(zNames, i)
+    //    {
+    //        label zoneI = fZones.findZoneID(zNames[i]);
+    //
+    //        if (zoneI == -1)
+    //        {
+    //            FatalErrorIn
+    //            ("decompositionMethod::decompose(const polyMesh&)")
+    //                << "Unknown preserveFaceZone " << zNames[i]
+    //                << endl << "Valid faceZones are " << fZones.names()
+    //                << exit(FatalError);
+    //        }
+    //
+    //        const faceZone& fz = fZones[zoneI];
+    //
+    //        forAll(fz, i)
+    //        {
+    //            sameProcFaces.insert(fz[i]);
+    //        }
+    //    }
+    //}
+    //
+    //
+    //// Specified processor for group of cells connected to faces
+    //
+    ////- Sets of faces to move together
+    //PtrList<labelList> specifiedProcessorFaces;
+    ////- Destination processor
+    //labelList specifiedProcessor;
+    //
+    //label nProcSets = 0;
+    //if (decompositionDict_.found("singleProcessorFaceSets"))
+    //{
+    //    List<Tuple2<word, label> > zNameAndProcs
+    //    (
+    //        decompositionDict_.lookup("singleProcessorFaceSets")
+    //    );
+    //
+    //    specifiedProcessorFaces.setSize(zNameAndProcs.size());
+    //    specifiedProcessor.setSize(zNameAndProcs.size());
+    //
+    //    forAll(zNameAndProcs, setI)
+    //    {
+    //        Info<< "Keeping all cells connected to faceSet "
+    //            << zNameAndProcs[setI].first()
+    //            << " on processor " << zNameAndProcs[setI].second() << endl;
+    //
+    //        // Read faceSet
+    //        faceSet fz(mesh, zNameAndProcs[setI].first());
+    //
+    //        specifiedProcessorFaces.set(setI, new labelList(fz.sortedToc()));
+    //        specifiedProcessor[setI] = zNameAndProcs[setI].second();
+    //        nProcSets += fz.size();
+    //    }
+    //    reduce(nProcSets, sumOp<label>());
+    //}
+    //
+    //
+    //label nUnblocked = returnReduce(sameProcFaces.size(), sumOp<label>());
+    //
+    //if (nProcSets+nUnblocked > 0)
+    //{
+    //    Info<< "Constrained decomposition:" << endl
+    //        << "    faces with same owner and neighbour processor : "
+    //        << nUnblocked << endl
+    //        << "    faces all on same processor                   : "
+    //        << nProcSets << endl << endl;
+    //}
+    //
+    //
+    //// Faces where owner and neighbour are not 'connected' (= all except
+    //// sameProcFaces)
+    //boolList blockedFace(mesh.nFaces(), true);
+    //{
+    //    forAllConstIter(labelHashSet, sameProcFaces, iter)
+    //    {
+    //        blockedFace[iter.key()] = false;
+    //    }
+    //    syncTools::syncFaceList(mesh, blockedFace, andEqOp<bool>());
+    //
+    //    // Add all point connected faces
+    //    boolList procFacePoint(mesh.nPoints(), false);
+    //    forAll(specifiedProcessorFaces, setI)
+    //    {
+    //        const labelList& set = specifiedProcessorFaces[setI];
+    //        forAll(set, fI)
+    //        {
+    //            const face& f = mesh.faces()[set[fI]];
+    //            forAll(f, fp)
+    //            {
+    //                procFacePoint[f[fp]] = true;
+    //            }
+    //        }
+    //    }
+    //    syncTools::syncPointList(mesh, procFacePoint, orEqOp<bool>(), false);
+    //
+    //    forAll(procFacePoint, pointI)
+    //    {
+    //        if (procFacePoint[pointI])
+    //        {
+    //            const labelList& pFaces = mesh.pointFaces()[pointI];
+    //            forAll(pFaces, i)
+    //            {
+    //                blockedFace[pFaces[i]] = false;
+    //            }
+    //        }
+    //    }
+    //    syncTools::syncFaceList(mesh, blockedFace, andEqOp<bool>());
+    //}
+
+    boolList blockedFace;
+    PtrList<labelList> specifiedProcessorFaces;
+    labelList specifiedProcessor;
+    setConstraints
+    (
+        mesh,
+        blockedFace,
+        specifiedProcessorFaces,
+        specifiedProcessor
+    );
+
+
+    // Construct decomposition method and either do decomposition on
+    // cell centres or on agglomeration
+
+    labelList finalDecomp = decompose
+    (
+        mesh,
+        cellWeights,            // optional weights
+        blockedFace,            // any cells to be combined
+        specifiedProcessorFaces,// any whole cluster of cells to be kept
+        specifiedProcessor,
+        List<labelPair>()       // no baffles
+    );
+
+    return finalDecomp;
+}
+
+
 // ************************************************************************* //
diff --git a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H
index 2a1921778c..43e538dc4e 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H
+++ b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H
@@ -255,6 +255,34 @@ public:
             //    CompactListList<label>& cellCells
             //);
 
+
+            //- Helper: extract constraints
+            void setConstraints
+            (
+                const polyMesh& mesh,
+                boolList& blockedFace,
+                PtrList<labelList>& specifiedProcessorFaces,
+                labelList& specifiedProcessor
+            );
+
+            // Decompose a mesh with constraints:
+            // - blockedFace : whether owner and neighbour should be on
+            // - specifiedProcessorFaces, specifiedProcessor : sets of faces
+            //   that should go to same processor (as specified in
+            //   specifiedProcessor, can be -1)
+            // - explicitConnections : connections between baffle faces
+            // Set all to zero size to have unconstrained decomposition.
+            labelList decompose
+            (
+                const polyMesh& mesh,
+                const scalarField& cellWeights,
+                const boolList& blockedFace,
+                const PtrList<labelList>& specifiedProcessorFaces,
+                const labelList& specifiedProcessor,
+                const List<labelPair>& explicitConnections
+            );
+
+
             //- Decompose a mesh. Apply all constraints from decomposeParDict
             //  ('preserveFaceZones' etc). Calls either
             //  - no constraints, empty weights: