new decomposition options; fixed decomposePar preserveFaceZones

applications/utilities/parallelProcessing/decomposePar/distributeCells.C

@@ -28,6 +28,8 @@ License
 #include "decompositionMethod.H"
 #include "cpuTime.H"
 #include "cyclicPolyPatch.H"
+#include "cellSet.H"
+#include "regionSplit.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -127,117 +129,34 @@ void domainDecomposition::distributeCells()
     }
     else
     {
-        
+        // Faces where owner and neighbour are not 'connected' (= all except
-
+        // sameProcFaces)
-        // Work the faces whose neighbours need to be kept together into an
+        boolList blockedFace(nFaces(), true);
-        // agglomeration.
-
-        // Per cell the region/agglomeration it is in
-        labelList cellToRegion(nCells(), -1);
-
-        // Current region
-        label regionI = 0;
-
-        labelHashSet freeRegions;
 
         forAllConstIter(labelHashSet, sameProcFaces, iter)
         {
-            label patchI = boundaryMesh().whichPatch(iter.key());
+            blockedFace[iter.key()] = false;
-
-            label own = faceOwner()[iter.key()];
-            label nei = -1;
-
-            if (patchI == -1)
-            {
-                nei = faceNeighbour()[iter.key()];
-            }
-            else if (isA<cyclicPolyPatch>(boundaryMesh()[patchI]))
-            {
-                const cyclicPolyPatch& pp =
-                    refCast<const cyclicPolyPatch>(boundaryMesh()[patchI]);
-
-                nei = faceOwner()[pp.transformGlobalFace(iter.key())];
-            }
-
-            if (nei != -1)
-            {
-                label ownRegion = cellToRegion[own];
-                label neiRegion = cellToRegion[nei];
-
-                if (ownRegion == -1 && neiRegion == -1)
-                {
-                    // Allocate new agglomeration
-                    cellToRegion[own] = regionI;
-                    cellToRegion[nei] = regionI;
-                    regionI++;
-                }
-                else if (ownRegion != -1)
-                {
-                    // Owner already part of agglomeration. Add nei to it.
-                    cellToRegion[nei] = ownRegion;
-                }
-                else if (neiRegion != -1)
-                {
-                    // nei already part of agglomeration. Add own to it.
-                    cellToRegion[own] = neiRegion;
-                }
-                else if (ownRegion < neiRegion)
-                {
-                    // Renumber neiRegion
-                    forAll(cellToRegion, cellI)
-                    {
-                        if (cellToRegion[cellI] == neiRegion)
-                        {
-                            cellToRegion[cellI] = ownRegion;
-                        }
-                    }
-                    freeRegions.insert(neiRegion);
-                }
-                else if (ownRegion > neiRegion)
-                {
-                    // Renumber ownRegion
-                    forAll(cellToRegion, cellI)
-                    {
-                        if (cellToRegion[cellI] == ownRegion)
-                        {
-                            cellToRegion[cellI] = neiRegion;
-                        }
-                    }
-                    freeRegions.insert(ownRegion);
-                }
-            }
         }
 
+        // Connect coupled boundary faces
+        const polyBoundaryMesh& patches =  boundaryMesh();
 
-        // Do all other cells
+        forAll(patches, patchI)
-        forAll(cellToRegion, cellI)
         {
-            if (cellToRegion[cellI] == -1)
+            const polyPatch& pp = patches[patchI];
+
+            if (pp.coupled())
             {
-                cellToRegion[cellI] = regionI++;
+                forAll(pp, i)
-            }
-        }
-
-
-        // Compact out freeRegions
-        // ~~~~~~~~~~~~~~~~~~~~~~~
-
-        {
-            labelList compactRegion(regionI, -1);
-
-            regionI = 0;
-
-            forAll(compactRegion, i)
-            {
-                if (!freeRegions.found(compactRegion[i]))
                 {
-                    compactRegion[i] = regionI++;
+                    blockedFace[pp.start()+i] = false;
                 }
             }
-
-            inplaceRenumber(compactRegion, cellToRegion);
         }
 
+        // Determine global regions, separated by blockedFaces
+        regionSplit globalRegion(*this, blockedFace);
+
 
         // Determine region cell centres
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -249,11 +168,11 @@ void domainDecomposition::distributeCells()
 
         const point greatPoint(GREAT, GREAT, GREAT);
 
-        pointField regionCentres(regionI, greatPoint);
+        pointField regionCentres(globalRegion.nRegions(), greatPoint);
 
-        forAll(cellToRegion, cellI)
+        forAll(globalRegion, cellI)
         {
-            label regionI = cellToRegion[cellI];
+            label regionI = globalRegion[cellI];
 
             if (regionCentres[regionI] == greatPoint)
             {
@@ -261,10 +180,9 @@ void domainDecomposition::distributeCells()
             }
         }
 
-
         // Do decomposition on agglomeration
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        cellToProc_ = decomposePtr().decompose(cellToRegion, regionCentres);
+        cellToProc_ = decomposePtr().decompose(globalRegion, regionCentres);
     }
 
     Info<< "\nFinished decomposition in "

src/decompositionAgglomeration/decompositionMethods/decompositionMethod/decompositionMethod.C

@@ -125,4 +125,50 @@ Foam::labelList Foam::decompositionMethod::decompose
 }
 
 
+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].shrink());
+        dynCellCells[coarseI].clear();
+    }
+}
+
+
 // ************************************************************************* //

src/decompositionAgglomeration/decompositionMethods/decompositionMethod/decompositionMethod.H

@@ -36,9 +36,6 @@ SourceFiles
 #ifndef decompositionMethod_H
 #define decompositionMethod_H
 
-#include "typeInfo.H"
-#include "runTimeSelectionTables.H"
-#include "dictionary.H"
 #include "polyMesh.H"
 #include "pointField.H"
 
@@ -60,6 +57,15 @@ protected:
         label nProcessors_;
 
 
+        //- Helper: determine (non-parallel) cellCells from mesh agglomeration.
+        static void calcCellCells
+        (
+            const polyMesh& mesh,
+            const labelList& agglom,
+            const label nCoarse,
+            labelListList& cellCells
+        );
+
 private:
 
     // Private Member Functions
@@ -103,13 +109,13 @@ public:
 
     // Selectors
 
-        //- Return a reference to the selected turbulence model
+        //- Return a reference to the selected decomposition method
         static autoPtr<decompositionMethod> New
         (
             const dictionary& decompositionDict
         );
 
-        //- Return a reference to the selected turbulence model
+        //- Return a reference to the selected decomposition method
         static autoPtr<decompositionMethod> New
         (
             const dictionary& decompositionDict,
@@ -142,18 +148,35 @@ public:
         //  proc boundaries)
         virtual bool parallelAware() const = 0;
 
-        //- Return for every coordinate the wanted processor number
+        //- Return for every coordinate the wanted processor number. Use the
+        //  mesh connectivity (if needed)
         virtual labelList decompose(const pointField&) = 0;
 
         //- 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.
+        //  functionality natively. Coarse cells are local to the processor
+        //  (if in parallel). If you want to have coarse cells spanning
+        //  processors use the next function below instead.
         virtual labelList decompose
         (
-            const labelList& agglom,
+            const labelList& cellToRegion,
-            const pointField&
+            const pointField& regionPoints
         );
+
+        //- 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
+        virtual labelList decompose
+        (
+            const labelListList& globalCellCells,
+            const pointField& cc
+        ) = 0;
+
 };
 
 

src/decompositionAgglomeration/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H

@@ -154,7 +154,25 @@ public:
             return true;
         }
 
+        //- Return for every coordinate the wanted processor number. Use the
+        //  mesh connectivity (if needed)
         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
+        virtual labelList decompose
+        (
+            const labelListList& globalCellCells,
+            const pointField& cc
+        )
+        {
+            return decompose(cc);
+        }
 };
 
 

src/decompositionAgglomeration/decompositionMethods/manualDecomp/manualDecomp.H

@@ -97,7 +97,25 @@ public:
             return true;
         }
 
+        //- Return for every coordinate the wanted processor number. Use the
+        //  mesh connectivity (if needed)
         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
+        virtual labelList decompose
+        (
+            const labelListList& globalCellCells,
+            const pointField& cc
+        )
+        {
+            return decompose(cc);
+        }
 };
 
 

src/decompositionAgglomeration/decompositionMethods/metisDecomp/metisDecomp.C

@@ -29,7 +29,7 @@ License
 #include "floatScalar.H"
 #include "IFstream.H"
 #include "Time.H"
-#include "coupledPolyPatch.H"
+#include "cyclicPolyPatch.H"
 
 extern "C"
 {
@@ -52,113 +52,18 @@ namespace Foam
     );
 }
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::metisDecomp::metisDecomp
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Call Metis with options from dictionary.
+Foam::label Foam::metisDecomp::decompose
 (
-    const dictionary& decompositionDict,
+    const List<int>& adjncy,
-    const polyMesh& mesh
+    const List<int>& xadj,
+
+    List<int>& finalDecomp
 )
-:
-    decompositionMethod(decompositionDict),
-    mesh_(mesh)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::labelList Foam::metisDecomp::decompose(const pointField& points)
 {
-    // Make Metis CSR (Compressed Storage Format) storage
-    //   adjncy      : contains neighbours (= edges in graph)
-    //   xadj(celli) : start of information in adjncy for celli
-
-    List<int> xadj(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<coupledPolyPatch>(pbm[patchi]))
-        {
-            nInternalFaces += pbm[patchi].size();
-        }
-    }
-
-    // Create the adjncy array the size of the total number of internal and
-    // coupled faces
-    List<int> adjncy(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<coupledPolyPatch>
-                (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
-    forAll(pbm, patchi)
-    {
-        if (isA<coupledPolyPatch>(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;
-            }
-        }
-    }
-
-
     // C style numbering
     int numFlag = 0;
 
@@ -259,68 +164,69 @@ Foam::labelList Foam::metisDecomp::decompose(const pointField& points)
             );
             cellWeights.transfer(cellIOWeights);
 
-            if (cellWeights.size() != mesh_.nCells())
+            if (cellWeights.size() != xadj.size()-1)
             {
                 FatalErrorIn("metisDecomp::decompose(const pointField&)")
                     << "Number of cell weights " << cellWeights.size()
-                    << " does not equal number of cells " << mesh_.nCells()
+                    << " does not equal number of cells " << xadj.size()-1
                     << exit(FatalError);
             }
         }
 
-        if (metisDecompCoeffs.found("faceWeightsFile"))
+        //- faceWeights disabled. Only makes sense for cellCells from mesh.
-        {
+        //if (metisDecompCoeffs.found("faceWeightsFile"))
-            Info<< "metisDecomp : Using face-based weights." << endl;
+        //{
-
+        //    Info<< "metisDecomp : Using face-based weights." << endl;
-            word faceWeightsFile
+        //
-            (
+        //    word faceWeightsFile
-                metisDecompCoeffs.lookup("faceWeightsFile")
+        //    (
-            );
+        //        metisDecompCoeffs.lookup("faceWeightsFile")
-
+        //    );
-            IOList<int> weights
+        //
-            (
+        //    IOList<int> weights
-                IOobject
+        //    (
-                (
+        //        IOobject
-                    faceWeightsFile,
+        //        (
-                    mesh_.time().timeName(),
+        //            faceWeightsFile,
-                    mesh_,
+        //            mesh_.time().timeName(),
-                    IOobject::MUST_READ,
+        //            mesh_,
-                    IOobject::AUTO_WRITE
+        //            IOobject::MUST_READ,
-                )
+        //            IOobject::AUTO_WRITE
-            );
+        //        )
-
+        //    );
-            if (weights.size() != mesh_.nInternalFaces())
+        //
-            {
+        //    if (weights.size() != adjncy.size()/2)
-                FatalErrorIn("metisDecomp::decompose(const pointField&)")
+        //    {
-                    << "Number of face weights " << weights.size()
+        //        FatalErrorIn("metisDecomp::decompose(const pointField&)")
-                    << " does not equal number of internal faces "
+        //            << "Number of face weights " << weights.size()
-                    << mesh_.nInternalFaces()
+        //            << " does not equal number of internal faces "
-                    << exit(FatalError);
+        //            << adjncy.size()/2
-            }
+        //            << exit(FatalError);
-
+        //    }
-            // Assume symmetric weights. Keep same ordering as adjncy.
+        //
-            faceWeights.setSize(2*mesh_.nInternalFaces());
+        //    // Assume symmetric weights. Keep same ordering as adjncy.
-
+        //    faceWeights.setSize(adjncy.size());
-            labelList nFacesPerCell(mesh_.nCells(), 0);
+        //
-
+        //    labelList nFacesPerCell(mesh_.nCells(), 0);
-            for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+        //
-            {
+        //    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
-                label w = weights[faceI];
+        //    {
-
+        //        label w = weights[faceI];
-                label own = mesh_.faceOwner()[faceI];
+        //
-                label nei = mesh_.faceNeighbour()[faceI];
+        //        label own = mesh_.faceOwner()[faceI];
-
+        //        label nei = mesh_.faceNeighbour()[faceI];
-                faceWeights[xadj[own] + nFacesPerCell[own]++] = w;
+        //
-                faceWeights[xadj[nei] + nFacesPerCell[nei]++] = w;
+        //        faceWeights[xadj[own] + nFacesPerCell[own]++] = w;
-            }
+        //        faceWeights[xadj[nei] + nFacesPerCell[nei]++] = w;
-        }
+        //    }
+        //}
     }
 
-    int numCells = mesh_.nCells();
+    int numCells = xadj.size()-1;
     int nProcs = nProcessors_;
 
     // output: cell -> processor addressing
-    List<int> finalDecomp(mesh_.nCells());
+    finalDecomp.setSize(numCells);
 
     // output: number of cut edges
     int edgeCut = 0;
@@ -348,8 +254,8 @@ Foam::labelList Foam::metisDecomp::decompose(const pointField& points)
             METIS_WPartGraphRecursive
             (
                 &numCells,         // num vertices in graph
-                xadj.begin(),      // indexing into adjncy
+                const_cast<List<int>&>(xadj).begin(),   // indexing into adjncy
-                adjncy.begin(),    // neighbour info
+                const_cast<List<int>&>(adjncy).begin(), // neighbour info
                 vwgtPtr,           // vertexweights
                 adjwgtPtr,         // no edgeweights
                 &wgtFlag,
@@ -366,8 +272,8 @@ Foam::labelList Foam::metisDecomp::decompose(const pointField& points)
             METIS_PartGraphRecursive
             (
                 &numCells,         // num vertices in graph
-                xadj.begin(),      // indexing into adjncy
+                const_cast<List<int>&>(xadj).begin(),   // indexing into adjncy
-                adjncy.begin(),    // neighbour info
+                const_cast<List<int>&>(adjncy).begin(), // neighbour info
                 vwgtPtr,           // vertexweights
                 adjwgtPtr,         // no edgeweights
                 &wgtFlag,
@@ -386,8 +292,8 @@ Foam::labelList Foam::metisDecomp::decompose(const pointField& points)
             METIS_WPartGraphKway
             (
                 &numCells,         // num vertices in graph
-                xadj.begin(),      // indexing into adjncy
+                const_cast<List<int>&>(xadj).begin(),   // indexing into adjncy
-                adjncy.begin(),    // neighbour info
+                const_cast<List<int>&>(adjncy).begin(), // neighbour info
                 vwgtPtr,           // vertexweights
                 adjwgtPtr,         // no edgeweights
                 &wgtFlag,
@@ -404,8 +310,8 @@ Foam::labelList Foam::metisDecomp::decompose(const pointField& points)
             METIS_PartGraphKway
             (
                 &numCells,         // num vertices in graph
-                xadj.begin(),      // indexing into adjncy
+                const_cast<List<int>&>(xadj).begin(),   // indexing into adjncy
-                adjncy.begin(),    // neighbour info
+                const_cast<List<int>&>(adjncy).begin(), // neighbour info
                 vwgtPtr,           // vertexweights
                 adjwgtPtr,         // no edgeweights
                 &wgtFlag,
@@ -418,6 +324,131 @@ Foam::labelList Foam::metisDecomp::decompose(const pointField& points)
         }
     }
 
+    return edgeCut;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::metisDecomp::metisDecomp
+(
+    const dictionary& decompositionDict,
+    const polyMesh& mesh
+)
+:
+    decompositionMethod(decompositionDict),
+    mesh_(mesh)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::labelList Foam::metisDecomp::decompose(const pointField& points)
+{
+    if (points.size() != mesh_.nCells())
+    {
+        FatalErrorIn("metisDecomp::decompose(const pointField&)")
+            << "Can use this decomposition method only for the whole mesh"
+            << endl
+            << "and supply one coordinate (cellCentre) for every cell." << endl
+            << "The number of coordinates " << points.size() << endl
+            << "The number of cells in the mesh " << mesh_.nCells()
+            << exit(FatalError);
+    }
+
+    // Make Metis CSR (Compressed Storage Format) storage
+    //   adjncy      : contains neighbours (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
+
+    List<int> xadj(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
+    List<int> adjncy(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;
+            }
+        }
+    }
+
+    // Decompose using default weights
+    List<int> finalDecomp;
+    decompose(adjncy, xadj, finalDecomp);
+
+    // Copy back to labelList
     labelList decomp(finalDecomp.size());
     forAll(decomp, i)
     {
@@ -427,68 +458,28 @@ Foam::labelList Foam::metisDecomp::decompose(const pointField& points)
 }
 
 
-Foam::labelList Foam::metisDecomp::decompose
+// From cell-cell connections to Metis format (like CompactListList)
+void Foam::metisDecomp::calcMetisCSR
 (
-    const labelList& agglom,
+    const labelListList& cellCells,
-    const pointField& agglomPoints
+    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
-
-    List<int> xadj(agglomPoints.size()+1);
-
-    // Get cellCells on coarse mesh.
-    labelListList cellCells(agglomPoints.size());
-    {
-        List<DynamicList<label> > dynCellCells(cellCells.size());
-
-        forAll(mesh_.faceNeighbour(), faceI)
-        {
-            label own = agglom[mesh_.faceOwner()[faceI]];
-            label nei = agglom[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);
-                }
-            }
-        }
-
-        forAll(dynCellCells, coarseI)
-        {
-            cellCells[coarseI].transfer(dynCellCells[coarseI].shrink());
-            dynCellCells[coarseI].clear();
-        }
-    }
-
-
     // Count number of internal faces
-    label nInternalFaces = 0;
+    label nConnections = 0;
 
     forAll(cellCells, coarseI)
     {
-        const labelList& cCells = cellCells[coarseI];
+        nConnections += cellCells[coarseI].size();
-
-        forAll(cCells, i)
-        {
-            if (cCells[i] > coarseI)
-            {
-                nInternalFaces++;
-            }
-        }
     }
 
     // Create the adjncy array as twice the size of the total number of
     // internal faces
-    List<int> adjncy(2*nInternalFaces);
+    adjncy.setSize(nConnections);
+
+    xadj.setSize(cellCells.size()+1);
+
 
     // Fill in xadj
     // ~~~~~~~~~~~~
@@ -506,211 +497,47 @@ Foam::labelList Foam::metisDecomp::decompose
         }
     }
     xadj[cellCells.size()] = freeAdj;
+}
 
 
-    // C style numbering
+Foam::labelList Foam::metisDecomp::decompose
-    int numFlag = 0;
+(
-
+    const labelList& agglom,
-    // Method of decomposition
+    const pointField& agglomPoints
-    // recursive: multi-level recursive bisection (default)
+)
-    // k-way: multi-level k-way 
+{
-    word method("k-way");
+    if (agglom.size() != mesh_.nCells())
-
-    // decomposition options. 0 = use defaults
-    List<int> options(5, 0);
-
-    // processor weights initialised with no size, only used if specified in
-    // a file
-    Field<floatScalar> processorWeights;
-
-    // cell weights (so on the vertices of the dual)
-    List<int> cellWeights;
-
-    // Check for user supplied weights and decomp options
-    if (decompositionDict_.found("metisCoeffs"))
     {
-        dictionary metisDecompCoeffs
+        FatalErrorIn
         (
-            decompositionDict_.subDict("metisCoeffs")
+            "parMetisDecomp::decompose(const labelList&, const pointField&)"
+        )   << "Size of cell-to-coarse map " << agglom.size()
+            << " differs from number of cells in mesh " << mesh_.nCells()
+            << exit(FatalError);
+    }
+
+    // Make Metis CSR (Compressed Storage Format) storage
+    //   adjncy      : contains neighbours (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
+    List<int> adjncy;
+    List<int> xadj;
+    {
+        // Get cellCells on coarse mesh.
+        labelListList cellCells;
+        calcCellCells
+        (
+            mesh_,
+            agglom,
+            agglomPoints.size(),
+            cellCells
         );
 
-        if (metisDecompCoeffs.found("method"))
+        calcMetisCSR(cellCells, adjncy, xadj);
-        {
-            metisDecompCoeffs.lookup("method") >> method;
-
-            if (method != "recursive" && method != "k-way")
-            {
-                FatalErrorIn("metisDecomp::decompose()")
-                    << "Method " << method << " in metisCoeffs in dictionary : "
-                    << decompositionDict_.name()
-                    << " should be 'recursive' or 'k-way'"
-                    << exit(FatalError);
-            }
-
-            Info<< "metisDecomp : Using Metis options     " << options
-                << endl << endl;
-        }
-
-        if (metisDecompCoeffs.found("options"))
-        {
-            metisDecompCoeffs.lookup("options") >> options;
-
-            if (options.size() != 5)
-            {
-                FatalErrorIn("metisDecomp::decompose()")
-                    << "Number of options in metisCoeffs in dictionary : "
-                    << decompositionDict_.name()
-                    << " should be 5"
-                    << exit(FatalError);
-            }
-
-            Info<< "metisDecomp : Using Metis options     " << options
-                << endl << endl;
-        }
-
-        if (metisDecompCoeffs.found("processorWeights"))
-        {
-            metisDecompCoeffs.lookup("processorWeights") >> processorWeights;
-            processorWeights /= sum(processorWeights);
-
-            if (processorWeights.size() != nProcessors_)
-            {
-                FatalErrorIn("metisDecomp::decompose(const pointField&)")
-                    << "Number of processor weights "
-                    << processorWeights.size()
-                    << " does not equal number of domains " << nProcessors_
-                    << exit(FatalError);
-            }
-        }
-
-        if (metisDecompCoeffs.found("cellWeightsFile"))
-        {
-            Info<< "metisDecomp : Using cell-based weights." << endl;
-
-            word cellWeightsFile
-            (
-                metisDecompCoeffs.lookup("cellWeightsFile")
-            );
-
-            IOList<int> cellIOWeights
-            (
-                IOobject
-                (
-                    cellWeightsFile,
-                    mesh_.time().timeName(),
-                    mesh_,
-                    IOobject::MUST_READ,
-                    IOobject::AUTO_WRITE
-                )
-            );
-            cellWeights.transfer(cellIOWeights);
-
-            if (cellWeights.size() != cellCells.size())
-            {
-                FatalErrorIn("metisDecomp::decompose(const pointField&)")
-                    << "Number of cell weights " << cellWeights.size()
-                    << " does not equal number of agglomerated cells "
-                    << cellCells.size() << exit(FatalError);
-            }
-        }
     }
 
-    int numCells = cellCells.size();
+    // Decompose using default weights
-    int nProcs = nProcessors_;
+    List<int> finalDecomp;
-
+    decompose(adjncy, xadj, finalDecomp);
-    // output: cell -> processor addressing
-    List<int> finalDecomp(cellCells.size());
-
-    // output: number of cut edges
-    int edgeCut = 0;
-
-    // Vertex weight info
-    int wgtFlag = 0;
-    int* vwgtPtr = NULL;
-    int* adjwgtPtr = NULL;
-
-    if (cellWeights.size() > 0)
-    {
-        vwgtPtr = cellWeights.begin();
-        wgtFlag += 2;       // Weights on vertices
-    }
-
-    if (method == "recursive")
-    {
-        if (processorWeights.size())
-        {
-            METIS_WPartGraphRecursive
-            (
-                &numCells,         // num vertices in graph
-                xadj.begin(),      // indexing into adjncy
-                adjncy.begin(),    // neighbour info
-                vwgtPtr,           // vertexweights
-                adjwgtPtr,         // no edgeweights
-                &wgtFlag,
-                &numFlag,
-                &nProcs,
-                processorWeights.begin(),
-                options.begin(),
-                &edgeCut,
-                finalDecomp.begin()
-            );
-        }
-        else
-        {
-            METIS_PartGraphRecursive
-            (
-                &numCells,         // num vertices in graph
-                xadj.begin(),      // indexing into adjncy
-                adjncy.begin(),    // neighbour info
-                vwgtPtr,           // vertexweights
-                adjwgtPtr,         // no edgeweights
-                &wgtFlag,
-                &numFlag,
-                &nProcs,
-                options.begin(),
-                &edgeCut,
-                finalDecomp.begin()
-            );
-        }
-    }
-    else
-    {
-        if (processorWeights.size())
-        {
-            METIS_WPartGraphKway
-            (
-                &numCells,         // num vertices in graph
-                xadj.begin(),      // indexing into adjncy
-                adjncy.begin(),    // neighbour info
-                vwgtPtr,           // vertexweights
-                adjwgtPtr,         // no edgeweights
-                &wgtFlag,
-                &numFlag,
-                &nProcs,
-                processorWeights.begin(),
-                options.begin(),
-                &edgeCut,
-                finalDecomp.begin()
-            );
-        }
-        else
-        {
-            METIS_PartGraphKway
-            (
-                &numCells,         // num vertices in graph
-                xadj.begin(),      // indexing into adjncy
-                adjncy.begin(),    // neighbour info
-                vwgtPtr,           // vertexweights
-                adjwgtPtr,         // no edgeweights
-                &wgtFlag,
-                &numFlag,
-                &nProcs,
-                options.begin(),
-                &edgeCut,
-                finalDecomp.begin()
-            );
-        }
-    }
 
 
     // Rework back into decomposition for original mesh_
@@ -725,4 +552,44 @@ Foam::labelList Foam::metisDecomp::decompose
 }
 
 
+Foam::labelList Foam::metisDecomp::decompose
+(
+    const labelListList& globalCellCells,
+    const pointField& cellCentres
+)
+{
+    if (cellCentres.size() != globalCellCells.size())
+    {
+        FatalErrorIn
+        (
+            "metisDecomp::decompose(const pointField&, const labelListList&)"
+        )   << "Inconsistent number of cells (" << globalCellCells.size()
+            << ") and number of cell centres (" << cellCentres.size()
+            << ")." << exit(FatalError);
+    }
+
+
+    // Make Metis CSR (Compressed Storage Format) storage
+    //   adjncy      : contains neighbours (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
+
+    List<int> adjncy;
+    List<int> xadj;
+    calcMetisCSR(globalCellCells, adjncy, xadj);
+
+
+    // Decompose using default weights
+    List<int> finalDecomp;
+    decompose(adjncy, xadj, finalDecomp);
+
+    // Copy back to labelList
+    labelList decomp(finalDecomp.size());
+    forAll(decomp, i)
+    {
+        decomp[i] = finalDecomp[i];
+    }
+    return decomp;
+}
+
+
 // ************************************************************************* //

src/decompositionAgglomeration/decompositionMethods/metisDecomp/metisDecomp.H

@@ -55,6 +55,13 @@ class metisDecomp
 
     // Private Member Functions
 
+        label decompose
+        (
+            const List<int>& adjncy,
+            const List<int>& xadj,
+            List<int>& finalDecomp
+        );
+
         //- Disallow default bitwise copy construct and assignment
         void operator=(const metisDecomp&);
         metisDecomp(const metisDecomp&);
@@ -90,9 +97,40 @@ public:
             return false;
         }
 
+        //- Return for every coordinate the wanted processor number. Use the
+        //  mesh connectivity (if needed)
         virtual labelList decompose(const pointField&);
 
-        virtual labelList decompose(const labelList& agglom, 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.
+        virtual labelList decompose
+        (
+            const labelList& agglom,
+            const pointField&
+        );
+
+        //- Return for every coordinate the wanted processor number. Explicitly
+        //  provided mesh connectivity.
+        //  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
+        virtual labelList decompose
+        (
+            const labelListList& globalCellCells,
+            const pointField& cc
+        );
+
+        //- Helper to convert cellcells into Metis storage
+        static void calcMetisCSR
+        (
+            const labelListList& globalCellCells,
+            List<int>& adjncy,
+            List<int>& xadj
+        );
 
 };
 

src/decompositionAgglomeration/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H

@@ -92,6 +92,16 @@ public:
         }
 
         virtual labelList decompose(const pointField&);
+
+        //- Explicitly provided connectivity
+        virtual labelList decompose
+        (
+            const labelListList& globalCellCells,
+            const pointField& cc
+        )
+        {
+            return decompose(cc);
+        }
 };
 
 

src/decompositionAgglomeration/parMetisDecomp/parMetisDecomp.C

@@ -32,6 +32,7 @@ License
 #include "polyMesh.H"
 #include "Time.H"
 #include "labelIOField.H"
+#include "globalIndex.H"
 
 #include <mpi.h>
 
@@ -56,6 +57,267 @@ namespace Foam
 }
 
 
+//- Does prevention of 0 cell domains and calls parmetis.
+Foam::label Foam::parMetisDecomp::decompose
+(
+    Field<int>& xadj,
+    Field<int>& adjncy,
+    const pointField& cellCentres,
+    Field<int>& cellWeights,
+    Field<int>& faceWeights,
+    const List<int>& options,
+
+    List<int>& finalDecomp
+)
+{
+    // C style numbering
+    int numFlag = 0;
+
+    // Number of dimensions
+    int nDims = 3;
+
+    // Get number of cells on all processors
+    List<int> nLocalCells(Pstream::nProcs());
+    nLocalCells[Pstream::myProcNo()] = xadj.size()-1;
+    Pstream::gatherList(nLocalCells);
+    Pstream::scatterList(nLocalCells);
+
+    // Get cell offsets.
+    List<int> cellOffsets(Pstream::nProcs()+1);
+    int nGlobalCells = 0;
+    forAll(nLocalCells, procI)
+    {
+        cellOffsets[procI] = nGlobalCells;
+        nGlobalCells += nLocalCells[procI];
+    }
+    cellOffsets[Pstream::nProcs()] = nGlobalCells;
+
+    // Convert pointField into float
+    Field<floatScalar> xyz(3*cellCentres.size());
+    int compI = 0;
+    forAll(cellCentres, cellI)
+    {
+        const point& cc = cellCentres[cellI];
+        xyz[compI++] = float(cc.x());
+        xyz[compI++] = float(cc.y());
+        xyz[compI++] = float(cc.z());
+    }
+
+    // Make sure every domain has at least one cell
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // (Metis falls over with zero sized domains)
+    // Trickle cells from processors that have them down to those that
+    // don't.
+
+
+    // Number of cells to send down (is same as number of cells next processor
+    // has to receive)
+    List<int> nSendCells(Pstream::nProcs(), 0);
+
+    for (label procI = nLocalCells.size()-1; procI >=1; procI--)
+    {
+        if (nLocalCells[procI]-nSendCells[procI] < 1)
+        {
+            nSendCells[procI-1] = nSendCells[procI]-nLocalCells[procI]+1;
+        }
+    }
+
+    // First receive (so increasing the sizes of all arrays)
+
+    if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0)
+    {
+        // Receive cells from previous processor
+        IPstream fromPrevProc(Pstream::blocking, Pstream::myProcNo()-1);
+
+        Field<int> prevXadj(fromPrevProc);
+        Field<int> prevAdjncy(fromPrevProc);
+        Field<floatScalar> prevXyz(fromPrevProc);
+        Field<int> prevCellWeights(fromPrevProc);
+        Field<int> prevFaceWeights(fromPrevProc);
+
+        // Insert adjncy
+        prepend(prevAdjncy, adjncy);
+        // Adapt offsets and prepend xadj
+        xadj += prevAdjncy.size();
+        prepend(prevXadj, xadj);
+        // Coords
+        prepend(prevXyz, xyz);
+        // Weights
+        prepend(prevCellWeights, cellWeights);
+        prepend(prevFaceWeights, faceWeights);
+    }
+
+
+    // Send to my next processor
+
+    if (nSendCells[Pstream::myProcNo()] > 0)
+    {
+        // Send cells to next processor
+        OPstream toNextProc(Pstream::blocking, Pstream::myProcNo()+1);
+
+        int nCells = nSendCells[Pstream::myProcNo()];
+        int startCell = xadj.size()-1 - nCells;
+        int startFace = xadj[startCell];
+        int nFaces = adjncy.size()-startFace;
+
+        // Send for all cell data: last nCells elements
+        // Send for all face data: last nFaces elements
+        toNextProc
+            << Field<int>::subField(xadj, nCells, startCell)-startFace
+            << Field<int>::subField(adjncy, nFaces, startFace)
+            << SubField<floatScalar>(xyz, nDims*nCells, nDims*startCell)
+            <<
+            (
+                (cellWeights.size() > 0)
+              ? static_cast<const Field<int>&>
+                (
+                    Field<int>::subField(cellWeights, nCells, startCell)
+                )
+              : Field<int>(0)
+            )
+            <<
+            (
+                (faceWeights.size() > 0)
+              ? static_cast<const Field<int>&>
+                (
+                    Field<int>::subField(faceWeights, nFaces, startFace)
+                )
+              : Field<int>(0)
+            );
+
+        // Remove data that has been sent
+        if (faceWeights.size() > 0)
+        {
+            faceWeights.setSize(faceWeights.size()-nFaces);
+        }
+        if (cellWeights.size() > 0)
+        {
+            cellWeights.setSize(cellWeights.size()-nCells);
+        }
+        xyz.setSize(xyz.size()-nDims*nCells);
+        adjncy.setSize(adjncy.size()-nFaces);
+        xadj.setSize(xadj.size() - nCells);
+    }
+
+
+
+    // Adapt number of cells
+    forAll(nSendCells, procI)
+    {
+        // Sent cells
+        nLocalCells[procI] -= nSendCells[procI];
+
+        if (procI >= 1)
+        {
+            // Received cells
+            nLocalCells[procI] += nSendCells[procI-1];
+        }
+    }
+    // Adapt cellOffsets
+    nGlobalCells = 0;
+    forAll(nLocalCells, procI)
+    {
+        cellOffsets[procI] = nGlobalCells;
+        nGlobalCells += nLocalCells[procI];
+    }
+
+
+    // Weight info
+    int wgtFlag = 0;
+    int* vwgtPtr = NULL;
+    int* adjwgtPtr = NULL;
+
+    if (cellWeights.size() > 0)
+    {
+        vwgtPtr = cellWeights.begin();
+        wgtFlag += 2;       // Weights on vertices
+    }
+    if (faceWeights.size() > 0)
+    {
+        adjwgtPtr = faceWeights.begin();
+        wgtFlag += 1;       // Weights on edges
+    }
+
+
+    // Number of weights or balance constraints
+    int nCon = 1;
+    // Per processor, per constraint the weight
+    Field<floatScalar> tpwgts(nCon*nProcessors_, 1./nProcessors_);
+    // Imbalance tolerance
+    Field<floatScalar> ubvec(nCon, 1.02);
+    if (nProcessors_ == 1)
+    {
+        // If only one processor there is no imbalance.
+        ubvec[0] = 1;
+    }
+
+    MPI_Comm comm = MPI_COMM_WORLD;
+
+    // output: cell -> processor addressing
+    finalDecomp.setSize(nLocalCells[Pstream::myProcNo()]);
+
+    // output: number of cut edges
+    int edgeCut = 0;
+
+
+    ParMETIS_V3_PartGeomKway
+    (
+        cellOffsets.begin(),    // vtxDist
+        xadj.begin(),
+        adjncy.begin(),
+        vwgtPtr,                // vertexweights
+        adjwgtPtr,              // edgeweights
+        &wgtFlag,
+        &numFlag,
+        &nDims,
+        xyz.begin(),
+        &nCon,
+        &nProcessors_,          // nParts
+        tpwgts.begin(),
+        ubvec.begin(),
+        const_cast<List<int>&>(options).begin(),
+        &edgeCut,
+        finalDecomp.begin(),
+        &comm
+    );
+
+
+    // If we sent cells across make sure we undo it
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Receive back from next processor if I sent something
+    if (nSendCells[Pstream::myProcNo()] > 0)
+    {
+        IPstream fromNextProc(Pstream::blocking, Pstream::myProcNo()+1);
+
+        List<int> nextFinalDecomp(fromNextProc);
+
+        append(nextFinalDecomp, finalDecomp);
+    }
+
+    // Send back to previous processor.
+    if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0)
+    {
+        OPstream toPrevProc(Pstream::blocking, Pstream::myProcNo()-1);
+
+        int nToPrevious = nSendCells[Pstream::myProcNo()-1];
+
+        toPrevProc <<
+            SubList<int>
+            (
+                finalDecomp,
+                nToPrevious,
+                finalDecomp.size()-nToPrevious
+            );
+
+        // Remove locally what has been sent
+        finalDecomp.setSize(finalDecomp.size()-nToPrevious);
+    }
+
+    return edgeCut;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::parMetisDecomp::parMetisDecomp
@@ -73,30 +335,35 @@ Foam::parMetisDecomp::parMetisDecomp
 
 Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
 {
+    if (points.size() != mesh_.nCells())
+    {
+        FatalErrorIn("parMetisDecomp::decompose(const pointField&)")
+            << "Can use this decomposition method only for the whole mesh"
+            << endl
+            << "and supply one coordinate (cellCentre) for every cell." << endl
+            << "The number of coordinates " << points.size() << endl
+            << "The number of cells in the mesh " << mesh_.nCells()
+            << exit(FatalError);
+    }
+
     // For running sequential ...
     if (Pstream::nProcs() <= 1)
     {
         return metisDecomp(decompositionDict_, mesh_).decompose(points);
     }
 
-    //
-    // Make Metis Distributed CSR (Compressed Storage Format) storage
-    //   adjncy      : contains cellCells (= edges in graph)
-    //   xadj(celli) : start of information in adjncy for celli
-    //
-
     // Create global cell numbers
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Get number of cells on all processors
-    labelList nLocalCells(Pstream::nProcs());
+    List<int> nLocalCells(Pstream::nProcs());
     nLocalCells[Pstream::myProcNo()] = mesh_.nCells();
     Pstream::gatherList(nLocalCells);
     Pstream::scatterList(nLocalCells);
 
     // Get cell offsets.
-    labelList cellOffsets(Pstream::nProcs()+1);
+    List<int> cellOffsets(Pstream::nProcs()+1);
-    label nGlobalCells = 0;
+    int nGlobalCells = 0;
     forAll(nLocalCells, procI)
     {
         cellOffsets[procI] = nGlobalCells;
@@ -104,7 +371,14 @@ Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
     }
     cellOffsets[Pstream::nProcs()] = nGlobalCells;
 
-    label myOffset = cellOffsets[Pstream::myProcNo()];
+    int myOffset = cellOffsets[Pstream::myProcNo()];
+
+
+    //
+    // Make Metis Distributed CSR (Compressed Storage Format) storage
+    //   adjncy      : contains cellCells (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
+    //
 
 
 
@@ -116,7 +390,7 @@ Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
     // Get renumbered owner on other side of coupled faces
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    labelList globalNeighbour(mesh_.nFaces()-mesh_.nInternalFaces());
+    List<int> globalNeighbour(mesh_.nFaces()-mesh_.nInternalFaces());
 
     forAll(patches, patchI)
     {
@@ -142,7 +416,7 @@ Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Number of faces per cell
-    labelList nFacesPerCell(mesh_.nCells(), 0);
+    List<int> nFacesPerCell(mesh_.nCells(), 0);
 
     // Number of coupled faces
     label nCoupledFaces = 0;
@@ -167,15 +441,15 @@ Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
                 nFacesPerCell[faceOwner[faceI++]]++;
             }
         }
-    }            
+    }
 
 
     // Fill in xadj
     // ~~~~~~~~~~~~
 
-    labelField xadj(mesh_.nCells()+1, -1);
+    Field<int> xadj(mesh_.nCells()+1, -1);
 
-    label freeAdj = 0;
+    int freeAdj = 0;
 
     for (label cellI = 0; cellI < mesh_.nCells(); cellI++)
     {
@@ -190,7 +464,7 @@ Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
     // Fill in adjncy
     // ~~~~~~~~~~~~~~
 
-    labelField adjncy(2*mesh_.nInternalFaces() + nCoupledFaces, -1);
+    Field<int> adjncy(2*mesh_.nInternalFaces() + nCoupledFaces, -1);
 
     nFacesPerCell = 0;
 
@@ -223,41 +497,21 @@ Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
                 bFaceI++;
             }
         }
-    }            
-
-
-
-    // C style numbering
-    int numFlag = 0;
-
-    // Number of dimensions
-    int nDims = 3;
-
-    // cell centres
-    Field<floatScalar> xyz(nDims*mesh_.nCells());
-    const pointField& cellCentres = mesh_.cellCentres();
-    label compI = 0;
-    forAll(cellCentres, cellI)
-    {
-        const point& cc = cellCentres[cellI];
-        xyz[compI++] = float(cc.x());
-        xyz[compI++] = float(cc.y());
-        xyz[compI++] = float(cc.z());
     }
 
 
+
     // decomposition options. 0 = use defaults
-    labelList options(3, 0);
+    List<int> options(3, 0);
     //options[0] = 1;     // don't use defaults but use values below
     //options[1] = -1;    // full debug info
     //options[2] = 15;    // random number seed
 
-
     // cell weights (so on the vertices of the dual)
-    labelField cellWeights;
+    Field<int> cellWeights;
 
     // face weights (so on the edges of the dual)
-    labelField faceWeights;
+    Field<int> faceWeights;
 
     // Check for user supplied weights and decomp options
     if (decompositionDict_.found("metisCoeffs"))
@@ -364,7 +618,7 @@ Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
                         faceI++;
                     }
                 }
-            }            
+            }
         }
 
         if (parMetisDecompCoeffs.found("options"))
@@ -384,219 +638,353 @@ Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
     }
 
 
-
+    // Do actual decomposition
-    // Make sure every domain has at least one cell
+    List<int> finalDecomp;
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    decompose
-    // (Metis falls over with zero sized domains)
-    // Trickle cells from processors that have them down to those that
-    // don't.
-
-
-    // Number of cells to send down (is same as number of cells next processor
-    // has to receive)
-    labelList nSendCells(Pstream::nProcs(), 0);
-
-    for (label procI = nLocalCells.size()-1; procI >=1; procI--)
-    {
-        if (nLocalCells[procI]-nSendCells[procI] < 1)
-        {
-            nSendCells[procI-1] = nSendCells[procI]-nLocalCells[procI]+1;
-        }
-    }
-
-    // First receive (so increasing the sizes of all arrays)
-
-    if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0)
-    {
-        // Receive cells from previous processor
-        IPstream fromPrevProc(Pstream::blocking, Pstream::myProcNo()-1);
-
-        labelField prevXadj(fromPrevProc);
-        labelField prevAdjncy(fromPrevProc);
-        Field<floatScalar> prevXyz(fromPrevProc);
-        labelField prevCellWeights(fromPrevProc);
-        labelField prevFaceWeights(fromPrevProc);
-
-        // Insert adjncy
-        prepend(prevAdjncy, adjncy);
-        // Adapt offsets and prepend xadj
-        xadj += prevAdjncy.size();
-        prepend(prevXadj, xadj);
-        // Coords
-        prepend(prevXyz, xyz);
-        // Weights
-        prepend(prevCellWeights, cellWeights);
-        prepend(prevFaceWeights, faceWeights);
-    }
-
-
-    // Send to my next processor
-
-    if (nSendCells[Pstream::myProcNo()] > 0)
-    {
-        // Send cells to next processor
-        OPstream toNextProc(Pstream::blocking, Pstream::myProcNo()+1);
-
-        label nCells = nSendCells[Pstream::myProcNo()];
-        label startCell = xadj.size()-1 - nCells;
-        label startFace = xadj[startCell];
-        label nFaces = adjncy.size()-startFace;
-
-        // Send for all cell data: last nCells elements
-        // Send for all face data: last nFaces elements
-        toNextProc
-            << labelField::subField(xadj, nCells, startCell)-startFace
-            << labelField::subField(adjncy, nFaces, startFace)
-            << SubField<floatScalar>(xyz, nDims*nCells, nDims*startCell)
-            <<
-            (
-                (cellWeights.size() > 0)
-              ? static_cast<const labelField&>
-                (
-                    labelField::subField(cellWeights, nCells, startCell)
-                )
-              : labelField(0)
-            )
-            <<
-            (
-                (faceWeights.size() > 0)
-              ? static_cast<const labelField&>
-                (
-                    labelField::subField(faceWeights, nFaces, startFace)
-                )
-              : labelField(0)
-            );
-
-        // Remove data that has been sent
-        if (faceWeights.size() > 0)
-        {
-            faceWeights.setSize(faceWeights.size()-nFaces);
-        }
-        if (cellWeights.size() > 0)
-        {
-            cellWeights.setSize(cellWeights.size()-nCells);
-        }
-        xyz.setSize(xyz.size()-nDims*nCells);
-        adjncy.setSize(adjncy.size()-nFaces);
-        xadj.setSize(xadj.size() - nCells);
-    }
-
-
-
-    // Adapt number of cells
-    forAll(nSendCells, procI)
-    {
-        // Sent cells
-        nLocalCells[procI] -= nSendCells[procI];
-
-        if (procI >= 1)
-        {
-            // Received cells
-            nLocalCells[procI] += nSendCells[procI-1];
-        }
-    }
-    // Adapt cellOffsets
-    nGlobalCells = 0;
-    forAll(nLocalCells, procI)
-    {
-        cellOffsets[procI] = nGlobalCells;
-        nGlobalCells += nLocalCells[procI];
-    }
-
-
-    // Weight info
-    int wgtFlag = 0;
-    label* vwgtPtr = NULL;
-    label* adjwgtPtr = NULL;
-
-    if (cellWeights.size() > 0)
-    {
-        vwgtPtr = cellWeights.begin();
-        wgtFlag += 2;       // Weights on vertices
-    }
-    if (faceWeights.size() > 0)
-    {
-        adjwgtPtr = faceWeights.begin();
-        wgtFlag += 1;       // Weights on edges
-    }
-
-
-    // Number of weights or balance constraints
-    int nCon = 1;
-    // Per processor, per constraint the weight
-    Field<floatScalar> tpwgts(nCon*nProcessors_, 1./nProcessors_);
-    // Imbalance tolerance
-    Field<floatScalar> ubvec(nCon, 1.02);
-    if (nProcessors_ == 1)
-    {
-        // If only one processor there is no imbalance.
-        ubvec[0] = 1;
-    }
-
-    MPI_Comm comm = MPI_COMM_WORLD;
-
-    // output: cell -> processor addressing
-    labelList finalDecomp(nLocalCells[Pstream::myProcNo()]);
-
-    // output: number of cut edges
-    int edgeCut = 0;
-
-
-    ParMETIS_V3_PartGeomKway
     (
-        cellOffsets.begin(),    // vtxDist
+        xadj,
-        xadj.begin(),
+        adjncy,
-        adjncy.begin(),
+        points,
-        vwgtPtr,                // vertexweights
+        cellWeights,
-        adjwgtPtr,              // edgeweights
+        faceWeights,
-        &wgtFlag,
+        options,
-        &numFlag,
+
-        &nDims,
+        finalDecomp
-        xyz.begin(),
-        &nCon,
-        &nProcessors_,          // nParts
-        tpwgts.begin(),
-        ubvec.begin(),
-        options.begin(),
-        &edgeCut,
-        finalDecomp.begin(),
-        &comm
     );
 
-
+    // Copy back to labelList
-    // If we sent cells across make sure we undo it
+    labelList decomp(finalDecomp.size());
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    forAll(decomp, i)
-
-    // Receive back from next processor if I sent something
-    if (nSendCells[Pstream::myProcNo()] > 0)
     {
-        IPstream fromNextProc(Pstream::blocking, Pstream::myProcNo()+1);
+        decomp[i] = finalDecomp[i];
+    }
+    return decomp;
+}
 
-        labelList nextFinalDecomp(fromNextProc);
 
-        append(nextFinalDecomp, finalDecomp);
+Foam::labelList Foam::parMetisDecomp::decompose
+(
+    const labelList& cellToRegion,
+    const pointField& regionPoints
+)
+{
+    const labelList& faceOwner = mesh_.faceOwner();
+    const labelList& faceNeighbour = mesh_.faceNeighbour();
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    if (cellToRegion.size() != mesh_.nCells())
+    {
+        FatalErrorIn
+        (
+            "parMetisDecomp::decompose(const labelList&, const pointField&)"
+        )   << "Size of cell-to-coarse map " << cellToRegion.size()
+            << " differs from number of cells in mesh " << mesh_.nCells()
+            << exit(FatalError);
     }
 
-    // Send back to previous processor.
+
-    if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0)
+    // Global region numbering engine
+    globalIndex globalRegions(regionPoints.size());
+
+
+    // Get renumbered owner region on other side of coupled faces
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    List<int> globalNeighbour(mesh_.nFaces()-mesh_.nInternalFaces());
+
+    forAll(patches, patchI)
     {
-        OPstream toPrevProc(Pstream::blocking, Pstream::myProcNo()-1);
+        const polyPatch& pp = patches[patchI];
 
-        label nToPrevious = nSendCells[Pstream::myProcNo()-1];
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+            label bFaceI = pp.start() - mesh_.nInternalFaces();
 
-        toPrevProc <<
+            forAll(pp, i)
-            SubList<label>
+            {
+                label ownRegion = cellToRegion[faceOwner[faceI]];
+                globalNeighbour[bFaceI++] = globalRegions.toGlobal(ownRegion);
+                faceI++;
+            }
+        }
+    }
+
+    // Get the cell on the other side of coupled patches
+    syncTools::swapBoundaryFaceList(mesh_, globalNeighbour, false);
+
+
+    // Get globalCellCells on coarse mesh
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelListList globalRegionRegions;
+    {
+        List<DynamicList<label> > dynRegionRegions(regionPoints.size());
+
+        // Internal faces first
+        forAll(faceNeighbour, faceI)
+        {
+            label ownRegion = cellToRegion[faceOwner[faceI]];
+            label neiRegion = cellToRegion[faceNeighbour[faceI]];
+
+            if (ownRegion != neiRegion)
+            {
+                label globalOwn = globalRegions.toGlobal(ownRegion);
+                label globalNei = globalRegions.toGlobal(neiRegion);
+
+                if (findIndex(dynRegionRegions[ownRegion], globalNei) == -1)
+                {
+                    dynRegionRegions[ownRegion].append(globalNei);
+                }
+                if (findIndex(dynRegionRegions[neiRegion], globalOwn) == -1)
+                {
+                    dynRegionRegions[neiRegion].append(globalOwn);
+                }
+            }
+        }
+
+        // Coupled boundary faces
+        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 ownRegion = cellToRegion[faceOwner[faceI]];
+                    label globalNei = globalNeighbour[bFaceI++];
+                    faceI++;
+
+                    if (findIndex(dynRegionRegions[ownRegion], globalNei) == -1)
+                    {
+                        dynRegionRegions[ownRegion].append(globalNei);
+                    }
+                }
+            }
+        }
+
+        globalRegionRegions.setSize(dynRegionRegions.size());
+        forAll(dynRegionRegions, i)
+        {
+            globalRegionRegions[i].transfer(dynRegionRegions[i].shrink());
+            dynRegionRegions[i].clear();
+        }
+    }
+
+    labelList regionDecomp(decompose(globalRegionRegions, regionPoints));
+
+    // Rework back into decomposition for original mesh_
+    labelList cellDistribution(cellToRegion.size());
+
+    forAll(cellDistribution, cellI)
+    {
+        cellDistribution[cellI] = regionDecomp[cellToRegion[cellI]];
+    }
+    return cellDistribution;
+}
+
+
+Foam::labelList Foam::parMetisDecomp::decompose
+(
+    const labelListList& globalCellCells,
+    const pointField& cellCentres
+)
+{
+    if (cellCentres.size() != globalCellCells.size())
+    {
+        FatalErrorIn
+        (
+            "parMetisDecomp::decompose(const labelListList&, const pointField&)"
+        )   << "Inconsistent number of cells (" << globalCellCells.size()
+            << ") and number of cell centres (" << cellCentres.size()
+            << ")." << exit(FatalError);
+    }
+
+    // For running sequential ...
+    if (Pstream::nProcs() <= 1)
+    {
+        return metisDecomp(decompositionDict_, mesh_)
+            .decompose(globalCellCells, cellCentres);
+    }
+
+
+    // Make Metis Distributed CSR (Compressed Storage Format) storage
+
+    // Connections
+    Field<int> adjncy;
+    // Offsets into adjncy
+    Field<int> xadj;
+    metisDecomp::calcMetisCSR(globalCellCells, adjncy, xadj);
+
+    // decomposition options. 0 = use defaults
+    List<int> options(3, 0);
+    //options[0] = 1;     // don't use defaults but use values below
+    //options[1] = -1;    // full debug info
+    //options[2] = 15;    // random number seed
+
+    // cell weights (so on the vertices of the dual)
+    Field<int> cellWeights;
+
+    // face weights (so on the edges of the dual)
+    Field<int> faceWeights;
+
+    // Check for user supplied weights and decomp options
+    if (decompositionDict_.found("metisCoeffs"))
+    {
+        dictionary parMetisDecompCoeffs
+        (
+            decompositionDict_.subDict("metisCoeffs")
+        );
+
+        if (parMetisDecompCoeffs.found("cellWeightsFile"))
+        {
+            word cellWeightsFile
             (
-                finalDecomp,
+                parMetisDecompCoeffs.lookup("cellWeightsFile")
-                nToPrevious,
-                finalDecomp.size()-nToPrevious
             );
 
-        // Remove locally what has been sent
+            Info<< "parMetisDecomp : Using cell-based weights read from "
-        finalDecomp.setSize(finalDecomp.size()-nToPrevious);
+                << cellWeightsFile << endl;
+
+            labelIOField cellIOWeights
+            (
+                IOobject
+                (
+                    cellWeightsFile,
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                )
+            );
+            cellWeights.transfer(cellIOWeights);
+
+            if (cellWeights.size() != cellCentres.size())
+            {
+                FatalErrorIn
+                (
+                    "parMetisDecomp::decompose"
+                    "(const labelListList&, const pointField&)"
+                )   << "Number of cell weights " << cellWeights.size()
+                    << " read from " << cellIOWeights.objectPath()
+                    << " does not equal number of cells " << cellCentres.size()
+                    << exit(FatalError);
+            }
+        }
+
+        //- faceWeights disabled. Only makes sense for cellCells from mesh.
+        //if (parMetisDecompCoeffs.found("faceWeightsFile"))
+        //{
+        //    word faceWeightsFile
+        //    (
+        //        parMetisDecompCoeffs.lookup("faceWeightsFile")
+        //    );
+        //
+        //    Info<< "parMetisDecomp : Using face-based weights read from "
+        //        << faceWeightsFile << endl;
+        //
+        //    labelIOField weights
+        //    (
+        //        IOobject
+        //        (
+        //            faceWeightsFile,
+        //            mesh_.time().timeName(),
+        //            mesh_,
+        //            IOobject::MUST_READ,
+        //            IOobject::AUTO_WRITE
+        //        )
+        //    );
+        //
+        //    if (weights.size() != mesh_.nFaces())
+        //    {
+        //        FatalErrorIn("parMetisDecomp::decompose(const pointField&)")
+        //            << "Number of face weights " << weights.size()
+        //            << " does not equal number of internal and boundary faces "
+        //            << mesh_.nFaces()
+        //            << exit(FatalError);
+        //    }
+        //
+        //    faceWeights.setSize(2*mesh_.nInternalFaces()+nCoupledFaces);
+        //
+        //    // Assume symmetric weights. Keep same ordering as adjncy.
+        //    nFacesPerCell = 0;
+        //
+        //    // Handle internal faces
+        //    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+        //    {
+        //        label w = weights[faceI];
+        //
+        //        label own = faceOwner[faceI];
+        //        label nei = faceNeighbour[faceI];
+        //
+        //        faceWeights[xadj[own] + nFacesPerCell[own]++] = w;
+        //        faceWeights[xadj[nei] + nFacesPerCell[nei]++] = w;
+        //    }
+        //    // Coupled boundary faces
+        //    forAll(patches, patchI)
+        //    {
+        //        const polyPatch& pp = patches[patchI];
+        //
+        //        if (pp.coupled())
+        //        {
+        //            label faceI = pp.start();
+        //
+        //            forAll(pp, i)
+        //            {
+        //                label w = weights[faceI];
+        //                label own = faceOwner[faceI];
+        //                adjncy[xadj[own] + nFacesPerCell[own]++] = w;
+        //                faceI++;
+        //            }
+        //        }
+        //    }
+        //}
+
+        if (parMetisDecompCoeffs.found("options"))
+        {
+            parMetisDecompCoeffs.lookup("options") >> options;
+
+            Info<< "Using Metis options     " << options
+                << endl << endl;
+
+            if (options.size() != 3)
+            {
+                FatalErrorIn
+                (
+                    "parMetisDecomp::decompose"
+                    "(const labelListList&, const pointField&)"
+                )   << "Number of options " << options.size()
+                    << " should be three." << exit(FatalError);
+            }
+        }
     }
 
-    return finalDecomp;
+
+    // Do actual decomposition
+    List<int> finalDecomp;
+    decompose
+    (
+        xadj,
+        adjncy,
+        cellCentres,
+        cellWeights,
+        faceWeights,
+        options,
+
+        finalDecomp
+    );
+
+    // Copy back to labelList
+    labelList decomp(finalDecomp.size());
+    forAll(decomp, i)
+    {
+        decomp[i] = finalDecomp[i];
+    }
+    return decomp;
 }
 
 

src/decompositionAgglomeration/parMetisDecomp/parMetisDecomp.H

@@ -62,6 +62,18 @@ class parMetisDecomp
         template<class Type>
         static void append(const UList<Type>&, List<Type>&);
 
+        label decompose
+        (
+            Field<int>& xadj,
+            Field<int>& adjncy,
+            const pointField& cellCentres,
+            Field<int>& cellWeights,
+            Field<int>& faceWeights,
+            const List<int>& options,
+
+            List<int>& finalDecomp
+        );
+
 
         //- Disallow default bitwise copy construct and assignment
         void operator=(const parMetisDecomp&);
@@ -93,12 +105,37 @@ public:
     // Member Functions
 
         //- parMetis handles Foam processor boundaries
-        bool parallelAware() const
+        virtual bool parallelAware() const
         {
             return true;
         }
 
-        labelList decompose(const pointField&);
+        //- Return for every coordinate the wanted processor number. Use the
+        //  mesh connectivity (if needed)
+        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.
+        virtual labelList decompose
+        (
+            const labelList& agglom,
+            const pointField&
+        );
+
+        //- Return for every coordinate the wanted processor number. Explicitly
+        //  provided mesh connectivity.
+        //  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
+        virtual labelList decompose
+        (
+            const labelListList& globalCellCells,
+            const pointField& cc
+        );
 };