diff --git a/applications/utilities/parallelProcessing/decomposePar/distributeCells.C b/applications/utilities/parallelProcessing/decomposePar/distributeCells.C index 8477a7d087..ef7a613a11 100644 --- a/applications/utilities/parallelProcessing/decomposePar/distributeCells.C +++ b/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 { - - - // Work the faces whose neighbours need to be kept together into an - // agglomeration. - - // Per cell the region/agglomeration it is in - labelList cellToRegion(nCells(), -1); - - // Current region - label regionI = 0; - - labelHashSet freeRegions; + // Faces where owner and neighbour are not 'connected' (= all except + // sameProcFaces) + boolList blockedFace(nFaces(), true); forAllConstIter(labelHashSet, sameProcFaces, iter) { - label patchI = boundaryMesh().whichPatch(iter.key()); - - label own = faceOwner()[iter.key()]; - label nei = -1; - - if (patchI == -1) - { - nei = faceNeighbour()[iter.key()]; - } - else if (isA(boundaryMesh()[patchI])) - { - const cyclicPolyPatch& pp = - refCast(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); - } - } + blockedFace[iter.key()] = false; } + // Connect coupled boundary faces + const polyBoundaryMesh& patches = boundaryMesh(); - // Do all other cells - forAll(cellToRegion, cellI) + forAll(patches, patchI) { - if (cellToRegion[cellI] == -1) + const polyPatch& pp = patches[patchI]; + + if (pp.coupled()) { - cellToRegion[cellI] = regionI++; - } - } - - - // Compact out freeRegions - // ~~~~~~~~~~~~~~~~~~~~~~~ - - { - labelList compactRegion(regionI, -1); - - regionI = 0; - - forAll(compactRegion, i) - { - if (!freeRegions.found(compactRegion[i])) + forAll(pp, 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 " diff --git a/src/decompositionAgglomeration/decompositionMethods/decompositionMethod/decompositionMethod.C b/src/decompositionAgglomeration/decompositionMethods/decompositionMethod/decompositionMethod.C index 5d8453a1a2..b2a507ce42 100644 --- a/src/decompositionAgglomeration/decompositionMethods/decompositionMethod/decompositionMethod.C +++ b/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 > 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(); + } +} + + // ************************************************************************* // diff --git a/src/decompositionAgglomeration/decompositionMethods/decompositionMethod/decompositionMethod.H b/src/decompositionAgglomeration/decompositionMethods/decompositionMethod/decompositionMethod.H index 7f2804fe83..9d2095da44 100644 --- a/src/decompositionAgglomeration/decompositionMethods/decompositionMethod/decompositionMethod.H +++ b/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 New ( const dictionary& decompositionDict ); - //- Return a reference to the selected turbulence model + //- Return a reference to the selected decomposition method static autoPtr 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 pointField& + const labelList& cellToRegion, + 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; + }; diff --git a/src/decompositionAgglomeration/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H b/src/decompositionAgglomeration/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H index 7b1a3953cc..f7a3798f5f 100644 --- a/src/decompositionAgglomeration/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H +++ b/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); + } }; diff --git a/src/decompositionAgglomeration/decompositionMethods/manualDecomp/manualDecomp.H b/src/decompositionAgglomeration/decompositionMethods/manualDecomp/manualDecomp.H index 28bafc1007..d10888762a 100644 --- a/src/decompositionAgglomeration/decompositionMethods/manualDecomp/manualDecomp.H +++ b/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); + } }; diff --git a/src/decompositionAgglomeration/decompositionMethods/metisDecomp/metisDecomp.C b/src/decompositionAgglomeration/decompositionMethods/metisDecomp/metisDecomp.C index 884542e25b..ce44affbaa 100644 --- a/src/decompositionAgglomeration/decompositionMethods/metisDecomp/metisDecomp.C +++ b/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 polyMesh& mesh + const List& adjncy, + const List& xadj, + + List& 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 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(pbm[patchi])) - { - nInternalFaces += pbm[patchi].size(); - } - } - - // Create the adjncy array the size of the total number of internal and - // coupled faces - List 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 - (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(pbm[patchi])) - { - const unallocLabelList& faceCells = pbm[patchi].faceCells(); - - label sizeby2 = faceCells.size()/2; - - for (label facei=0; facei weights - ( - IOobject - ( - faceWeightsFile, - mesh_.time().timeName(), - mesh_, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ) - ); - - if (weights.size() != mesh_.nInternalFaces()) - { - FatalErrorIn("metisDecomp::decompose(const pointField&)") - << "Number of face weights " << weights.size() - << " does not equal number of internal faces " - << mesh_.nInternalFaces() - << exit(FatalError); - } - - // Assume symmetric weights. Keep same ordering as adjncy. - faceWeights.setSize(2*mesh_.nInternalFaces()); - - labelList nFacesPerCell(mesh_.nCells(), 0); - - for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++) - { - label w = weights[faceI]; - - label own = mesh_.faceOwner()[faceI]; - label nei = mesh_.faceNeighbour()[faceI]; - - faceWeights[xadj[own] + nFacesPerCell[own]++] = w; - faceWeights[xadj[nei] + nFacesPerCell[nei]++] = w; - } - } + //- faceWeights disabled. Only makes sense for cellCells from mesh. + //if (metisDecompCoeffs.found("faceWeightsFile")) + //{ + // Info<< "metisDecomp : Using face-based weights." << endl; + // + // word faceWeightsFile + // ( + // metisDecompCoeffs.lookup("faceWeightsFile") + // ); + // + // IOList weights + // ( + // IOobject + // ( + // faceWeightsFile, + // mesh_.time().timeName(), + // mesh_, + // IOobject::MUST_READ, + // IOobject::AUTO_WRITE + // ) + // ); + // + // if (weights.size() != adjncy.size()/2) + // { + // FatalErrorIn("metisDecomp::decompose(const pointField&)") + // << "Number of face weights " << weights.size() + // << " does not equal number of internal faces " + // << adjncy.size()/2 + // << exit(FatalError); + // } + // + // // Assume symmetric weights. Keep same ordering as adjncy. + // faceWeights.setSize(adjncy.size()); + // + // labelList nFacesPerCell(mesh_.nCells(), 0); + // + // for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++) + // { + // label w = weights[faceI]; + // + // label own = mesh_.faceOwner()[faceI]; + // label nei = mesh_.faceNeighbour()[faceI]; + // + // 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 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 - adjncy.begin(), // neighbour info + const_cast&>(xadj).begin(), // indexing into adjncy + const_cast&>(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 - adjncy.begin(), // neighbour info + const_cast&>(xadj).begin(), // indexing into adjncy + const_cast&>(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 - adjncy.begin(), // neighbour info + const_cast&>(xadj).begin(), // indexing into adjncy + const_cast&>(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 - adjncy.begin(), // neighbour info + const_cast&>(xadj).begin(), // indexing into adjncy + const_cast&>(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 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(pbm[patchi])) + { + nInternalFaces += pbm[patchi].size(); + } + } + + // Create the adjncy array the size of the total number of internal and + // coupled faces + List 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(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(pbm[patchi])) + { + const unallocLabelList& faceCells = pbm[patchi].faceCells(); + + label sizeby2 = faceCells.size()/2; + + for (label facei=0; facei 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 pointField& agglomPoints + const labelListList& cellCells, + List& adjncy, + List& xadj ) { - // Make Metis CSR (Compressed Storage Format) storage - // adjncy : contains neighbours (= edges in graph) - // xadj(celli) : start of information in adjncy for celli - - List xadj(agglomPoints.size()+1); - - // Get cellCells on coarse mesh. - labelListList cellCells(agglomPoints.size()); - { - List > 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]; - - forAll(cCells, i) - { - if (cCells[i] > coarseI) - { - nInternalFaces++; - } - } + nConnections += cellCells[coarseI].size(); } // Create the adjncy array as twice the size of the total number of // internal faces - List 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 - int numFlag = 0; - - // Method of decomposition - // recursive: multi-level recursive bisection (default) - // k-way: multi-level k-way - word method("k-way"); - - // decomposition options. 0 = use defaults - List options(5, 0); - - // processor weights initialised with no size, only used if specified in - // a file - Field processorWeights; - - // cell weights (so on the vertices of the dual) - List cellWeights; - - // Check for user supplied weights and decomp options - if (decompositionDict_.found("metisCoeffs")) +Foam::labelList Foam::metisDecomp::decompose +( + const labelList& agglom, + const pointField& agglomPoints +) +{ + if (agglom.size() != mesh_.nCells()) { - 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 adjncy; + List xadj; + { + // Get cellCells on coarse mesh. + labelListList cellCells; + calcCellCells + ( + mesh_, + agglom, + agglomPoints.size(), + cellCells ); - if (metisDecompCoeffs.found("method")) - { - 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 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); - } - } + calcMetisCSR(cellCells, adjncy, xadj); } - int numCells = cellCells.size(); - int nProcs = nProcessors_; - - // output: cell -> processor addressing - List 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() - ); - } - } + // Decompose using default weights + List finalDecomp; + decompose(adjncy, xadj, finalDecomp); // 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 adjncy; + List xadj; + calcMetisCSR(globalCellCells, adjncy, xadj); + + + // Decompose using default weights + List finalDecomp; + decompose(adjncy, xadj, finalDecomp); + + // Copy back to labelList + labelList decomp(finalDecomp.size()); + forAll(decomp, i) + { + decomp[i] = finalDecomp[i]; + } + return decomp; +} + + // ************************************************************************* // diff --git a/src/decompositionAgglomeration/decompositionMethods/metisDecomp/metisDecomp.H b/src/decompositionAgglomeration/decompositionMethods/metisDecomp/metisDecomp.H index b4aaef7a92..0dea688dc1 100644 --- a/src/decompositionAgglomeration/decompositionMethods/metisDecomp/metisDecomp.H +++ b/src/decompositionAgglomeration/decompositionMethods/metisDecomp/metisDecomp.H @@ -55,6 +55,13 @@ class metisDecomp // Private Member Functions + label decompose + ( + const List& adjncy, + const List& xadj, + List& 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& adjncy, + List& xadj + ); }; diff --git a/src/decompositionAgglomeration/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H b/src/decompositionAgglomeration/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H index 6e26f50257..f807edaf32 100644 --- a/src/decompositionAgglomeration/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H +++ b/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); + } }; diff --git a/src/decompositionAgglomeration/parMetisDecomp/parMetisDecomp.C b/src/decompositionAgglomeration/parMetisDecomp/parMetisDecomp.C index 155f4a7b6d..3178356cd1 100644 --- a/src/decompositionAgglomeration/parMetisDecomp/parMetisDecomp.C +++ b/src/decompositionAgglomeration/parMetisDecomp/parMetisDecomp.C @@ -32,6 +32,7 @@ License #include "polyMesh.H" #include "Time.H" #include "labelIOField.H" +#include "globalIndex.H" #include @@ -56,6 +57,267 @@ namespace Foam } +//- Does prevention of 0 cell domains and calls parmetis. +Foam::label Foam::parMetisDecomp::decompose +( + Field& xadj, + Field& adjncy, + const pointField& cellCentres, + Field& cellWeights, + Field& faceWeights, + const List& options, + + List& finalDecomp +) +{ + // C style numbering + int numFlag = 0; + + // Number of dimensions + int nDims = 3; + + // Get number of cells on all processors + List nLocalCells(Pstream::nProcs()); + nLocalCells[Pstream::myProcNo()] = xadj.size()-1; + Pstream::gatherList(nLocalCells); + Pstream::scatterList(nLocalCells); + + // Get cell offsets. + List 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 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 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 prevXadj(fromPrevProc); + Field prevAdjncy(fromPrevProc); + Field prevXyz(fromPrevProc); + Field prevCellWeights(fromPrevProc); + Field 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::subField(xadj, nCells, startCell)-startFace + << Field::subField(adjncy, nFaces, startFace) + << SubField(xyz, nDims*nCells, nDims*startCell) + << + ( + (cellWeights.size() > 0) + ? static_cast&> + ( + Field::subField(cellWeights, nCells, startCell) + ) + : Field(0) + ) + << + ( + (faceWeights.size() > 0) + ? static_cast&> + ( + Field::subField(faceWeights, nFaces, startFace) + ) + : Field(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 tpwgts(nCon*nProcessors_, 1./nProcessors_); + // Imbalance tolerance + Field 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&>(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 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 + ( + 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 nLocalCells(Pstream::nProcs()); nLocalCells[Pstream::myProcNo()] = mesh_.nCells(); Pstream::gatherList(nLocalCells); Pstream::scatterList(nLocalCells); // Get cell offsets. - labelList cellOffsets(Pstream::nProcs()+1); - label nGlobalCells = 0; + List cellOffsets(Pstream::nProcs()+1); + 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 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 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 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 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 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 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 cellWeights; // face weights (so on the edges of the dual) - labelField faceWeights; + Field 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) } - - // 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) - 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 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(xyz, nDims*nCells, nDims*startCell) - << - ( - (cellWeights.size() > 0) - ? static_cast - ( - labelField::subField(cellWeights, nCells, startCell) - ) - : labelField(0) - ) - << - ( - (faceWeights.size() > 0) - ? static_cast - ( - 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 tpwgts(nCon*nProcessors_, 1./nProcessors_); - // Imbalance tolerance - Field 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 + // Do actual decomposition + List finalDecomp; + decompose ( - cellOffsets.begin(), // vtxDist - xadj.begin(), - adjncy.begin(), - vwgtPtr, // vertexweights - adjwgtPtr, // edgeweights - &wgtFlag, - &numFlag, - &nDims, - xyz.begin(), - &nCon, - &nProcessors_, // nParts - tpwgts.begin(), - ubvec.begin(), - options.begin(), - &edgeCut, - finalDecomp.begin(), - &comm + xadj, + adjncy, + points, + cellWeights, + faceWeights, + options, + + finalDecomp ); - - // If we sent cells across make sure we undo it - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - // Receive back from next processor if I sent something - if (nSendCells[Pstream::myProcNo()] > 0) + // Copy back to labelList + labelList decomp(finalDecomp.size()); + forAll(decomp, i) { - 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 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 << - SubList