diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposeParDict b/applications/utilities/parallelProcessing/decomposePar/decomposeParDict
index e60757fffd..c42bdc1c9a 100644
--- a/applications/utilities/parallelProcessing/decomposePar/decomposeParDict
+++ b/applications/utilities/parallelProcessing/decomposePar/decomposeParDict
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-numberOfSubdomains  8;
+numberOfSubdomains  2;
 
 //- Keep owner and neighbour on same processor for faces in zones:
 // preserveFaceZones (heater solid1 solid3);
@@ -24,13 +24,22 @@ numberOfSubdomains  8;
 //  (makes sense only for cyclic patches)
 //preservePatches (cyclic_half0 cyclic_half1);
 
+//- Keep all of faceZone on a single processor. This puts all cells
+//  connected with a point, edge or face on the same processor.
+//  (just having face connected cells might not guarantee a balanced
+//  decomposition)
+// The processor can be explicitly provided or -1 to have
+// decompositionMethod choose.
+//singleProcessorFaceZones ((f0 -1));
+
+
 //- Use the volScalarField named here as a weight for each cell in the
 //  decomposition.  For example, use a particle population field to decompose
 //  for a balanced number of particles in a lagrangian simulation.
 // weightField dsmcRhoNMean;
 
 method          scotch;
-// method          hierarchical;
+//method          hierarchical;
 // method          simple;
 // method          metis;
 // method          manual;
@@ -70,7 +79,7 @@ simpleCoeffs
 
 hierarchicalCoeffs
 {
-    n           (2 2 1);
+    n           (1 2 1);
     delta       0.001;
     order       xyz;
 }
diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C
index 4ee3c7555e..0d40f8b232 100644
--- a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C
+++ b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C
@@ -28,6 +28,7 @@ License
 #include "cpuTime.H"
 #include "cellSet.H"
 #include "regionSplit.H"
+#include "Tuple2.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,7 +48,8 @@ void Foam::domainDecomposition::distributeCells()
     {
         wordList pNames(decompositionDict_.lookup("preservePatches"));
 
-        Info<< "Keeping owner of faces in patches " << pNames
+        Info<< nl
+            << "Keeping owner of faces in patches " << pNames
             << " on same processor. This only makes sense for cyclics." << endl;
 
         const polyBoundaryMesh& patches = boundaryMesh();
@@ -76,7 +78,8 @@ void Foam::domainDecomposition::distributeCells()
     {
         wordList zNames(decompositionDict_.lookup("preserveFaceZones"));
 
-        Info<< "Keeping owner and neighbour of faces in zones " << zNames
+        Info<< nl
+            << "Keeping owner and neighbour of faces in zones " << zNames
             << " on same processor" << endl;
 
         const faceZoneMesh& fZones = faceZones();
@@ -103,6 +106,67 @@ void Foam::domainDecomposition::distributeCells()
     }
 
 
+    // Specified processor for owner and neighbour of faces
+    Map<label> specifiedProcessorFaces;
+
+    if (decompositionDict_.found("singleProcessorFaceZones"))
+    {
+        List<Tuple2<word, label> > zNameAndProcs
+        (
+            decompositionDict_.lookup("singleProcessorFaceZones")
+        );
+
+        const faceZoneMesh& fZones = faceZones();
+
+        label nCells = 0;
+
+        Info<< endl;
+
+        forAll(zNameAndProcs, i)
+        {
+            Info<< "Keeping all cells connected to faceZone "
+                << zNameAndProcs[i].first()
+                << " on processor " << zNameAndProcs[i].second() << endl;
+
+            label zoneI = fZones.findZoneID(zNameAndProcs[i].first());
+
+            if (zoneI == -1)
+            {
+                FatalErrorIn("domainDecomposition::distributeCells()")
+                    << "Unknown singleProcessorFaceZone "
+                    << zNameAndProcs[i].first()
+                    << endl << "Valid faceZones are " << fZones.names()
+                    << exit(FatalError);
+            }
+
+            const faceZone& fz = fZones[zoneI];
+
+            nCells += fz.size();
+        }
+
+
+        // Size
+        specifiedProcessorFaces.resize(2*nCells);
+
+
+        // Fill
+        forAll(zNameAndProcs, i)
+        {
+            label zoneI = fZones.findZoneID(zNameAndProcs[i].first());
+            const faceZone& fz = fZones[zoneI];
+
+            label procI = zNameAndProcs[i].second();
+
+            forAll(fz, fzI)
+            {
+                label faceI = fz[fzI];
+
+                specifiedProcessorFaces.insert(faceI, procI);
+            }
+        }
+    }
+
+
     // Construct decomposition method and either do decomposition on
     // cell centres or on agglomeration
 
@@ -112,7 +176,8 @@ void Foam::domainDecomposition::distributeCells()
         decompositionDict_
     );
 
-    if (sameProcFaces.empty())
+
+    if (sameProcFaces.empty() && specifiedProcessorFaces.empty())
     {
         if (decompositionDict_.found("weightField"))
         {
@@ -146,9 +211,11 @@ void Foam::domainDecomposition::distributeCells()
     }
     else
     {
-        Info<< "Selected " << sameProcFaces.size()
-            << " faces whose owner and neighbour cell should be kept on the"
-            << " same processor" << endl;
+        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)
@@ -159,6 +226,24 @@ void Foam::domainDecomposition::distributeCells()
             blockedFace[iter.key()] = false;
         }
 
+
+        // For specifiedProcessorFaces add all point connected faces
+        {
+            forAllConstIter(Map<label>, specifiedProcessorFaces, iter)
+            {
+                const face& f = faces()[iter.key()];
+                forAll(f, fp)
+                {
+                    const labelList& pFaces = pointFaces()[f[fp]];
+                    forAll(pFaces, i)
+                    {
+                        blockedFace[pFaces[i]] = false;
+                    }
+                }
+            }
+        }
+
+
         // Connect coupled boundary faces
         const polyBoundaryMesh& patches =  boundaryMesh();
 
@@ -201,10 +286,11 @@ void Foam::domainDecomposition::distributeCells()
 
         // Do decomposition on agglomeration
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        scalarField regionWeights(globalRegion.nRegions(), 0);
+
         if (decompositionDict_.found("weightField"))
         {
-            scalarField regionWeights(globalRegion.nRegions(), 0);
-
             word weightName = decompositionDict_.lookup("weightField");
 
             volScalarField weights
@@ -226,23 +312,46 @@ void Foam::domainDecomposition::distributeCells()
 
                 regionWeights[regionI] += weights.internalField()[cellI];
             }
-
-            cellToProc_ = decomposePtr().decompose
-            (
-                *this,
-                globalRegion,
-                regionCentres,
-                regionWeights
-            );
         }
         else
         {
-            cellToProc_ = decomposePtr().decompose
-            (
-                *this,
-                globalRegion,
-                regionCentres
-            );
+            forAll(globalRegion, cellI)
+            {
+                label regionI = globalRegion[cellI];
+
+                regionWeights[regionI] += 1.0;
+            }
+        }
+
+        cellToProc_ = decomposePtr().decompose
+        (
+            *this,
+            globalRegion,
+            regionCentres,
+            regionWeights
+        );
+
+
+        // For specifiedProcessorFaces rework the cellToProc. Note that
+        // this makes the decomposition unbalanced.
+        if (specifiedProcessorFaces.size())
+        {
+            labelList procMap(identity(cellToProc_.size()));
+
+            forAllConstIter(Map<label>, specifiedProcessorFaces, iter)
+            {
+                label faceI = iter.key();
+                label wantedProcI = iter();
+
+                if (wantedProcI != -1)
+                {
+                    cellToProc_[faceOwner()[faceI]] = wantedProcI;
+                    if (isInternalFace(faceI))
+                    {
+                        cellToProc_[faceNeighbour()[faceI]] = wantedProcI;
+                    }
+                }
+            }
         }
     }
 
diff --git a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
index 663f5b003b..1808fc84ab 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
+++ b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
@@ -403,6 +403,7 @@ void Foam::decompositionMethod::calcCellCells
     forAll(cellCells, cellI)
     {
         nbrCells.clear();
+        nbrCells.insert(cellI);
 
         label& endIndex = cellCells.offsets()[cellI+1];