From 8ca9d01c74aadc8cb5196ce385d44ebff8bdfc94 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Thu, 12 Apr 2018 17:47:18 +0100
Subject: [PATCH] Modified decomposition tools: passive patches

---
 .../decomposeReconstruct/Make/files           |   3 +
 .../finiteVolume/decomposeMesh.C              | 359 ++++++++++++++----
 .../finiteVolume/distributeCells.C            |  43 ++-
 .../finiteVolume/domainDecomposition.C        |  72 ++--
 .../finiteVolume/domainDecomposition.H        |  27 +-
 5 files changed, 392 insertions(+), 112 deletions(-)

diff --git a/src/decompositionMethods/decomposeReconstruct/Make/files b/src/decompositionMethods/decomposeReconstruct/Make/files
index 061c4f182..0d3f4f821 100644
--- a/src/decompositionMethods/decomposeReconstruct/Make/files
+++ b/src/decompositionMethods/decomposeReconstruct/Make/files
@@ -1,3 +1,6 @@
+passivePatches/passiveProcessorPolyPatch/passiveProcessorPolyPatch.C
+passivePatches/passiveProcessorFvPatch/passiveProcessorFvPatch.C
+
 decomposeTools/finiteVolume/domainDecomposition.C
 decomposeTools/finiteVolume/distributeCells.C
 decomposeTools/finiteVolume/decomposeMesh.C
diff --git a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/decomposeMesh.C b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/decomposeMesh.C
index 189c8c558..56bbea30d 100644
--- a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/decomposeMesh.C
+++ b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/decomposeMesh.C
@@ -36,7 +36,8 @@ Description
 #include "boolList.H"
 #include "cellList.H"
 #include "primitiveMesh.H"
-#include "cyclicPolyPatch.H"
+#include "processorFvPatch.H"
+#include "cyclicFvPatch.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -50,11 +51,11 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
     // calculate the addressing information for the original mesh
     Info<< "\nCalculating original mesh data" << endl;
 
-    // set references to the original mesh
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    // Set references to the original mesh
+    const fvBoundaryMesh& patches = mesh_.boundary();
 
     // Access all faces to grab the zones
-    const faceList& fcs = mesh_.allFaces();
+    const faceList& allFaces = mesh_.allFaces();
     const labelList& owner = mesh_.faceOwner();
     const labelList& neighbour = mesh_.faceNeighbour();
 
@@ -100,12 +101,12 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
     {
         List<SLList<label> > procFaceList(nProcs_);
 
-        forAll (neighbour, facei)
+        forAll (neighbour, faceI)
         {
-            if (cellToProc_[owner[facei]] == cellToProc_[neighbour[facei]])
+            if (cellToProc_[owner[faceI]] == cellToProc_[neighbour[faceI]])
             {
                 // Face internal to processor
-                procFaceList[cellToProc_[owner[facei]]].append(facei);
+                procFaceList[cellToProc_[owner[faceI]]].append(faceI);
             }
         }
 
@@ -125,17 +126,19 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
 
         List<SLList<label> > procPatchIndex(nProcs_);
 
-        forAll (neighbour, facei)
+        forAll (neighbour, faceI)
         {
-            if (cellToProc_[owner[facei]] != cellToProc_[neighbour[facei]])
+            if (cellToProc_[owner[faceI]] != cellToProc_[neighbour[faceI]])
             {
-                // inter - processor patch face found. Go through the list of
+                // Inter - processor patch face found. Go through the list of
                 // inside boundaries for the owner processor and try to find
                 // this inter-processor patch.
 
-                label ownerProc = cellToProc_[owner[facei]];
-                label neighbourProc = cellToProc_[neighbour[facei]];
+                const label ownerProc = cellToProc_[owner[faceI]];
+                const label neighbourProc = cellToProc_[neighbour[faceI]];
 
+                // Search algorithm repeated in processor patches.  Reconsider
+                // HJ, 11/Apr/2018
                 SLList<label>::iterator curInterProcBdrsOwnIter =
                     interProcBoundaries[ownerProc].begin();
 
@@ -158,10 +161,10 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                 {
                     if (curInterProcBdrsOwnIter() == neighbourProc)
                     {
-                        // the inter - processor boundary exists. Add the face
+                        // Inter - processor boundary exists. Add the face
                         interProcBouFound = true;
 
-                        curInterProcBFacesOwnIter().append(facei);
+                        curInterProcBFacesOwnIter().append(faceI);
 
                         SLList<label>::iterator curInterProcBdrsNeiIter =
                             interProcBoundaries[neighbourProc].begin();
@@ -190,7 +193,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                                 // boundary found. Add the face
                                 neighbourFound = true;
 
-                                curInterProcBFacesNeiIter().append(facei);
+                                curInterProcBFacesNeiIter().append(faceI);
                             }
 
                             if (neighbourFound) break;
@@ -220,66 +223,264 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
 
                     // owner
                     interProcBoundaries[ownerProc].append(neighbourProc);
-                    interProcBFaces[ownerProc].append(SLList<label>(facei));
+                    interProcBFaces[ownerProc].append(SLList<label>(faceI));
 
                     // neighbour
                     interProcBoundaries[neighbourProc].append(ownerProc);
                     interProcBFaces[neighbourProc].append
                     (
-                        SLList<label>(facei)
+                        SLList<label>(faceI)
                     );
                 }
             }
         }
 
+        // Rewrite:
+        // Handling of coupled patches is changed.  HJ, 11/Apr/2018
+
+        // Prepare collection of patch faces
+
+        // For all processors, set the size of start index and patch size
+        // lists to the number of patches in the mesh
+        forAll (procPatchSize_, procI)
+        {
+            procPatchSize_[procI].setSize(patches.size());
+            procPatchStartIndex_[procI].setSize(patches.size());
+        }
+
+        forAll (patches, patchI)
+        {
+            // Reset size and start index for all processors
+            forAll (procPatchSize_, procI)
+            {
+                procPatchSize_[procI][patchI] = 0;
+                procPatchStartIndex_[procI][patchI] =
+                    procFaceList[procI].size();
+            }
+        }
+
+
+        // Algorithm:
+        // When running the decomposition in parallel, it is assumed that
+        // the result will be used in dynamic load balancing
+        // For load balancing, the processor patches need to match in order for
+        // the patch-to-patch matching to work propely on mesh reconstruction
+        // Therefore, all processor patches need to be split into the matching
+        // and non-matching part be examining the cellToProc_ data on
+        // the neighbour side.  This has been prepared in patchNbrCellToProc_
+        // HJ, 11/Apr/2018
+
+        forAll (patches, patchI)
+        {
+            // Check the processor patch for which neighbour data exists
+            if
+            (
+                isA<processorFvPatch>(patches[patchI])
+             && !patchNbrCellToProc_[patchI].empty()
+            )
+            {
+                // Get patch start
+                const label patchStart = patches[patchI].patch().start();
+
+                // Get faceCells
+                const labelList& fc = patches[patchI].faceCells();
+
+                // Get neighbour cellToProc addressing across the interface
+                const labelList& curNbrPtc = patchNbrCellToProc_[patchI];
+
+                forAll (fc, patchFaceI)
+                {
+                    // Local owner proc is looked up using faceCells
+                    const label ownerProc = cellToProc_[fc[patchFaceI]];
+
+                    // Neighbour proc is looked up directly
+                    const label neighbourProc = curNbrPtc[patchFaceI];
+
+                    // Check change in processor type across the processor
+                    // boundary
+                    // If ownerProc and neighbourProc are the same, the
+                    // processor face will be merged, meaning that it remains
+                    // in the (old) processor patch
+                    // If ownerProc and neighbourProc are different, this will
+                    // be a new processor boundary created from the existing
+                    // processor face and added afterwards
+                    if (ownerProc != neighbourProc)
+                    {
+                        Pout<< "New proc from old proc[" << patchFaceI << "]: "
+                            << ownerProc << " " << neighbourProc << endl;
+                        // Search algorithm repeated in processor patches.
+                        // HJ, 11/Apr/2018
+                        SLList<label>::iterator curInterProcBdrsOwnIter =
+                            interProcBoundaries[ownerProc].begin();
+
+                        SLList<SLList<label> >::iterator
+                            curInterProcBFacesOwnIter =
+                            interProcBFaces[ownerProc].begin();
+
+                        bool interProcBouFound = false;
+
+                        // WARNING: Synchronous SLList iterators
+
+                        for
+                        (
+                            ;
+                            curInterProcBdrsOwnIter
+                         != interProcBoundaries[ownerProc].end()
+                         && curInterProcBFacesOwnIter
+                         != interProcBFaces[ownerProc].end();
+                            ++curInterProcBdrsOwnIter,
+                            ++curInterProcBFacesOwnIter
+                        )
+                        {
+                            if (curInterProcBdrsOwnIter() == neighbourProc)
+                            {
+                                // Inter - processor boundary exists.
+                                // Add the face
+                                interProcBouFound = true;
+
+                                curInterProcBFacesOwnIter().append(patchStart + patchFaceI);
+
+                                SLList<label>::iterator
+                                    curInterProcBdrsNeiIter =
+                                    interProcBoundaries[neighbourProc].begin();
+
+                                SLList<SLList<label> >::iterator
+                                    curInterProcBFacesNeiIter =
+                                    interProcBFaces[neighbourProc].begin();
+
+                                bool neighbourFound = false;
+
+                                // WARNING: Synchronous SLList iterators
+
+                                for
+                                (
+                                    ;
+                                    curInterProcBdrsNeiIter !=
+                                    interProcBoundaries[neighbourProc].end()
+                                 && curInterProcBFacesNeiIter !=
+                                    interProcBFaces[neighbourProc].end();
+                                    ++curInterProcBdrsNeiIter,
+                                    ++curInterProcBFacesNeiIter
+                                )
+                                {
+                                    if (curInterProcBdrsNeiIter() == ownerProc)
+                                    {
+                                        // Boundary found. Add the face
+                                        neighbourFound = true;
+
+                                        curInterProcBFacesNeiIter().append
+                                        (
+                                            patchStart + patchFaceI
+                                        );
+                                    }
+
+                                    if (neighbourFound) break;
+                                }
+
+                                if (interProcBouFound && !neighbourFound)
+                                {
+                                    FatalErrorIn
+                                    (
+                                        "domainDecomposition::decomposeMesh()"
+                                    )   << "Inconsistency in inter - "
+                                        << "processor boundary lists for "
+                                        << "processors "
+                                        << ownerProc << " and "
+                                        << neighbourProc
+                                        << abort(FatalError);
+                                }
+                            }
+
+                            if (interProcBouFound) break;
+                        }
+
+                        if (!interProcBouFound)
+                        {
+                            // inter - processor boundaries do not exist and
+                            // need to be created
+
+                            // set the new addressing information
+
+                            // owner
+                            interProcBoundaries[ownerProc].append
+                            (
+                                neighbourProc
+                            );
+                            interProcBFaces[ownerProc].append
+                            (
+                                SLList<label>(patchStart + patchFaceI)
+                            );
+
+                            // neighbour
+                            interProcBoundaries[neighbourProc].append
+                            (
+                                ownerProc
+                            );
+                            interProcBFaces[neighbourProc].append
+                            (
+                                SLList<label>(patchStart + patchFaceI)
+                            );
+                        }
+                    }
+                    else
+                    {
+                        // Owner and neighbour processor index are the same
+                        // across the processor boundary.  This face will be
+                        // re-merged into internal faces in load balancing
+                        // and therefore remains in the processor patch
+                        Pout<< "Preserved proc[" << patchFaceI << "]: "
+                            << ownerProc << " " << neighbourProc << endl;
+
+                        // add the face
+                        procFaceList[ownerProc].append(patchStart + patchFaceI);
+
+                        // increment the number of faces for this patch
+                        procPatchSize_[ownerProc][patchI]++;
+                    }
+                }
+            }
+        }
+
         // Loop through patches. For cyclic boundaries detect inter-processor
         // faces; for all other, add faces to the face list and remember start
         // and size of all patches.
 
-        // for all processors, set the size of start index and patch size
-        // lists to the number of patches in the mesh
-        forAll (procPatchSize_, procI)
+        forAll (patches, patchI)
         {
-            procPatchSize_[procI].setSize(patches.size());
-            procPatchStartIndex_[procI].setSize(patches.size());
-        }
+            const label patchStart = patches[patchI].patch().start();
 
-        forAll (patches, patchi)
-        {
-            // Reset size and start index for all processors
-            forAll (procPatchSize_, procI)
-            {
-                procPatchSize_[procI][patchi] = 0;
-                procPatchStartIndex_[procI][patchi] =
-                    procFaceList[procI].size();
-            }
-
-            const label patchStart = patches[patchi].start();
-
-            if (!isA<cyclicPolyPatch>(patches[patchi]))
+            // Do normal patches.  Note that processor patches
+            // have already been done and need to be skipped
+            // For all other patches patchNbrCellToProc_ will be empty
+            if
+            (
+                !isA<cyclicFvPatch>(patches[patchI])
+             && patchNbrCellToProc_[patchI].empty()
+            )
             {
                 // Normal patch. Add faces to processor where the cell
                 // next to the face lives
 
                 const unallocLabelList& patchFaceCells =
-                    patches[patchi].faceCells();
+                    patches[patchI].faceCells();
 
-                forAll (patchFaceCells, facei)
+                forAll (patchFaceCells, patchFaceI)
                 {
-                    const label curProc = cellToProc_[patchFaceCells[facei]];
+                    const label curProc =
+                        cellToProc_[patchFaceCells[patchFaceI]];
 
                     // add the face
-                    procFaceList[curProc].append(patchStart + facei);
+                    procFaceList[curProc].append(patchStart + patchFaceI);
 
                     // increment the number of faces for this patch
-                    procPatchSize_[curProc][patchi]++;
+                    procPatchSize_[curProc][patchI]++;
                 }
             }
             else
             {
                 // Cyclic patch special treatment
 
-                const polyPatch& cPatch = patches[patchi];
+                const fvPatch& cPatch = patches[patchI];
 
                 const label cycOffset = cPatch.size()/2;
 
@@ -297,12 +498,12 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                     cycOffset
                 );
 
-                forAll (firstFaceCells, facei)
+                forAll (firstFaceCells, patchFaceI)
                 {
                     if
                     (
-                        cellToProc_[firstFaceCells[facei]]
-                     != cellToProc_[secondFaceCells[facei]]
+                        cellToProc_[firstFaceCells[patchFaceI]]
+                     != cellToProc_[secondFaceCells[patchFaceI]]
                     )
                     {
                         // This face becomes an inter-processor boundary face
@@ -313,9 +514,11 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
 
                         cyclicParallel_ = true;
 
-                        label ownerProc = cellToProc_[firstFaceCells[facei]];
+                        label ownerProc =
+                            cellToProc_[firstFaceCells[patchFaceI]];
+
                         label neighbourProc =
-                            cellToProc_[secondFaceCells[facei]];
+                            cellToProc_[secondFaceCells[patchFaceI]];
 
                         SLList<label>::iterator curInterProcBdrsOwnIter =
                             interProcBoundaries[ownerProc].begin();
@@ -346,7 +549,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                                 interProcBouFound = true;
 
                                 curInterProcBFacesOwnIter().append
-                                    (patchStart + facei);
+                                    (patchStart + patchFaceI);
 
                                 SLList<label>::iterator curInterProcBdrsNeiIter
                                    = interProcBoundaries[neighbourProc].begin();
@@ -380,7 +583,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                                             (
                                                 patchStart
                                               + cycOffset
-                                              + facei
+                                              + patchFaceI
                                             );
                                     }
 
@@ -415,7 +618,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                             interProcBoundaries[ownerProc]
                                 .append(neighbourProc);
                             interProcBFaces[ownerProc]
-                                .append(SLList<label>(patchStart + facei));
+                                .append(SLList<label>(patchStart + patchFaceI));
 
                             // neighbour
                             interProcBoundaries[neighbourProc]
@@ -427,7 +630,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                                     (
                                         patchStart
                                       + cycOffset
-                                      + facei
+                                      + patchFaceI
                                     )
                                 );
                         }
@@ -435,13 +638,14 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                     else
                     {
                         // This cyclic face remains on the processor
-                        label ownerProc = cellToProc_[firstFaceCells[facei]];
+                        label ownerProc =
+                            cellToProc_[firstFaceCells[patchFaceI]];
 
                         // add the face
-                        procFaceList[ownerProc].append(patchStart + facei);
+                        procFaceList[ownerProc].append(patchStart + patchFaceI);
 
                         // increment the number of faces for this patch
-                        procPatchSize_[ownerProc][patchi]++;
+                        procPatchSize_[ownerProc][patchI]++;
 
                         // Note: I cannot add the other side of the cyclic
                         // boundary here because this would violate the order.
@@ -453,23 +657,24 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                 // Ordering in cyclic boundaries is important.
                 // Add the other half of cyclic faces for cyclic boundaries
                 // that remain on the processor
-                forAll (secondFaceCells, facei)
+                forAll (secondFaceCells, patchFaceI)
                 {
                     if
                     (
-                        cellToProc_[firstFaceCells[facei]]
-                     == cellToProc_[secondFaceCells[facei]]
+                        cellToProc_[firstFaceCells[patchFaceI]]
+                     == cellToProc_[secondFaceCells[patchFaceI]]
                     )
                     {
                         // This cyclic face remains on the processor
-                        label ownerProc = cellToProc_[firstFaceCells[facei]];
+                        label ownerProc =
+                            cellToProc_[firstFaceCells[patchFaceI]];
 
-                        // add the second face
+                        // Add the second face
                         procFaceList[ownerProc].append
-                            (patchStart + cycOffset + facei);
+                            (patchStart + cycOffset + patchFaceI);
 
-                        // increment the number of faces for this patch
-                        procPatchSize_[ownerProc][patchi]++;
+                        // Increment the number of faces for this patch
+                        procPatchSize_[ownerProc][patchI]++;
                     }
                 }
             }
@@ -557,7 +762,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
             labelList& curProcProcessorPatchStartIndex =
                 procProcessorPatchStartIndex_[procI];
 
-            // calculate the size
+            // Calculate the size
             label nFacesOnProcessor = curProcFaces.size();
 
             for
@@ -727,15 +932,15 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
 
         label nPatches = 0;
 
-        forAll (oldPatchSizes, patchi)
+        forAll (oldPatchSizes, patchI)
         {
-            if (!filterEmptyPatches || oldPatchSizes[patchi] > 0)
+            if (!filterEmptyPatches || oldPatchSizes[patchI] > 0)
             {
-                curBoundaryAddressing[nPatches] = patchi;
+                curBoundaryAddressing[nPatches] = patchI;
 
-                curPatchSizes[nPatches] = oldPatchSizes[patchi];
+                curPatchSizes[nPatches] = oldPatchSizes[patchI];
 
-                curPatchStarts[nPatches] = oldPatchStarts[patchi];
+                curPatchStarts[nPatches] = oldPatchStarts[patchI];
 
                 nPatches++;
             }
@@ -788,7 +993,8 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
         for (label faceI = 0; faceI < nLiveProcFaces_[procI]; faceI++)
         {
             // Because of the turning index, some labels may be negative
-            const labelList& facePoints = fcs[mag(procFaceLabels[faceI]) - 1];
+            const labelList& facePoints =
+                allFaces[mag(procFaceLabels[faceI]) - 1];
 
             forAll (facePoints, pointI)
             {
@@ -823,7 +1029,8 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
         )
         {
             // Because of the turning index, some labels may be negative
-            const labelList& facePoints = fcs[mag(procFaceLabels[faceI]) - 1];
+            const labelList& facePoints =
+                allFaces[mag(procFaceLabels[faceI]) - 1];
 
             forAll (facePoints, pointI)
             {
@@ -882,10 +1089,10 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
             // Reset the lookup list
             pointsUsage = 0;
 
-            forAll (curProcessorPatchStarts, patchi)
+            forAll (curProcessorPatchStarts, patchI)
             {
-                const label curStart = curProcessorPatchStarts[patchi];
-                const label curEnd = curStart + curProcessorPatchSizes[patchi];
+                const label curStart = curProcessorPatchStarts[patchI];
+                const label curEnd = curStart + curProcessorPatchSizes[patchI];
 
                 for
                 (
@@ -899,16 +1106,16 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                     // HJ, 5/Dec/2001
                     const label curF = mag(curFaceLabels[faceI]) - 1;
 
-                    const face& f = fcs[curF];
+                    const face& f = allFaces[curF];
 
                     forAll (f, pointI)
                     {
                         if (pointsUsage[f[pointI]] == 0)
                         {
                             // Point not previously used
-                            pointsUsage[f[pointI]] = patchi + 1;
+                            pointsUsage[f[pointI]] = patchI + 1;
                         }
-                        else if (pointsUsage[f[pointI]] != patchi + 1)
+                        else if (pointsUsage[f[pointI]] != patchI + 1)
                         {
                             // Point used by some other patch = global point!
                             gSharedPoints.insert(f[pointI]);
diff --git a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/distributeCells.C b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/distributeCells.C
index 4503bcf2b..8fc092847 100644
--- a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/distributeCells.C
+++ b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/distributeCells.C
@@ -26,7 +26,7 @@ License
 #include "domainDecomposition.H"
 #include "decompositionMethod.H"
 #include "cpuTime.H"
-#include "cyclicPolyPatch.H"
+#include "processorFvPatch.H"
 #include "cellSet.H"
 #include "regionSplit.H"
 
@@ -38,7 +38,6 @@ void Foam::domainDecomposition::distributeCells()
 
     cpuTime decompositionTime;
 
-
     // See if any faces need to have owner and neighbour on same processor
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -179,6 +178,46 @@ void Foam::domainDecomposition::distributeCells()
         cellToProc_ = decomposePtr().decompose(globalRegion, regionCentres);
     }
 
+    // If running in parallel, sync cellToProc_ across coupled boundaries
+    // Initialise transfer of restrict addressing on the interface
+    if (Pstream::parRun())
+    {
+        const fvBoundaryMesh& patches = mesh_.boundary();
+
+        forAll (patches, patchI)
+        {
+            if (isA<processorFvPatch>(patches[patchI]))
+            // if (patches[patchI].coupled())
+            {
+                const lduInterface& cpPatch =
+                    refCast<const lduInterface>(patches[patchI]);
+
+                cpPatch.initInternalFieldTransfer
+                (
+                    Pstream::blocking,
+                    cellToProc_
+                );
+            }
+        }
+
+        forAll (patches, patchI)
+        {
+            if (isA<processorFvPatch>(patches[patchI]))
+            // if (patches[patchI].coupled())
+            {
+                const lduInterface& cpPatch =
+                    refCast<const lduInterface>(patches[patchI]);
+
+                patchNbrCellToProc_[patchI] =
+                    cpPatch.internalFieldTransfer
+                    (
+                        Pstream::blocking,
+                        cellToProc_
+                    );
+            }
+        }
+    }
+    
     Info<< "\nFinished decomposition in "
         << decompositionTime.elapsedCpuTime()
         << " s" << endl;
diff --git a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.C b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.C
index 6dae788ed..c43246b63 100644
--- a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.C
+++ b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.C
@@ -28,6 +28,7 @@ License
 #include "dictionary.H"
 #include "labelIOList.H"
 #include "processorPolyPatch.H"
+#include "passiveProcessorPolyPatch.H"
 #include "fvMesh.H"
 #include "OSspecific.H"
 #include "Map.H"
@@ -49,8 +50,8 @@ Foam::domainDecomposition::domainDecomposition
     mesh_(mesh),
     decompositionDict_(dict),
     nProcs_(readInt(decompositionDict_.lookup("numberOfSubdomains"))),
-    distributed_(false),
     cellToProc_(mesh_.nCells()),
+    patchNbrCellToProc_(mesh_.boundaryMesh().size()),
     procPointAddressing_(nProcs_),
     procFaceAddressing_(nProcs_),
     nInternalProcFaces_(nProcs_),
@@ -64,12 +65,7 @@ Foam::domainDecomposition::domainDecomposition
     procProcessorPatchStartIndex_(nProcs_),
     globallySharedPoints_(0),
     cyclicParallel_(false)
-{
-    if (decompositionDict_.found("distributed"))
-    {
-        distributed_ = Switch(decompositionDict_.lookup("distributed"));
-    }
-}
+{}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -84,7 +80,8 @@ Foam::autoPtr<Foam::fvMesh> Foam::domainDecomposition::processorMesh
 (
     const label procI,
     const Time& processorDb,
-    const word& regionName
+    const word& regionName,
+    const bool createPassiveProcPatches
 ) const
 {
     // Create processor points
@@ -265,27 +262,54 @@ Foam::autoPtr<Foam::fvMesh> Foam::domainDecomposition::processorMesh
         nPatches++;
     }
 
-    forAll (curProcessorPatchSizes, procPatchI)
+    if (createPassiveProcPatches)
     {
-        procPatches[nPatches] =
-            new processorPolyPatch
-            (
-                word("procBoundary") + Foam::name(procI)
-              + word("to")
-              + Foam::name(curNeighbourProcessors[procPatchI]),
-                curProcessorPatchSizes[procPatchI],
-                curProcessorPatchStarts[procPatchI],
-                nPatches,
-                procMesh.boundaryMesh(),
-                procI,
-                curNeighbourProcessors[procPatchI]
-        );
+        forAll (curProcessorPatchSizes, procPatchI)
+        {
+            procPatches[nPatches] =
+                new passiveProcessorPolyPatch
+                (
+                    word("procBoundary") + Foam::name(procI)
+                  + word("to")
+                  + Foam::name(curNeighbourProcessors[procPatchI]),
+                    curProcessorPatchSizes[procPatchI],
+                    curProcessorPatchStarts[procPatchI],
+                    nPatches,
+                    procMesh.boundaryMesh(),
+                    procI,
+                    curNeighbourProcessors[procPatchI]
+                );
 
-        nPatches++;
+            nPatches++;
+        }
+    }
+    else
+    {
+        forAll (curProcessorPatchSizes, procPatchI)
+        {
+            procPatches[nPatches] =
+                new processorPolyPatch
+                (
+                    word("procBoundary") + Foam::name(procI)
+                  + word("to")
+                  + Foam::name(curNeighbourProcessors[procPatchI]),
+                    curProcessorPatchSizes[procPatchI],
+                    curProcessorPatchStarts[procPatchI],
+                    nPatches,
+                    procMesh.boundaryMesh(),
+                    procI,
+                    curNeighbourProcessors[procPatchI]
+                );
+
+            nPatches++;
+        }
     }
 
+
     // Add boundary patches to polyMesh and fvMesh
-    procMesh.addFvPatches(procPatches);
+    // Note: Mark boundary as invalid to disable analysis
+    // due to the presence of old/new patches
+    procMesh.addFvPatches(procPatches, false);
 
     // Create and add zones
 
diff --git a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.H b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.H
index f31e7bc43..5b3e42a39 100644
--- a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.H
+++ b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.H
@@ -65,12 +65,13 @@ class domainDecomposition
         //- Number of processors in decomposition
         label nProcs_;
 
-        //- Is the decomposition data to be distributed for each processor
-        bool distributed_;
-
         //- Processor label for each cell
         labelList cellToProc_;
 
+        //- Processor label for neighbour cell for each processor boundary
+        //  Data is used when running decomposition in parallel
+        labelListList patchNbrCellToProc_;
+
         //- Labels of points for each processor
         labelListList procPointAddressing_;
 
@@ -159,18 +160,19 @@ public:
             return nProcs_;
         }
 
-        //- Is the decomposition data to be distributed for each processor
-        bool distributed() const
-        {
-            return distributed_;
-        }
-
         //- Return cell-processor decomposition labels
         const labelList& cellToProc() const
         {
             return cellToProc_;
         }
 
+        //- Return cell-processor decomposition labels for cells across
+        //  coupled boundaries when running decomposition in parallel
+        const labelListList& patchNbrCellToProc() const
+        {
+            return patchNbrCellToProc_;
+        }
+
         //- Decompose mesh. Optionally remove zero-sized patches.
         void decomposeMesh(const bool filterEmptyPatches);
 
@@ -178,11 +180,16 @@ public:
         // Decomposed mesh and addressing
 
             //- Create a decomposed mesh for a given processor index
+            //  Note: at the point of construction, the boundary is marked
+            //  as invalid.  If the mesh should be used immediately upon
+            //  creation, initialise the boundary patches before use
+            //  HJ, 12/Apr/2018
             autoPtr<fvMesh> processorMesh
             (
                 const label procI,
                 const Time& procDb,
-                const word& regionName
+                const word& regionName,
+                const bool createPassiveProcPatches = false
             ) const;
 
             //- Return processor point addressing