From 927c94c6d858cea9caa81c0c13ab20c30f6adae9 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Fri, 14 Oct 2011 13:00:53 +0100
Subject: [PATCH] Patch constrained decomposition

---
 .../decompositionMethods/Make/files           |   3 +-
 .../decompositionMethod/decompositionMethod.C | 437 ++++++++++++++----
 .../decompositionMethod/decompositionMethod.H |  90 ++--
 .../geomDecomp/geomDecomp.C                   |   3 +-
 .../geomDecomp/geomDecomp.H                   |  19 +-
 .../hierarchGeomDecomp/hierarchGeomDecomp.C   |  21 +-
 .../hierarchGeomDecomp/hierarchGeomDecomp.H   |  43 +-
 .../manualDecomp/manualDecomp.C               |  32 +-
 .../manualDecomp/manualDecomp.H               |  21 +-
 .../patchConstrainedDecomp.C                  | 120 +++++
 .../patchConstrainedDecomp.H                  | 142 ++++++
 .../simpleGeomDecomp/simpleGeomDecomp.C       |   8 +-
 .../simpleGeomDecomp/simpleGeomDecomp.H       |  10 +-
 13 files changed, 753 insertions(+), 196 deletions(-)
 create mode 100644 src/decompositionMethods/decompositionMethods/patchConstrainedDecomp/patchConstrainedDecomp.C
 create mode 100644 src/decompositionMethods/decompositionMethods/patchConstrainedDecomp/patchConstrainedDecomp.H

diff --git a/src/decompositionMethods/decompositionMethods/Make/files b/src/decompositionMethods/decompositionMethods/Make/files
index 43bb9a25c..28ba31a12 100644
--- a/src/decompositionMethods/decompositionMethods/Make/files
+++ b/src/decompositionMethods/decompositionMethods/Make/files
@@ -1,7 +1,8 @@
 decompositionMethod/decompositionMethod.C
+manualDecomp/manualDecomp.C
 geomDecomp/geomDecomp.C
 simpleGeomDecomp/simpleGeomDecomp.C
 hierarchGeomDecomp/hierarchGeomDecomp.C
-manualDecomp/manualDecomp.C
+patchConstrainedDecomp/patchConstrainedDecomp.C
 
 LIB = $(FOAM_LIBBIN)/libdecompositionMethods
diff --git a/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.C b/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.C
index 05a206020..5adf0f20c 100644
--- a/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.C
+++ b/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.C
@@ -28,6 +28,9 @@ InClass
 \*---------------------------------------------------------------------------*/
 
 #include "decompositionMethod.H"
+#include "cyclicPolyPatch.H"
+#include "syncTools.H"
+#include "globalIndex.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -38,7 +41,336 @@ namespace Foam
     defineRunTimeSelectionTable(decompositionMethod, dictionaryMesh);
 }
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::decompositionMethod::calcCSR
+(
+    const labelListList& cellCells,
+    List<int>& adjncy,
+    List<int>& xadj
+)
+{
+    // Count number of internal faces
+    label nConnections = 0;
+
+    forAll(cellCells, coarseI)
+    {
+        nConnections += cellCells[coarseI].size();
+    }
+
+    // Create the adjncy array as twice the size of the total number of
+    // internal faces
+    adjncy.setSize(nConnections);
+
+    xadj.setSize(cellCells.size() + 1);
+
+
+    // Fill in xadj
+    // ~~~~~~~~~~~~
+    label freeAdj = 0;
+
+    forAll(cellCells, coarseI)
+    {
+        xadj[coarseI] = freeAdj;
+
+        const labelList& cCells = cellCells[coarseI];
+
+        forAll(cCells, i)
+        {
+            adjncy[freeAdj++] = cCells[i];
+        }
+    }
+    xadj[cellCells.size()] = freeAdj;
+}
+
+
+
+void Foam::decompositionMethod::calcCSR
+(
+    const polyMesh& mesh,
+    List<int>& adjncy,
+    List<int>& xadj
+)
+{
+    // Make Metis CSR (Compressed Storage Format) storage
+    //   adjncy      : contains neighbours (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
+
+    xadj.setSize(mesh.nCells() + 1);
+
+    // Initialise the number of internal faces of the cells to twice the
+    // number of internal faces
+    label nInternalFaces = 2*mesh.nInternalFaces();
+
+    // Check the boundary for coupled patches and add to the number of
+    // internal faces
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+    forAll(pbm, patchi)
+    {
+        if (isA<cyclicPolyPatch>(pbm[patchi]))
+        {
+            nInternalFaces += pbm[patchi].size();
+        }
+    }
+
+    // Create the adjncy array the size of the total number of internal and
+    // coupled faces
+    adjncy.setSize(nInternalFaces);
+
+    // Fill in xadj
+    // ~~~~~~~~~~~~
+    label freeAdj = 0;
+
+    for (label cellI = 0; cellI < mesh.nCells(); cellI++)
+    {
+        xadj[cellI] = freeAdj;
+
+        const labelList& cFaces = mesh.cells()[cellI];
+
+        forAll(cFaces, i)
+        {
+            label faceI = cFaces[i];
+
+            if
+            (
+                mesh.isInternalFace(faceI)
+             || isA<cyclicPolyPatch>(pbm[pbm.whichPatch(faceI)])
+            )
+            {
+                freeAdj++;
+            }
+        }
+    }
+    xadj[mesh.nCells()] = freeAdj;
+
+
+    // Fill in adjncy
+    // ~~~~~~~~~~~~~~
+
+    labelList nFacesPerCell(mesh.nCells(), 0);
+
+    // Internal faces
+    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+    {
+        label own = mesh.faceOwner()[faceI];
+        label nei = mesh.faceNeighbour()[faceI];
+
+        adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
+        adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
+    }
+
+    // Coupled faces. Only cyclics done.
+    forAll(pbm, patchi)
+    {
+        if (isA<cyclicPolyPatch>(pbm[patchi]))
+        {
+            const unallocLabelList& faceCells = pbm[patchi].faceCells();
+
+            label sizeby2 = faceCells.size()/2;
+
+            for (label facei=0; facei<sizeby2; facei++)
+            {
+                label own = faceCells[facei];
+                label nei = faceCells[facei + sizeby2];
+
+                adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
+                adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
+            }
+        }
+    }
+}
+
+
+void Foam::decompositionMethod::calcDistributedCSR
+(
+    const polyMesh& mesh,
+    List<int>& adjncy,
+    List<int>& xadj
+)
+{
+    // Create global cell numbers
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    globalIndex globalCells(mesh.nCells());
+
+
+    //
+    // Make Metis Distributed CSR (Compressed Storage Format) storage
+    //   adjncy      : contains cellCells (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
+    //
+
+
+    const labelList& faceOwner = mesh.faceOwner();
+    const labelList& faceNeighbour = mesh.faceNeighbour();
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+
+    // Get renumbered owner on other side of coupled faces
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    List<int> globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+            label bFaceI = pp.start() - mesh.nInternalFaces();
+
+            forAll(pp, i)
+            {
+                globalNeighbour[bFaceI++] = globalCells.toGlobal
+                (
+                    faceOwner[faceI++]
+                );
+            }
+        }
+    }
+
+    // Get the cell on the other side of coupled patches
+    syncTools::swapBoundaryFaceList(mesh, globalNeighbour, false);
+
+
+    // Count number of faces (internal + coupled)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Number of faces per cell
+    List<int> nFacesPerCell(mesh.nCells(), 0);
+
+    // Number of coupled faces
+    label nCoupledFaces = 0;
+
+    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+    {
+        nFacesPerCell[faceOwner[faceI]]++;
+        nFacesPerCell[faceNeighbour[faceI]]++;
+    }
+    // Handle coupled faces
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+
+            forAll(pp, i)
+            {
+                nCoupledFaces++;
+                nFacesPerCell[faceOwner[faceI++]]++;
+            }
+        }
+    }
+
+
+    // Fill in xadj
+    // ~~~~~~~~~~~~
+
+    xadj.setSize(mesh.nCells() + 1);
+
+    int freeAdj = 0;
+
+    for (label cellI = 0; cellI < mesh.nCells(); cellI++)
+    {
+        xadj[cellI] = freeAdj;
+
+        freeAdj += nFacesPerCell[cellI];
+    }
+    xadj[mesh.nCells()] = freeAdj;
+
+
+
+    // Fill in adjncy
+    // ~~~~~~~~~~~~~~
+
+    adjncy.setSize(2*mesh.nInternalFaces() + nCoupledFaces);
+
+    nFacesPerCell = 0;
+
+    // For internal faces is just offsetted owner and neighbour
+    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+    {
+        label own = faceOwner[faceI];
+        label nei = faceNeighbour[faceI];
+
+        adjncy[xadj[own] + nFacesPerCell[own]++] = globalCells.toGlobal(nei);
+        adjncy[xadj[nei] + nFacesPerCell[nei]++] = globalCells.toGlobal(own);
+    }
+    // For boundary faces is offsetted coupled neighbour
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+            label bFaceI = pp.start()-mesh.nInternalFaces();
+
+            forAll(pp, i)
+            {
+                label own = faceOwner[faceI];
+                adjncy[xadj[own] + nFacesPerCell[own]++] =
+                    globalNeighbour[bFaceI];
+
+                faceI++;
+                bFaceI++;
+            }
+        }
+    }
+}
+
+
+void Foam::decompositionMethod::calcCellCells
+(
+    const polyMesh& mesh,
+    const labelList& fineToCoarse,
+    const label nCoarse,
+    labelListList& cellCells
+)
+{
+    if (fineToCoarse.size() != mesh.nCells())
+    {
+        FatalErrorIn
+        (
+            "decompositionMethod::calcCellCells"
+            "(const labelList&, labelListList&) const"
+        )   << "Only valid for mesh agglomeration." << exit(FatalError);
+    }
+
+    List<DynamicList<label> > dynCellCells(nCoarse);
+
+    forAll(mesh.faceNeighbour(), faceI)
+    {
+        label own = fineToCoarse[mesh.faceOwner()[faceI]];
+        label nei = fineToCoarse[mesh.faceNeighbour()[faceI]];
+
+        if (own != nei)
+        {
+            if (findIndex(dynCellCells[own], nei) == -1)
+            {
+                dynCellCells[own].append(nei);
+            }
+            if (findIndex(dynCellCells[nei], own) == -1)
+            {
+                dynCellCells[nei].append(own);
+            }
+        }
+    }
+
+    cellCells.setSize(dynCellCells.size());
+    forAll(dynCellCells, coarseI)
+    {
+        cellCells[coarseI].transfer(dynCellCells[coarseI]);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 Foam::autoPtr<Foam::decompositionMethod> Foam::decompositionMethod::New
 (
@@ -102,17 +434,40 @@ Foam::autoPtr<Foam::decompositionMethod> Foam::decompositionMethod::New
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
 Foam::labelList Foam::decompositionMethod::decompose
 (
     const pointField& points
 )
 {
-    scalarField weights(0);
+    scalarField weights(points.size(), 1);
 
     return decompose(points, weights);
 }
 
 
+Foam::labelList Foam::decompositionMethod::decompose
+(
+    const labelList& fineToCoarse,
+    const pointField& coarsePoints
+)
+{
+    // Decompose based on agglomerated points
+    labelList coarseDistribution(decompose(coarsePoints));
+
+    // Rework back into decomposition for original mesh_
+    labelList fineDistribution(fineToCoarse.size());
+
+    forAll(fineDistribution, i)
+    {
+        fineDistribution[i] = coarseDistribution[fineToCoarse[i]];
+    }
+
+    return fineDistribution;
+}
+
+
 Foam::labelList Foam::decompositionMethod::decompose
 (
     const labelList& fineToCoarse,
@@ -135,82 +490,4 @@ Foam::labelList Foam::decompositionMethod::decompose
 }
 
 
-Foam::labelList Foam::decompositionMethod::decompose
-(
-    const labelList& fineToCoarse,
-    const pointField& coarsePoints
-)
-{
-    // Decompose based on agglomerated points
-    labelList coarseDistribution(decompose(coarsePoints));
-
-    // Rework back into decomposition for original mesh_
-    labelList fineDistribution(fineToCoarse.size());
-
-    forAll(fineDistribution, i)
-    {
-        fineDistribution[i] = coarseDistribution[fineToCoarse[i]];
-    }
-
-    return fineDistribution;
-}
-
-
-void Foam::decompositionMethod::calcCellCells
-(
-    const polyMesh& mesh,
-    const labelList& fineToCoarse,
-    const label nCoarse,
-    labelListList& cellCells
-)
-{
-    if (fineToCoarse.size() != mesh.nCells())
-    {
-        FatalErrorIn
-        (
-            "decompositionMethod::calcCellCells"
-            "(const labelList&, labelListList&) const"
-        )   << "Only valid for mesh agglomeration." << exit(FatalError);
-    }
-
-    List<DynamicList<label> > dynCellCells(nCoarse);
-
-    forAll(mesh.faceNeighbour(), faceI)
-    {
-        label own = fineToCoarse[mesh.faceOwner()[faceI]];
-        label nei = fineToCoarse[mesh.faceNeighbour()[faceI]];
-
-        if (own != nei)
-        {
-            if (findIndex(dynCellCells[own], nei) == -1)
-            {
-                dynCellCells[own].append(nei);
-            }
-            if (findIndex(dynCellCells[nei], own) == -1)
-            {
-                dynCellCells[nei].append(own);
-            }
-        }
-    }
-
-    cellCells.setSize(dynCellCells.size());
-    forAll(dynCellCells, coarseI)
-    {
-        cellCells[coarseI].transfer(dynCellCells[coarseI]);
-    }
-}
-
-
-Foam::labelList Foam::decompositionMethod::decompose
-(
-    const labelListList& globalCellCells,
-    const pointField& cc
-)
-{
-    scalarField cWeights(0);
-
-    return decompose(globalCellCells, cc, cWeights);
-}
-
-
 // ************************************************************************* //
diff --git a/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.H b/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.H
index e241ffd63..7957a96a6 100644
--- a/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.H
+++ b/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.H
@@ -48,30 +48,64 @@ namespace Foam
 
 class decompositionMethod
 {
-
 protected:
 
     // Protected data
 
+        //- Decomposition dictionary
         const dictionary& decompositionDict_;
+
+        //- Number of processors in decomposition
         label nProcessors_;
 
 
-        //- Helper: determine (non-parallel) cellCells from mesh agglomeration.
+        //- Helper to convert connectivity supplied as cellCells into
+        //  simple CSR (Metis, scotch) storage
+        static void calcCSR
+        (
+            const labelListList& globalCellCells,
+            List<int>& adjncy,
+            List<int>& xadj
+        );
+
+        //- Helper: convert local connectivity from the mesh
+        //  into CSR (Metis, scotch) storage
+        //  Treats cyclics as coupled, but not processor patches
+        static void calcCSR
+        (
+            const polyMesh& mesh,
+            List<int>& adjncy,
+            List<int>& xadj
+        );
+
+        //- Helper: convert mesh connectivity into distributed CSR
+        //  Very dubious coding.  HJ, 1/Mar/2011
+        static void calcDistributedCSR
+        (
+            const polyMesh&,
+            List<int>& adjncy,
+            List<int>& xadj
+        );
+
+        //- Helper: determine (non-parallel) cellCells from mesh
+        // with agglomeration
         static void calcCellCells
         (
             const polyMesh& mesh,
-            const labelList& agglom,
+            const labelList& fineToCoarse,
             const label nCoarse,
             labelListList& cellCells
         );
 
+
 private:
 
     // Private Member Functions
 
-        //- Disallow default bitwise copy construct and assignment
+        //- Disallow default bitwise copy construct
         decompositionMethod(const decompositionMethod&);
+
+        //- Disallow default bitwise assignment
         void operator=(const decompositionMethod&);
 
 
@@ -109,7 +143,7 @@ public:
 
     // Selectors
 
-        //- Return a reference to the selected decomposition method
+        //- Return a pointer to the selected decomposition method
         static autoPtr<decompositionMethod> New
         (
             const dictionary& decompositionDict
@@ -148,46 +182,37 @@ public:
         //  proc boundaries)
         virtual bool parallelAware() const = 0;
 
-        //- Return for every coordinate the wanted processor number. Use the
-        //  mesh connectivity (if needed)
+        //- Decompose cells
+        //  If needed, use connectivity directly from the mesh
+        //  Calls decompose (below) with uniform weights
+        virtual labelList decompose(const pointField&);
+
+        //- Decompose cells with weights
         virtual labelList decompose
         (
             const pointField& points,
             const scalarField& pointWeights
         ) = 0;
 
-        //- Like decompose but with uniform weights on the points
-        virtual labelList decompose(const pointField&);
 
-
-        //- Return for every coordinate the wanted processor number. Gets
-        //  passed agglomeration map (from fine to coarse cells) and coarse cell
-        //  location. Can be overridden by decomposers that provide this
-        //  functionality natively. Coarse cells are local to the processor
-        //  (if in parallel). If you want to have coarse cells spanning
-        //  processors use the globalCellCells instead.
+        //- Decompose cell clusters
+        //  Calls decompose (below) with uniform weights
         virtual labelList decompose
         (
-            const labelList& cellToRegion,
-            const pointField& regionPoints,
-            const scalarField& regionWeights
+            const labelList& fineToCoarse,
+            const pointField& coarsePoints
         );
 
-        //- Like decompose but with uniform weights on the regions
+        //- Decompose cell clusters with weights on clusters
         virtual labelList decompose
         (
-            const labelList& cellToRegion,
-            const pointField& regionPoints
+            const labelList& fineToCoarse,
+            const pointField& coarsePoints,
+            const scalarField& coarseWeights
         );
 
 
-        //- Return for every coordinate the wanted processor number. Explicitly
-        //  provided connectivity - does not use mesh_.
-        //  The connectivity is equal to mesh.cellCells() except for
-        //  - in parallel the cell numbers are global cell numbers (starting
-        //    from 0 at processor0 and then incrementing all through the
-        //    processors)
-        //  - the connections are across coupled patches
+        //- Decompose cells with weights with explicitly provided connectivity
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
@@ -195,13 +220,6 @@ public:
             const scalarField& cWeights
         ) = 0;
 
-        //- Like decompose but with uniform weights on the cells
-        virtual labelList decompose
-        (
-            const labelListList& globalCellCells,
-            const pointField& cc
-        );
-
 };
 
 
diff --git a/src/decompositionMethods/decompositionMethods/geomDecomp/geomDecomp.C b/src/decompositionMethods/decompositionMethods/geomDecomp/geomDecomp.C
index 4202c9b2a..bc6697a55 100644
--- a/src/decompositionMethods/decompositionMethods/geomDecomp/geomDecomp.C
+++ b/src/decompositionMethods/decompositionMethods/geomDecomp/geomDecomp.C
@@ -40,8 +40,7 @@ Foam::geomDecomp::geomDecomp
     delta_(readScalar(geomDecomDict_.lookup("delta"))),
     rotDelta_(I)
 {
-    // check that the case makes sense :
-
+    // Check that the case makes sense
     if (nProcessors_ != n_.x()*n_.y()*n_.z())
     {
         FatalErrorIn
diff --git a/src/decompositionMethods/decompositionMethods/geomDecomp/geomDecomp.H b/src/decompositionMethods/decompositionMethods/geomDecomp/geomDecomp.H
index dee0ea40f..060898e94 100644
--- a/src/decompositionMethods/decompositionMethods/geomDecomp/geomDecomp.H
+++ b/src/decompositionMethods/decompositionMethods/geomDecomp/geomDecomp.H
@@ -37,7 +37,7 @@ SourceFiles
 #define geomDecomp_H
 
 #include "decompositionMethod.H"
-#include "Vector.H"
+#include "labelVector.H"
 
 namespace Foam
 {
@@ -50,17 +50,32 @@ class geomDecomp
 :
     public decompositionMethod
 {
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        geomDecomp(const geomDecomp&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const geomDecomp&);
+
 
 protected:
 
     // Protected data
 
+        //- Geometric decomposition dictionary
         const dictionary& geomDecomDict_;
 
-        Vector<label> n_;
+        //- Number of splits in the (x y z) direction
+        labelVector n_;
+
+        //- Small rotation to achieve smooth geometric decomposition
         scalar delta_;
+
+        //- Small rotation tensor, calculated from delta
         tensor rotDelta_;
 
+
 public:
 
     // Constructors
diff --git a/src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C b/src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
index 3b46977ea..9ed680ca4 100644
--- a/src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
+++ b/src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
@@ -50,6 +50,7 @@ namespace Foam
     );
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 void Foam::hierarchGeomDecomp::setDecompOrder()
@@ -157,7 +158,7 @@ void Foam::hierarchGeomDecomp::calculateSortedWeightedSizes
 {
     // Evaluate cumulative weights.
     sortedWeightedSizes[0] = 0;
-    forAll(current, i)
+    forAll (current, i)
     {
         label pointI = current[indices[i]];
         sortedWeightedSizes[i + 1] = sortedWeightedSizes[i] + weights[pointI];
@@ -328,7 +329,7 @@ void Foam::hierarchGeomDecomp::findBinary
 }
 
 
-// Sort points into bins according to one component. Recurses to next component.
+// Sort points into bins according to one component. Recurses to next component
 void Foam::hierarchGeomDecomp::sortComponent
 (
     const label sizeTol,
@@ -351,7 +352,7 @@ void Foam::hierarchGeomDecomp::sortComponent
     // Storage for sorted component compI
     SortableList<scalar> sortedCoord(current.size());
 
-    forAll(current, i)
+    forAll (current, i)
     {
         label pointI = current[i];
 
@@ -458,7 +459,7 @@ void Foam::hierarchGeomDecomp::sortComponent
         // Copy localSize elements starting from leftIndex.
         labelList slice(localSize);
 
-        forAll(slice, i)
+        forAll (slice, i)
         {
             label pointI = current[sortedCoord.indices()[leftIndex+i]];
 
@@ -502,7 +503,7 @@ void Foam::hierarchGeomDecomp::sortComponent
 }
 
 
-// Sort points into bins according to one component. Recurses to next component.
+// Sort points into bins according to one component. Recurses to next component
 void Foam::hierarchGeomDecomp::sortComponent
 (
     const label sizeTol,
@@ -526,7 +527,7 @@ void Foam::hierarchGeomDecomp::sortComponent
     // Storage for sorted component compI
     SortableList<scalar> sortedCoord(current.size());
 
-    forAll(current, i)
+    forAll (current, i)
     {
         label pointI = current[i];
 
@@ -641,7 +642,7 @@ void Foam::hierarchGeomDecomp::sortComponent
         // Copy localSize elements starting from leftIndex.
         labelList slice(localSize);
 
-        forAll(slice, i)
+        forAll (slice, i)
         {
             label pointI = current[sortedCoord.indices()[leftIndex+i]];
 
@@ -725,7 +726,8 @@ Foam::labelList Foam::hierarchGeomDecomp::decompose
 
     // Start off with every point sorted onto itself.
     labelList slice(points.size());
-    forAll(slice, i)
+
+    forAll (slice, i)
     {
         slice[i] = i;
     }
@@ -767,7 +769,8 @@ Foam::labelList Foam::hierarchGeomDecomp::decompose
 
     // Start off with every point sorted onto itself.
     labelList slice(points.size());
-    forAll(slice, i)
+
+    forAll (slice, i)
     {
         slice[i] = i;
     }
diff --git a/src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H b/src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H
index 7304b45c6..797341745 100644
--- a/src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H
+++ b/src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H
@@ -77,6 +77,13 @@ class hierarchGeomDecomp
 
     // Private Member Functions
 
+        //- Disallow default bitwise copy construct
+        hierarchGeomDecomp(const hierarchGeomDecomp&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const hierarchGeomDecomp&);
+
+
         //- Convert ordering string ("xyz") into list of components.
         void setDecompOrder();
 
@@ -105,13 +112,13 @@ class hierarchGeomDecomp
         //  wantedSize. Binary search.
         static void findBinary
         (
-            const label sizeTol,        // size difference considered acceptible
+            const label sizeTol,        // Acceptable size difference
             const List<scalar>&,
             const label leftIndex,      // index of previous value
             const scalar leftValue,     // value at leftIndex
             const scalar maxValue,      // global max of values
             const scalar wantedSize,    // wanted size
-            label& mid,                 // index where size of bin is wantedSize
+            label& mid,                 // index where bin size is wantedSize
             scalar& midValue            // value at mid
         );
 
@@ -120,14 +127,14 @@ class hierarchGeomDecomp
         //  wantedSize. Binary search.
         static void findBinary
         (
-            const label sizeTol,        // size difference considered acceptible
+            const label sizeTol,        // Acceptable size difference
             const List<scalar>& sortedWeightedSizes,
             const List<scalar>&,
             const label leftIndex,      // index of previous value
             const scalar leftValue,     // value at leftIndex
             const scalar maxValue,      // global max of values
             const scalar wantedSize,    // wanted size
-            label& mid,                 // index where size of bin is wantedSize
+            label& mid,                 // index where bin size is wantedSize
             scalar& midValue            // value at mid
         );
 
@@ -156,11 +163,6 @@ class hierarchGeomDecomp
         );
 
 
-        //- Disallow default bitwise copy construct and assignment
-        void operator=(const hierarchGeomDecomp&);
-        hierarchGeomDecomp(const hierarchGeomDecomp&);
-
-
 public:
 
     //- Runtime type information
@@ -194,25 +196,19 @@ public:
             return true;
         }
 
-        //- Return for every coordinate the wanted processor number. Use the
-        //  mesh connectivity (if needed)
+        //- Decompose cells without weights. Code for weighted decomposition
+        //  is a bit complex so it is kept separate for now
+        virtual labelList decompose(const pointField&);
+
+        //- Decompose cells with weights
         virtual labelList decompose
         (
-            const pointField&,
+            const pointField& points,
             const scalarField& weights
         );
 
-        //- Without weights. Code for weighted decomposition is a bit complex
-        //  so kept separate for now.
-        virtual labelList decompose(const pointField&);
-
-        //- Return for every coordinate the wanted processor number. Explicitly
-        //  provided connectivity - does not use mesh_.
-        //  The connectivity is equal to mesh.cellCells() except for
-        //  - in parallel the cell numbers are global cell numbers (starting
-        //    from 0 at processor0 and then incrementing all through the
-        //    processors)
-        //  - the connections are across coupled patches
+        //- Decompose cells with weights with explicitly provided connectivity
+        //  Does not use mesh for connectivity
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
@@ -223,6 +219,7 @@ public:
             return decompose(cc, cWeights);
         }
 
+        //- Decompose cells with weights with explicitly provided connectivity
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
diff --git a/src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.C b/src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.C
index dc2beb2de..da7d53dd4 100644
--- a/src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.C
+++ b/src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.C
@@ -38,13 +38,6 @@ namespace Foam
 {
     defineTypeNameAndDebug(manualDecomp, 0);
 
-    addToRunTimeSelectionTable
-    (
-        decompositionMethod,
-        manualDecomp,
-        dictionary
-    );
-
     addToRunTimeSelectionTable
     (
         decompositionMethod,
@@ -56,14 +49,6 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::manualDecomp::manualDecomp(const dictionary& decompositionDict)
-:
-    decompositionMethod(decompositionDict)
-{
-    notImplemented("manualDecomp(const dictionary&)");
-}
-
-
 Foam::manualDecomp::manualDecomp
 (
     const dictionary& decompositionDict,
@@ -71,11 +56,14 @@ Foam::manualDecomp::manualDecomp
 )
 :
     decompositionMethod(decompositionDict),
-    meshPtr_(&mesh),
+    mesh_(mesh),
     decompDataFile_
     (
-        decompositionDict.subDict(word(decompositionDict.lookup("method"))
-      + "Coeffs").lookup("dataFile")
+        decompositionDict.subDict
+        (
+            word(decompositionDict.lookup("method"))
+          + "Coeffs"
+        ).lookup("dataFile")
     )
 {}
 
@@ -88,22 +76,20 @@ Foam::labelList Foam::manualDecomp::decompose
     const scalarField& pointWeights
 )
 {
-    const polyMesh& mesh = *meshPtr_;
-
     labelIOList finalDecomp
     (
         IOobject
         (
             decompDataFile_,
-            mesh.facesInstance(),
-            mesh,
+            mesh_.facesInstance(),
+            mesh_,
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE,
             false
         )
     );
 
-    // check if the final decomposition is OK
+    // Check if the final decomposition is OK
 
     if (finalDecomp.size() != points.size())
     {
diff --git a/src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.H b/src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.H
index 7cbe732ed..23df67442 100644
--- a/src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.H
+++ b/src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.H
@@ -42,7 +42,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class manualDecomp Declaration
+                        Class manualDecomp Declaration
 \*---------------------------------------------------------------------------*/
 
 class manualDecomp
@@ -51,17 +51,19 @@ class manualDecomp
 {
     // Private data
 
-        const polyMesh* meshPtr_;
+        const polyMesh& mesh_;
 
         fileName decompDataFile_;
 
 
     // Private Member Functions
 
-        //- Disallow default bitwise copy construct and assignment
-        void operator=(const manualDecomp&);
+        //- Disallow default bitwise copy construct
         manualDecomp(const manualDecomp&);
 
+        //- Disallow default bitwise assignment
+        void operator=(const manualDecomp&);
+
 
 public:
 
@@ -97,21 +99,14 @@ public:
             return true;
         }
 
-        //- Return for every coordinate the wanted processor number. Use the
-        //  mesh connectivity (if needed)
+        //- Decompose cells with weights
         virtual labelList decompose
         (
             const pointField& points,
             const scalarField& pointWeights
         );
 
-        //- Return for every coordinate the wanted processor number. Explicitly
-        //  provided connectivity - does not use mesh_.
-        //  The connectivity is equal to mesh.cellCells() except for
-        //  - in parallel the cell numbers are global cell numbers (starting
-        //    from 0 at processor0 and then incrementing all through the
-        //    processors)
-        //  - the connections are across coupled patches
+        //- Decompose cells with weights with explicitly provided connectivity
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
diff --git a/src/decompositionMethods/decompositionMethods/patchConstrainedDecomp/patchConstrainedDecomp.C b/src/decompositionMethods/decompositionMethods/patchConstrainedDecomp/patchConstrainedDecomp.C
new file mode 100644
index 000000000..11a01bae0
--- /dev/null
+++ b/src/decompositionMethods/decompositionMethods/patchConstrainedDecomp/patchConstrainedDecomp.C
@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Decomposition given a cell-to-processor association in a file
+
+\*---------------------------------------------------------------------------*/
+
+#include "patchConstrainedDecomp.H"
+#include "addToRunTimeSelectionTable.H"
+#include "IFstream.H"
+#include "labelIOList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(patchConstrainedDecomp, 0);
+
+    addToRunTimeSelectionTable
+    (
+        decompositionMethod,
+        patchConstrainedDecomp,
+        dictionaryMesh
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::patchConstrainedDecomp::patchConstrainedDecomp
+(
+    const dictionary& decompositionDict,
+    const polyMesh& mesh
+)
+:
+    decompositionMethod(decompositionDict),
+    mesh_(mesh),
+    dict_
+    (
+        decompositionDict.subDict
+        (
+            typeName + "Coeffs"
+        )
+    ),
+    baseDecompPtr_
+    (
+        decompositionMethod::New(dict_, mesh)
+    ),
+    patchConstraints_(dict_.lookup("patchConstraints"))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::labelList Foam::patchConstrainedDecomp::decompose
+(
+    const pointField& points,
+    const scalarField& pointWeights
+)
+{
+
+    labelList finalDecomp = baseDecompPtr_->decompose(points, pointWeights);
+
+    // Impose the decomposition along patches
+    forAll (patchConstraints_, i)
+    {
+        const label patchID =
+            mesh_.boundaryMesh().findPatchID(patchConstraints_[i].first());
+
+        const label procID = patchConstraints_[i].second();
+
+        if (patchID < 0 || procID < 0 || procID > nProcessors_ - 1)
+        {
+            FatalErrorIn
+            (
+                "labelList patchConstrainedDecomp::decompose\n"
+                "(\n"
+                "    const pointField& points,\n"
+                "    const scalarField& pointWeights\n"
+                ")"
+            )   << "Incorrect patch constraint definition for "
+                << patchConstraints_[i]
+                << abort(FatalError);
+        }
+
+        const labelList fc = mesh_.boundaryMesh()[patchID].faceCells();
+
+        forAll (fc, fcI)
+        {
+            finalDecomp[fc[fcI]] = procID;
+        }
+    }
+
+    return finalDecomp;
+}
+
+
+// ************************************************************************* //
diff --git a/src/decompositionMethods/decompositionMethods/patchConstrainedDecomp/patchConstrainedDecomp.H b/src/decompositionMethods/decompositionMethods/patchConstrainedDecomp/patchConstrainedDecomp.H
new file mode 100644
index 000000000..877f1ab4d
--- /dev/null
+++ b/src/decompositionMethods/decompositionMethods/patchConstrainedDecomp/patchConstrainedDecomp.H
@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::patchConstrainedDecomp
+
+Description
+    Base decomposition type is modified by assigning patches to given
+    processors.  This is controlled by the patch entry, giving patch names
+    and processor index for them.
+
+SourceFiles
+    patchConstrainedDecomp.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef patchConstrainedDecomp_H
+#define patchConstrainedDecomp_H
+
+#include "decompositionMethod.H"
+#include "Tuple2.H"
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class patchConstrainedDecomp Declaration
+\*---------------------------------------------------------------------------*/
+
+class patchConstrainedDecomp
+:
+    public decompositionMethod
+{
+    // Private data types
+
+    typedef List<Tuple2<word, label> > procWordList;
+
+
+    // Private data
+
+        //- Mesh
+        const polyMesh& mesh_;
+
+        //- Decomposition dictionary
+        dictionary dict_;
+
+        //- Base decomposition method
+        autoPtr<decompositionMethod> baseDecompPtr_;
+
+        //- Patch names to processor association
+        procWordList patchConstraints_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        patchConstrainedDecomp(const patchConstrainedDecomp&);
+
+        //- Disallow default bitwise copy construct and assignment
+        void operator=(const patchConstrainedDecomp&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("patchConstrained");
+
+
+    // Constructors
+
+        //- Construct given the decomposition dictionary and mesh
+        patchConstrainedDecomp
+        (
+            const dictionary& decompositionDict,
+            const polyMesh& mesh
+        );
+
+
+    // Destructor
+
+        virtual ~patchConstrainedDecomp()
+        {}
+
+
+    // Member Functions
+
+        //- Return parallel aware of base method
+        virtual bool parallelAware() const
+        {
+            return baseDecompPtr_->parallelAware();
+        }
+
+        //- Decompose cells with weights
+        virtual labelList decompose
+        (
+            const pointField& points,
+            const scalarField& pointWeights
+        );
+
+        //- Decompose cells with weights with explicitly provided connectivity
+        virtual labelList decompose
+        (
+            const labelListList& globalCellCells,
+            const pointField& cc,
+            const scalarField& cWeights
+        )
+        {
+            return decompose(cc, cWeights);
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C b/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
index d8c2d955b..e6b89871a 100644
--- a/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
+++ b/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
@@ -72,9 +72,9 @@ void Foam::simpleGeomDecomp::assignToProcessorGroup
 
     // assign cells to the first few processor groups (those with
     // one extra cell each
-    for (j=0; j<fstProcessorGroup; j++)
+    for (j = 0; j < fstProcessorGroup; j++)
     {
-        for (register label k=0; k<jumpb; k++)
+        for (register label k = 0; k < jumpb; k++)
         {
             processorGroup[ind++] = j;
         }
@@ -83,7 +83,7 @@ void Foam::simpleGeomDecomp::assignToProcessorGroup
     // and now to the `normal' processor groups
     for (; j<nProcGroup; j++)
     {
-        for (register label k=0; k<jump; k++)
+        for (register label k = 0; k < jump; k++)
         {
             processorGroup[ind++] = j;
         }
@@ -316,4 +316,6 @@ Foam::labelList Foam::simpleGeomDecomp::decompose
 
     return finalDecomp;
 }
+
+
 // ************************************************************************* //
diff --git a/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H b/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H
index abb2cf648..ed9ddb12a 100644
--- a/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H
+++ b/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H
@@ -41,7 +41,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class simpleGeomDecomp Declaration
+                      Class simpleGeomDecomp Declaration
 \*---------------------------------------------------------------------------*/
 
 class simpleGeomDecomp
@@ -61,10 +61,12 @@ class simpleGeomDecomp
             const scalar summedWeights
         );
 
-        //- Disallow default bitwise copy construct and assignment
-        void operator=(const simpleGeomDecomp&);
+        //- Disallow default bitwise copy construct
         simpleGeomDecomp(const simpleGeomDecomp&);
 
+        //- Disallow default bitwise assignment
+        void operator=(const simpleGeomDecomp&);
+
 
 public:
 
@@ -75,7 +77,7 @@ public:
     // Constructors
 
         //- Construct given the decomposition dictionary
-        simpleGeomDecomp(const dictionary& decompositionDict);
+        explicit simpleGeomDecomp(const dictionary& decompositionDict);
 
         //- Construct given the decomposition dictionary and mesh
         simpleGeomDecomp