Load balancing development, 5

src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/decomposeMesh.C

@@ -250,12 +250,10 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
 
         forAll (patches, patchI)
         {
-            // Reset size and start index for all processors
+            // Reset size for all processors
             forAll (procPatchSize_, procI)
             {
                 procPatchSize_[procI][patchI] = 0;
-                procPatchStartIndex_[procI][patchI] =
-                    procFaceList[procI].size();
             }
         }
 
@@ -270,6 +268,8 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
         // the neighbour side.  This has been prepared in patchNbrCellToProc_
         // HJ, 11/Apr/2018
 
+        // Dump current processor patch faces into new processor patches
+        // that will be created in decomposition when running in parallel
         forAll (patches, patchI)
         {
             // Check the processor patch for which neighbour data exists
@@ -338,7 +338,10 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                                 // Add the face
                                 interProcBouFound = true;
 
-                                curInterProcBFacesOwnIter().append(patchStart + patchFaceI);
+                                curInterProcBFacesOwnIter().append
+                                (
+                                    patchStart + patchFaceI
+                                );
 
                                 SLList<label>::iterator
                                     curInterProcBdrsNeiIter =
@@ -447,6 +450,13 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
 
         forAll (patches, patchI)
         {
+            // New patch: record start index for all processors
+            forAll (procPatchSize_, procI)
+            {
+                procPatchStartIndex_[procI][patchI] =
+                    procFaceList[procI].size();
+            }
+
             const label patchStart = patches[patchI].patch().start();
 
             // Do normal patches.  Note that processor patches

src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.C

@@ -50,6 +50,7 @@ Foam::domainDecomposition::domainDecomposition
     mesh_(mesh),
     decompositionDict_(dict),
     nProcs_(readInt(decompositionDict_.lookup("numberOfSubdomains"))),
+    distributed_(false),
     cellToProc_(mesh_.nCells()),
     patchNbrCellToProc_(mesh_.boundaryMesh().size()),
     procPointAddressing_(nProcs_),
@@ -65,7 +66,12 @@ Foam::domainDecomposition::domainDecomposition
     procProcessorPatchStartIndex_(nProcs_),
     globallySharedPoints_(0),
     cyclicParallel_(false)
-{}
+{
+    if (decompositionDict_.found("distributed"))
+    {
+        distributed_ = Switch(decompositionDict_.lookup("distributed"));
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.H

@@ -65,6 +65,9 @@ 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_;
 
@@ -160,6 +163,12 @@ public:
             return nProcs_;
         }
 
+        //- Is the decomposition data to be distributed for each processor
+        inline bool distributed() const
+        {
+            return distributed_;
+        }
+
         //- Return cell-processor decomposition labels
         const labelList& cellToProc() const
         {

src/decompositionMethods/decomposeReconstruct/passivePatches/passiveProcessorFvPatch/passiveProcessorFvPatch.C

@@ -0,0 +1,40 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "passiveProcessorFvPatch.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(passiveProcessorFvPatch, 0);
+addToRunTimeSelectionTable(fvPatch, passiveProcessorFvPatch, polyPatch);
+
+}
+
+
+// ************************************************************************* //

src/decompositionMethods/decomposeReconstruct/passivePatches/passiveProcessorFvPatch/passiveProcessorFvPatch.H

@@ -0,0 +1,94 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::passiveProcessorFvPatch
+
+Description
+    Passive processor fvPatch.  See passiveProcessorPolyPatch
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
+
+SourceFiles
+    passiveProcessorFvPatch.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef passiveProcessorFvPatch_H
+#define passiveProcessorFvPatch_H
+
+#include "coupledFvPatch.H"
+#include "processorLduInterface.H"
+#include "passiveProcessorPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class passiveProcessorFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class passiveProcessorFvPatch
+:
+    public fvPatch
+{
+public:
+
+    //- Runtime type information
+    TypeName(passiveProcessorPolyPatch::typeName_());
+
+
+    // Constructors
+
+        //- Construct from components
+        passiveProcessorFvPatch
+        (
+            const polyPatch& patch,
+            const fvBoundaryMesh& bm
+        )
+        :
+            fvPatch(patch, bm)
+        {}
+
+
+    //- Destructor
+    virtual ~passiveProcessorFvPatch()
+    {}
+
+
+    // Member functions
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/decompositionMethods/decomposeReconstruct/passivePatches/passiveProcessorPolyPatch/passiveProcessorPolyPatch.C

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "passiveProcessorPolyPatch.H"
+#include "addToRunTimeSelectionTable.H"
+#include "dictionary.H"
+#include "SubField.H"
+#include "matchPoints.H"
+#include "OFstream.H"
+#include "polyBoundaryMesh.H"
+#include "polyMesh.H"
+#include "foamTime.H"
+#include "transformList.H"
+#include "demandDrivenData.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(passiveProcessorPolyPatch, 0);
+
+    addToRunTimeSelectionTable
+    (
+        polyPatch,
+        passiveProcessorPolyPatch,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::passiveProcessorPolyPatch::passiveProcessorPolyPatch
+(
+    const word& name,
+    const label size,
+    const label start,
+    const label index,
+    const polyBoundaryMesh& bm,
+    const int myProcNo,
+    const int neighbProcNo
+)
+:
+    polyPatch(name, size, start, index, bm),
+    myProcNo_(myProcNo),
+    neighbProcNo_(neighbProcNo)
+{}
+
+
+Foam::passiveProcessorPolyPatch::passiveProcessorPolyPatch
+(
+    const word& name,
+    const dictionary& dict,
+    const label index,
+    const polyBoundaryMesh& bm
+)
+:
+    polyPatch(name, dict, index, bm),
+    myProcNo_(readLabel(dict.lookup("myProcNo"))),
+    neighbProcNo_(readLabel(dict.lookup("neighbProcNo")))
+{}
+
+
+Foam::passiveProcessorPolyPatch::passiveProcessorPolyPatch
+(
+    const passiveProcessorPolyPatch& pp,
+    const polyBoundaryMesh& bm
+)
+:
+    polyPatch(pp, bm),
+    myProcNo_(pp.myProcNo_),
+    neighbProcNo_(pp.neighbProcNo_)
+{}
+
+
+Foam::passiveProcessorPolyPatch::passiveProcessorPolyPatch
+(
+    const passiveProcessorPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const label newSize,
+    const label newStart
+)
+:
+    polyPatch(pp, bm, index, newSize, newStart),
+    myProcNo_(pp.myProcNo_),
+    neighbProcNo_(pp.neighbProcNo_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::passiveProcessorPolyPatch::~passiveProcessorPolyPatch()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::passiveProcessorPolyPatch::write(Ostream& os) const
+{
+    polyPatch::write(os);
+    os.writeKeyword("myProcNo") << myProcNo_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("neighbProcNo") << neighbProcNo_
+        << token::END_STATEMENT << nl;
+}
+
+
+// ************************************************************************* //

src/decompositionMethods/decomposeReconstruct/passivePatches/passiveProcessorPolyPatch/passiveProcessorPolyPatch.H

@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::passiveProcessorPolyPatch
+
+Description
+    Passive processor polyPatch.
+
+    Passive processor polyPatch is created in load balancing when new processor
+    faces are built on each parallel decomposition.  It contains the data
+    which is needed to rebuild correct processor boundaries after load balancing
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
+
+SourceFiles
+    passiveProcessorPolyPatch.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef passiveProcessorPolyPatch_H
+#define passiveProcessorPolyPatch_H
+
+#include "polyBoundaryMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class passiveProcessorPolyPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class passiveProcessorPolyPatch
+:
+    public polyPatch
+{
+    // Private data
+
+        //- My processor number
+        int myProcNo_;
+
+        //- Neighbour processor number
+        int neighbProcNo_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("passiveProcessor");
+
+
+    // Constructors
+
+        //- Construct from components
+        passiveProcessorPolyPatch
+        (
+            const word& name,
+            const label size,
+            const label start,
+            const label index,
+            const polyBoundaryMesh& bm,
+            const int myProcNo,
+            const int neighbProcNo
+        );
+
+        //- Construct from dictionary
+        passiveProcessorPolyPatch
+        (
+            const word& name,
+            const dictionary& dict,
+            const label index,
+            const polyBoundaryMesh&
+        );
+
+        //- Construct as copy, resetting the boundary mesh
+        passiveProcessorPolyPatch
+        (
+            const passiveProcessorPolyPatch&,
+            const polyBoundaryMesh&
+        );
+
+        //- Construct as given the original patch and resetting the
+        //  face list and boundary mesh information
+        passiveProcessorPolyPatch
+        (
+            const passiveProcessorPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const label newSize,
+            const label newStart
+        );
+
+        //- Construct and return a clone, resetting the boundary mesh
+        virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
+        {
+            return autoPtr<polyPatch>(new passiveProcessorPolyPatch(*this, bm));
+        }
+
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const label newSize,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new passiveProcessorPolyPatch
+                (
+                    refCast<const passiveProcessorPolyPatch>(*this),
+                    bm,
+                    index,
+                    newSize,
+                    newStart
+                )
+            );
+        }
+
+
+    //- Destructor
+    virtual ~passiveProcessorPolyPatch();
+
+
+    // Member functions
+
+        //- Return processor number
+        int myProcNo() const
+        {
+            return myProcNo_;
+        }
+
+        //- Return neigbour processor number
+        int neighbProcNo() const
+        {
+            return neighbProcNo_;
+        }
+
+        //- Write the polyPatch data as a dictionary
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C

@@ -25,6 +25,7 @@ License
 
 #include "processorMeshesReconstructor.H"
 #include "processorPolyPatch.H"
+#include "passiveProcessorPolyPatch.H"
 #include "sharedPoints.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -45,7 +46,7 @@ Foam::label Foam::processorMeshesReconstructor::firstValidMesh() const
 
     return 0;
 }
-    
+
 
 bool Foam::processorMeshesReconstructor::readMapping()
 {
@@ -200,7 +201,7 @@ bool Foam::processorMeshesReconstructor::readMapping()
             }
         }
 
-        Info<< "Addressing from files: " << nl;
+        Info<< "Addressing from files: " << endl;
         forAll (meshes_, procI)
         {
             if (meshes_.set(procI))
@@ -359,6 +360,18 @@ void Foam::processorMeshesReconstructor::reconstructPoints(fvMesh& mesh) const
 Foam::autoPtr<Foam::fvMesh>
 Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
 {
+    // Note:
+    // In load balancing, there will exist a double set of processor patches
+    // One, created by moving cells adjacent to "old" processor boundaries
+    // (processorPolyPatch) and another, created by splitting up previously
+    // internal mesh faces into "new" processor patches
+    // In order to reconstruct the mesh correctly, the two sets of processor
+    // boundaries are kept separately
+    // In reconstruction, the order of processor patch faces needs to be
+    // preserved.  This is achieved by
+    // - first adding original processor patch faces
+    // - adding passiveProcessor patch faces in the processor order
+
     // Check for read
 
     if (readMapping())
@@ -407,34 +420,55 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
         }
     }
 
-    // Prepare patch reconstruction from processor 0
-    wordList reconPatchNames;
-    wordList reconPatchTypes;
+    // Prepare patch reconstruction
+    HashTable<label, word> patchNameLookup
+    (
+        meshes_[firstValidMesh()].boundaryMesh().size()
+    );
+    DynamicList<word> patchTypeLookup
+    (
+        meshes_[firstValidMesh()].boundaryMesh().size()
+    );
+
     label nReconPatches = 0;
 
-    // Grab names and types of patches, excluding processor patches
-    // Order must be identical on all processors
+    forAll (meshes_, procI)
     {
-        const polyBoundaryMesh& proc0Boundary =
-            meshes_[firstValidMesh()].boundaryMesh();
-
-        reconPatchNames.setSize(proc0Boundary.size());
-        reconPatchTypes.setSize(proc0Boundary.size());
-        nReconPatches = 0;
-
-        forAll (proc0Boundary, patchI)
+        if (meshes_.set(procI))
         {
-            if (!isA<processorPolyPatch>(proc0Boundary[patchI]))
+            const polyBoundaryMesh& procPatches = meshes_[procI].boundaryMesh();
+
+            forAll (procPatches, patchI)
             {
-                // Add patch name to list of patches
-                reconPatchNames[nReconPatches] = proc0Boundary[patchI].name();
-                reconPatchTypes[nReconPatches] = proc0Boundary[patchI].type();
-                nReconPatches++;
+                if (!isA<processorPolyPatch>(procPatches[patchI]))
+                {
+                    const word& patchName = procPatches[patchI].name();
+
+                    // Regular patch.  Try to find it in a list
+                    if (!patchNameLookup.found(patchName))
+                    {
+                        // Patch not found.  Add it
+                        patchNameLookup.insert(patchName, nReconPatches);
+                        patchTypeLookup.append(procPatches[patchI].type());
+                        nReconPatches++;
+                    }
+                }
             }
         }
+    }
 
-        reconPatchNames.setSize(nReconPatches);
-        reconPatchTypes.setSize(nReconPatches);
+    // Fill in patch names and types
+    wordList reconPatchNames(patchNameLookup.size());
+    wordList reconPatchTypes(patchNameLookup.size());
+
+    wordList patchNameToc = patchNameLookup.toc();
+
+    forAll (patchNameToc, pnI)
+    {
+        const label pnIndex = patchNameLookup.find(patchNameToc[pnI])();
+
+        reconPatchNames[pnIndex] = patchNameToc[pnI];
+        reconPatchTypes[pnIndex] = patchTypeLookup[pnIndex];
     }
 
     // Prepare point, face and patch reconstruction
@@ -442,6 +476,7 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
     label nReconFaces = 0;
     label nReconCells = 0;
     labelList reconPatchSizes(reconPatchTypes.size(), 0);
+    Pout<< "reconPatchNames: " << reconPatchNames << endl;
 
     forAll (meshes_, procI)
     {
@@ -454,40 +489,44 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
 
             const polyBoundaryMesh& procPatches = meshes_[procI].boundaryMesh();
 
-            nReconPatches = 0;
-
             forAll (procPatches, patchI)
             {
                 if (!isA<processorPolyPatch>(procPatches[patchI]))
                 {
+                    // Find processor patch index in reconstructed boundary
+                    const label pnIndex = patchNameLookup.find
+                    (
+                        procPatches[patchI].name()
+                    )();
+
                     // Check patch name and type
                     if
                     (
                         procPatches[patchI].name()
-                     != reconPatchNames[nReconPatches]
+                     != reconPatchNames[pnIndex]
                      || procPatches[patchI].type()
-                     != reconPatchTypes[nReconPatches]
+                     != reconPatchTypes[pnIndex]
                     )
                     {
                         FatalErrorIn
                         (
                             "autoPtr<fvMesh> processorMeshesReconstructor::"
                             "reconstructMesh(const Time& db)"
-                        )   << "Patch names, types or ordering does not match "
-                            << "across processors"
+                        )   << "Patch name and type does not match "
+                            << "across processors for patch "
+                            << procPatches[patchI].name() << " type: "
+                            << procPatches[patchI].type()
                             << abort(FatalError);
                     }
 
                     // Record number of faces in patch
-                    reconPatchSizes[nReconPatches] +=
-                        procPatches[patchI].size();
-                    nReconPatches++;
+                    reconPatchSizes[pnIndex] += procPatches[patchI].size();
                 }
             }
         }
     }
 
-    Info<< "Estimated max global mesh size (with duplicates): " << nl
+    Pout<< "Estimated max global mesh size (with duplicates): " << nl
         << "    nPoints = " << nReconPoints << nl
         << "    nFaces = " << nReconFaces << nl
         << "    nCells = " << nReconCells << nl
@@ -612,8 +651,9 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
 
     // Dump first valid mesh without checking
     {
+        Pout<< "Dump mesh 0" << endl;
         const label fvmId = firstValidMesh();
-        
+
         cellOffset[fvmId] = 0;
 
         const polyMesh& curMesh = meshes_[fvmId];
@@ -678,6 +718,9 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
                 const processorPolyPatch& procPatch =
                     refCast<const processorPolyPatch>(procPatches[patchI]);
 
+                // Record boundary-processor addressing: unmapped patch
+                bpAddr[patchI] = -1;
+                
                 for
                 (
                     label faceI = procPatch.start();
@@ -704,12 +747,19 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
             else
             {
                 // Regular patch: dump faces into patch face list
-                faceList& curRpFaces = reconPatchFaces[patchI];
-                labelList& curRpfOwner = reconPatchOwner[patchI];
-                label& nRpf = reconPatchSizes[patchI];
 
                 const polyPatch& curPatch = procPatches[patchI];
 
+                // Find processor patch index in reconstructed boundary
+                const label pnIndex = patchNameLookup.find(curPatch.name())();
+
+                // Record boundary-processor addressing: mapped patch
+                bpAddr[patchI] = pnIndex;
+                
+                faceList& curRpFaces = reconPatchFaces[pnIndex];
+                labelList& curRpfOwner = reconPatchOwner[pnIndex];
+                label& nRpf = reconPatchSizes[pnIndex];
+
                 for
                 (
                     label faceI = curPatch.start();
@@ -730,7 +780,6 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
                     // HJ, 16/Feb/2011
                     fpAddr[faceI] = nRpf + 1;
                     nRpf++;
-
                 }
             }
         }
@@ -740,15 +789,6 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
         {
             cpAddr[cellI] = cellI; // + cellOffset[firstValidMesh()];
         }
-
-        // Sort out boundary addressing: i for live patches, -1 for processor
-        bpAddr = -1;
-
-        // Note: loop over mapped patches
-        forAll (reconPatchSizes, patchI)
-        {
-            bpAddr[patchI] = patchI;
-        }
     }
 
 
@@ -756,6 +796,7 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
 
     for (label procI = 1; procI < meshes_.size(); procI++)
     {
+        Pout<< "Dump mesh " << procI << endl;
         // Grab cell offset from previous offset and mesh size
         cellOffset[procI] =
             cellOffset[procI - 1] + meshes_[procI - 1].nCells();
@@ -784,7 +825,7 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
 
             // Collect point-processor addressing for points on processor patches
 
-            // Go through all patches.  For neighbour patches, access
+            // Go through all processor patches.  For neighbour patches, access
             // owner addressing and dump into ppAddr
             forAll (procPatches, patchI)
             {
@@ -868,7 +909,7 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
                 }
             }
 
-            Info<< "Processor " << procI << " merged " << nMergedPoints
+            Pout<< "Processor " << procI << " merged " << nMergedPoints
                 << " points out of local " << curPoints.size()
                 << " and total " << nReconPoints << endl;
 
@@ -934,6 +975,9 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
                     const processorPolyPatch& procPatch =
                         refCast<const processorPolyPatch>(procPatches[patchI]);
 
+                    // Record boundary-processor addressing: unmapped patch
+                    bpAddr[patchI] = -1;
+                
                     // If patch is a master, drop the faces and fill the
                     // owner side addressing
                     if (procPatch.master())
@@ -1008,12 +1052,19 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
                 else
                 {
                     // Regular patch: dump faces into patch face list
-                    faceList& curRpFaces = reconPatchFaces[patchI];
-                    labelList& curRpfOwner = reconPatchOwner[patchI];
-                    label& nRpf = reconPatchSizes[patchI];
 
                     const polyPatch& curPatch = procPatches[patchI];
 
+                    // Find processor patch index in reconstructed boundary
+                    const label pnIndex = patchNameLookup.find(curPatch.name())();
+
+                    // Record boundary-processor addressing: mapped patch
+                    bpAddr[patchI] = pnIndex;
+                
+                    faceList& curRpFaces = reconPatchFaces[pnIndex];
+                    labelList& curRpfOwner = reconPatchOwner[pnIndex];
+                    label& nRpf = reconPatchSizes[pnIndex];
+
                     for
                     (
                         label faceI = curPatch.start();
@@ -1043,16 +1094,6 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
             {
                 cpAddr[cellI] = cellI + cellOffset[procI];
             }
-
-            // Sort out boundary addressing: i for live patches,
-            // -1 for processor
-            bpAddr = -1;
-
-            // Note: loop over mapped patches
-            forAll (reconPatchSizes, patchI)
-            {
-                bpAddr[patchI] = patchI;
-            }
         }
     }
 
@@ -1147,6 +1188,7 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
         }
     }
 
+    // Create global mesh with given region name
     autoPtr<fvMesh> globalMeshPtr
     (
         new fvMesh
@@ -1166,31 +1208,97 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
     );
     fvMesh& globalMesh = globalMeshPtr();
 
-    // Create patch list using mesh from processor 0
-    List<polyPatch*> reconPatches(nReconPatches);
+    // Create patch list by cloning meshes.  If all processors hold all live
+    // patches, it is sufficient to rebuilt the patches only from the first
+    // valid processor
+    // Note: 
+    List<polyPatch*> reconPatches(nReconPatches, NULL);
 
+    forAll (meshes_, procI)
     {
-        const polyBoundaryMesh& procPatches =
-            meshes_[firstValidMesh()].boundaryMesh();
-
-        reconPatches.setSize(reconPatchSizes.size());
-
-        forAll (reconPatchSizes, patchI)
+        if (meshes_.set(procI))
         {
-            reconPatches[patchI] =
-                procPatches[patchI].clone
-                (
-                    globalMesh.boundaryMesh(),
-                    patchI,
-                    reconPatchSizes[patchI],
-                    reconPatchStarts[patchI]
-                ).ptr();
+            const polyBoundaryMesh& procPatches = meshes_[procI].boundaryMesh();
+
+            forAll (procPatches, patchI)
+            {
+                // Processor patches have disappeared: skip them
+                if (!isA<processorPolyPatch>(procPatches[patchI]))
+                {
+                    // Find processor patch index in reconstructed boundary
+                    const label pnIndex = patchNameLookup.find
+                    (
+                        procPatches[patchI].name()
+                    )();
+
+                    // Check if the patch has already been set
+                    if (reconPatches[pnIndex] == NULL)
+                    {
+                        // Patch not set: clone it
+                        // Note: watch indices: setting pnIndex from patchI
+
+                        if (isA<passiveProcessorPolyPatch>(procPatches[patchI]))
+                        {
+                            // For a passive processor patch, create new
+                            // processor patch
+                            const passiveProcessorPolyPatch& ppPatch =
+                                refCast<const passiveProcessorPolyPatch>
+                                (
+                                    procPatches[patchI]
+                                );
+
+                            reconPatches[pnIndex] = new processorPolyPatch
+                            (
+                                ppPatch.name(),
+                                reconPatchSizes[pnIndex],
+                                reconPatchStarts[pnIndex],
+                                pnIndex,
+                                globalMesh.boundaryMesh(),
+                                Pstream::myProcNo(),   // Use correct local proc
+                                ppPatch.neighbProcNo()
+                            );
+                        }
+                        else
+                        {
+                            // Regular patch: clone
+                            reconPatches[pnIndex] =
+                                procPatches[patchI].clone
+                                (
+                                    globalMesh.boundaryMesh(),
+                                    patchI,
+                                    reconPatchSizes[pnIndex],
+                                    reconPatchStarts[pnIndex]
+                                ).ptr();
+                        }
+                    }
+                }
+            }
         }
     }
 
-    // Add both poly and fv boundary patches
-    globalMesh.addFvPatches(reconPatches);
-
+    // Check the list and fill in the missing slots
+    forAll (reconPatches, patchI)
+    {
+        if (reconPatches[patchI] == NULL)
+        {
+            // Patch not set.  Check its type
+            FatalErrorIn
+            (
+                "autoPtr<fvMesh> processorMeshesReconstructor::"
+                "reconstructMesh(const Time& db)"
+            )   << "Reconstructed patch " << patchI
+                << " name " << reconPatchNames[patchI]
+                << " type " << reconPatchTypes[patchI]
+                << " not set."
+                << abort(FatalError);
+        }
+    }
+    
+    // Add boundary patches to polyMesh and fvMesh
+    // Note: Mark boundary as invalid to disable analysis
+    // due to the presence of old/new patches
+    globalMesh.addFvPatches(reconPatches, false);
+    Pout<< "global boundary: " << globalMesh.boundaryMesh() << endl;
     // TODO: point, face and cell zones
 
     Info<< "Reconstructed addressing: " << nl;

src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.H

@@ -25,6 +25,14 @@ Class
     Foam::processorMeshesReconstructor
 
 Description
+    The processorMeshesReconstructor will take a list of meshes with
+    matching processor boundaries and build a single combined mesh by
+    matching processor patches to each other.
+
+    In order to 
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
 
 SourceFiles
     processorMeshesReconstructor.C

src/decompositionMethods/decompositionMethods/patchConstrainedDecomp/patchConstrainedDecomp.H

@@ -96,10 +96,9 @@ public:
         );
 
 
-    // Destructor
-
-        virtual ~patchConstrainedDecomp()
-        {}
+    //- Destructor
+    virtual ~patchConstrainedDecomp()
+    {}
 
 
     // Member Functions

src/dynamicMesh/loadBalanceFvMesh/loadBalanceFvMesh/loadBalanceFvMesh.C

@@ -136,12 +136,12 @@ bool Foam::loadBalanceFvMesh::update()
 
         return false;
     }
-            
+
     Info<< "Hello from loadBalanceFvMesh::update()" << endl;
 
     // Check imbalance.  Note: add run-time selection for the imbalance
     // weights criterion
-    
+
     // Create a parallel decomposition
     domainDecomposition meshDecomp
     (
@@ -152,7 +152,8 @@ bool Foam::loadBalanceFvMesh::update()
     // Decompose the mesh based on processor data
     meshDecomp.decomposeMesh(false);
 
-    // Analyse mesh decomposition to see how many cells are passed to which processor
+    // Analyse mesh decomposition to see how many cells are passed
+    // to which processor
     labelListList migratedCells(meshDecomp.nProcs());
 
     // Fill local list with number of local cells to be sent to each processor
@@ -171,10 +172,13 @@ bool Foam::loadBalanceFvMesh::update()
     // Now that each processor has filled in its own part, combine the data
     Pstream::gatherList(migratedCells);
     Pstream::scatterList(migratedCells);
-    Pout<< "migratedCells: " << migratedCells << endl;
+
     // Reading through second index now tells how many cells will arrive
     // from which processor
 
+    // Find out which processor faces will become local internal faces
+    // by comparing decomposition index from the other side
+
 
     // Split the mesh and fields over processors and send
     for (label procI = 0; procI < meshDecomp.nProcs(); procI++)
@@ -193,7 +197,8 @@ bool Foam::loadBalanceFvMesh::update()
                 (
                     procI,
                     time(),
-                    "processorPart" + Foam::name(procI)
+                    "processorPart" + Foam::name(procI),
+                    true    // Create passive processor patches
                 );
                 fvMesh& procMesh = procMeshPtr();
 
@@ -314,6 +319,76 @@ bool Foam::loadBalanceFvMesh::update()
         << reconMesh.boundary().size()
         << endl;
 
+    // Apply changes to the local mesh:
+    // - refactor the boundary to match new patches.  Note: processor
+    //   patch types may be added or removed
+    // - reset all primitives
+    Pout<< "Resetting mesh for load balancing.  Old boundary names: "
+        << boundaryMesh().names() << nl
+        << "New boundary names: "
+        << reconMesh.boundaryMesh().names() << nl
+        << endl;
+
+    polyBoundaryMesh& bMesh = const_cast<polyBoundaryMesh&>(boundaryMesh());
+
+    const polyBoundaryMesh& reconBMesh = reconMesh.boundaryMesh();
+
+    // Resize the existing boundary to match in the number of patches
+    if (bMesh.size() > reconBMesh.size())
+    {
+        // Decreasing in size: check patches that are due to disappear
+        Pout<< "New boundary is smaller: "
+            << bMesh.size() << " to " << reconBMesh.size() << endl;
+        // Check before resizing: only processor patches may be deleted
+        bool okToDelete = true;
+        
+        for (label patchI = reconBMesh.size(); patchI < bMesh.size(); patchI++)
+        {
+            if (!isA<processorPolyPatch>(bMesh[patchI]))
+            {
+                okToDelete = false;
+                break;
+            }
+        }
+
+        if (!okToDelete)
+        {
+            FatalErrorIn("bool loadBalanceFvMesh::update()")
+                << "Boundary resizing error: deleting non-processor patches "
+                << bMesh.size() << " to " << reconBMesh.size() << nl
+                << "old patch types: " << bMesh.types() << nl
+                << "new patch types: " << reconBMesh.types() << nl
+                << "This is not allowed."
+                << abort(FatalError);
+        }
+
+        // Resize the boundary
+        bMesh.setSize(reconBMesh.size());
+    }
+    else if (bMesh.size() < reconBMesh.size())
+    {
+        // Increasing in size: check patches that are due to disappear
+        Pout<< "New boundary is larger: "
+            << bMesh.size() << " to " << reconBMesh.size() << endl;
+
+        bMesh.setSize(reconBMesh.size());
+    }
+
+    // Delete processor patches from old boundary
+
+    
+    // Check alignment and types of old and new boundary
+    // forAll (bMesh, patchI)
+    // {
+    //     if (!bMesh[patchI].set())
+    //     {
+    //     }
+    // }
+
+    // To Do: reset mesh in polyMesh
+    // To Do: build a reconstructor from addressing data
+
+        
     // Create field reconstructor
     // fvFieldReconstructor fieldReconstructor
     // (
@@ -330,7 +405,7 @@ bool Foam::loadBalanceFvMesh::update()
     // Otherwise, they will be deleted.  HJ, 5/Mar/2018
     autoPtr<fvMesh> curMesh = procMeshes.set(Pstream::myProcNo(), NULL);
     curMesh.ptr();
-    
+
     return true;
 }