src/decompositionMethods/parMetisDecomp/parMetisDecomp.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     |
   5     \\  /    A nd           | For copyright notice see file Copyright
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "parMetisDecomp.H"
  27 #include "metisDecomp.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "floatScalar.H"
  30 #include "Time.H"
  31 #include "labelIOField.H"
  32 #include "syncTools.H"
  33 #include "globalIndex.H"
  34
  35 #include <mpi.h>
  36
  37 extern "C"
  38 {
  39 #   include "parmetis.h"
  40 }
  41
  42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  43
  44 namespace Foam
  45 {
  46     defineTypeNameAndDebug(parMetisDecomp, 0);
  47
  48     addToRunTimeSelectionTable
  49     (
  50         decompositionMethod,
  51         parMetisDecomp,
  52         dictionaryMesh
  53     );
  54 }
  55
  56
  57 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  58
  59 //- Does prevention of 0 cell domains and calls parmetis.
  60 Foam::label Foam::parMetisDecomp::decompose
  61 (
  62     Field<int>& xadj,
  63     Field<int>& adjncy,
  64     const pointField& cellCentres,
  65     Field<int>& cellWeights,
  66     Field<int>& faceWeights,
  67     const List<int>& options,
  68     List<int>& finalDecomp
  69 )
  70 {
  71     // C style numbering
  72     int numFlag = 0;
  73
  74     // Number of dimensions
  75     int nDims = 3;
  76
  77
  78     if (cellCentres.size() != xadj.size()-1)
  79     {
  80         FatalErrorIn("parMetisDecomp::decompose(..)")
  81             << "cellCentres:" << cellCentres.size()
  82             << " xadj:" << xadj.size()
  83             << abort(FatalError);
  84     }
  85
  86
  87     // Get number of cells on all processors
  88     List<int> nLocalCells(Pstream::nProcs());
  89     nLocalCells[Pstream::myProcNo()] = xadj.size()-1;
  90     Pstream::gatherList(nLocalCells);
  91     Pstream::scatterList(nLocalCells);
  92
  93     // Get cell offsets.
  94     List<int> cellOffsets(Pstream::nProcs()+1);
  95     int nGlobalCells = 0;
  96     forAll(nLocalCells, procI)
  97     {
  98         cellOffsets[procI] = nGlobalCells;
  99         nGlobalCells += nLocalCells[procI];
 100     }
 101     cellOffsets[Pstream::nProcs()] = nGlobalCells;
 102
 103     // Convert pointField into the data type parMetis expects (float or double)
 104     Field<real_t> xyz(3*cellCentres.size());
 105     int compI = 0;
 106     forAll(cellCentres, cellI)
 107     {
 108         const point& cc = cellCentres[cellI];
 109         xyz[compI++] = float(cc.x());
 110         xyz[compI++] = float(cc.y());
 111         xyz[compI++] = float(cc.z());
 112     }
 113
 114     // Make sure every domain has at least one cell
 115     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 116     // (Metis falls over with zero sized domains)
 117     // Trickle cells from processors that have them up to those that
 118     // don't.
 119
 120
 121     // Number of cells to send to the next processor
 122     // (is same as number of cells next processor has to receive)
 123     List<int> nSendCells(Pstream::nProcs(), 0);
 124
 125     for (label procI = nLocalCells.size()-1; procI >=1; procI--)
 126     {
 127         if (nLocalCells[procI]-nSendCells[procI] < 1)
 128         {
 129             nSendCells[procI-1] = nSendCells[procI]-nLocalCells[procI]+1;
 130         }
 131     }
 132
 133     // First receive (so increasing the sizes of all arrays)
 134
 135     if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0)
 136     {
 137         // Receive cells from previous processor
 138         IPstream fromPrevProc(Pstream::blocking, Pstream::myProcNo()-1);
 139
 140         Field<int> prevXadj(fromPrevProc);
 141         Field<int> prevAdjncy(fromPrevProc);
 142         Field<real_t> prevXyz(fromPrevProc);
 143         Field<int> prevCellWeights(fromPrevProc);
 144         Field<int> prevFaceWeights(fromPrevProc);
 145
 146         if (prevXadj.size() != nSendCells[Pstream::myProcNo()-1])
 147         {
 148             FatalErrorIn("parMetisDecomp::decompose(..)")
 149                 << "Expected from processor " << Pstream::myProcNo()-1
 150                 << " connectivity for " << nSendCells[Pstream::myProcNo()-1]
 151                 << " nCells but only received " << prevXadj.size()
 152                 << abort(FatalError);
 153         }
 154
 155         // Insert adjncy
 156         prepend(prevAdjncy, adjncy);
 157         // Adapt offsets and prepend xadj
 158         xadj += prevAdjncy.size();
 159         prepend(prevXadj, xadj);
 160         // Coords
 161         prepend(prevXyz, xyz);
 162         // Weights
 163         prepend(prevCellWeights, cellWeights);
 164         prepend(prevFaceWeights, faceWeights);
 165     }
 166
 167
 168     // Send to my next processor
 169
 170     if (nSendCells[Pstream::myProcNo()] > 0)
 171     {
 172         // Send cells to next processor
 173         OPstream toNextProc(Pstream::blocking, Pstream::myProcNo()+1);
 174
 175         int nCells = nSendCells[Pstream::myProcNo()];
 176         int startCell = xadj.size()-1 - nCells;
 177         int startFace = xadj[startCell];
 178         int nFaces = adjncy.size()-startFace;
 179
 180         // Send for all cell data: last nCells elements
 181         // Send for all face data: last nFaces elements
 182         toNextProc
 183             << Field<int>::subField(xadj, nCells, startCell)-startFace
 184             << Field<int>::subField(adjncy, nFaces, startFace)
 185             << SubField<real_t>(xyz, nDims*nCells, nDims*startCell)
 186             <<
 187             (
 188                 cellWeights.size()
 189               ? static_cast<const Field<int>&>
 190                 (
 191                     Field<int>::subField(cellWeights, nCells, startCell)
 192                 )
 193               : Field<int>(0)
 194             )
 195             <<
 196             (
 197                 faceWeights.size()
 198               ? static_cast<const Field<int>&>
 199                 (
 200                     Field<int>::subField(faceWeights, nFaces, startFace)
 201                 )
 202               : Field<int>(0)
 203             );
 204
 205         // Remove data that has been sent
 206         if (faceWeights.size())
 207         {
 208             faceWeights.setSize(faceWeights.size()-nFaces);
 209         }
 210         if (cellWeights.size())
 211         {
 212             cellWeights.setSize(cellWeights.size()-nCells);
 213         }
 214         xyz.setSize(xyz.size()-nDims*nCells);
 215         adjncy.setSize(adjncy.size()-nFaces);
 216         xadj.setSize(xadj.size() - nCells);
 217     }
 218
 219
 220
 221     // Adapt number of cells
 222     forAll(nSendCells, procI)
 223     {
 224         // Sent cells
 225         nLocalCells[procI] -= nSendCells[procI];
 226
 227         if (procI >= 1)
 228         {
 229             // Received cells
 230             nLocalCells[procI] += nSendCells[procI-1];
 231         }
 232     }
 233     // Adapt cellOffsets
 234     nGlobalCells = 0;
 235     forAll(nLocalCells, procI)
 236     {
 237         cellOffsets[procI] = nGlobalCells;
 238         nGlobalCells += nLocalCells[procI];
 239     }
 240
 241
 242     if (nLocalCells[Pstream::myProcNo()] != (xadj.size()-1))
 243     {
 244         FatalErrorIn("parMetisDecomp::decompose(..)")
 245             << "Have connectivity for " << xadj.size()-1
 246             << " cells but nLocalCells:" << nLocalCells[Pstream::myProcNo()]
 247             << abort(FatalError);
 248     }
 249
 250     // Weight info
 251     int wgtFlag = 0;
 252     int* vwgtPtr = NULL;
 253     int* adjwgtPtr = NULL;
 254
 255     if (cellWeights.size())
 256     {
 257         vwgtPtr = cellWeights.begin();
 258         wgtFlag += 2;       // Weights on vertices
 259     }
 260     if (faceWeights.size())
 261     {
 262         adjwgtPtr = faceWeights.begin();
 263         wgtFlag += 1;       // Weights on edges
 264     }
 265
 266
 267     // Number of weights or balance constraints
 268     int nCon = 1;
 269     // Per processor, per constraint the weight
 270     Field<real_t> tpwgts(nCon*nProcessors_, 1./nProcessors_);
 271     // Imbalance tolerance
 272     Field<real_t> ubvec(nCon, 1.02);
 273     if (nProcessors_ == 1)
 274     {
 275         // If only one processor there is no imbalance.
 276         ubvec[0] = 1;
 277     }
 278
 279     MPI_Comm comm = MPI_COMM_WORLD;
 280
 281     // output: cell -> processor addressing
 282     finalDecomp.setSize(nLocalCells[Pstream::myProcNo()]);
 283
 284     // output: number of cut edges
 285     int edgeCut = 0;
 286
 287
 288     ParMETIS_V3_PartGeomKway
 289     (
 290         cellOffsets.begin(),    // vtxDist
 291         xadj.begin(),
 292         adjncy.begin(),
 293         vwgtPtr,                // vertexweights
 294         adjwgtPtr,              // edgeweights
 295         &wgtFlag,
 296         &numFlag,
 297         &nDims,
 298         xyz.begin(),
 299         &nCon,
 300         &nProcessors_,          // nParts
 301         tpwgts.begin(),
 302         ubvec.begin(),
 303         const_cast<List<int>&>(options).begin(),
 304         &edgeCut,
 305         finalDecomp.begin(),
 306         &comm
 307     );
 308
 309
 310     // If we sent cells across make sure we undo it
 311     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 312
 313     // Receive back from next processor if I sent something
 314     if (nSendCells[Pstream::myProcNo()] > 0)
 315     {
 316         IPstream fromNextProc(Pstream::blocking, Pstream::myProcNo()+1);
 317
 318         List<int> nextFinalDecomp(fromNextProc);
 319
 320         if (nextFinalDecomp.size() != nSendCells[Pstream::myProcNo()])
 321         {
 322             FatalErrorIn("parMetisDecomp::decompose(..)")
 323                 << "Expected from processor " << Pstream::myProcNo()+1
 324                 << " decomposition for " << nSendCells[Pstream::myProcNo()]
 325                 << " nCells but only received " << nextFinalDecomp.size()
 326                 << abort(FatalError);
 327         }
 328
 329         append(nextFinalDecomp, finalDecomp);
 330     }
 331
 332     // Send back to previous processor.
 333     if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0)
 334     {
 335         OPstream toPrevProc(Pstream::blocking, Pstream::myProcNo()-1);
 336
 337         int nToPrevious = nSendCells[Pstream::myProcNo()-1];
 338
 339         toPrevProc <<
 340             SubList<int>
 341             (
 342                 finalDecomp,
 343                 nToPrevious,
 344                 finalDecomp.size()-nToPrevious
 345             );
 346
 347         // Remove locally what has been sent
 348         finalDecomp.setSize(finalDecomp.size()-nToPrevious);
 349     }
 350
 351     return edgeCut;
 352 }
 353
 354
 355 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 356
 357 Foam::parMetisDecomp::parMetisDecomp
 358 (
 359     const dictionary& decompositionDict,
 360     const polyMesh& mesh
 361 )
 362 :
 363     decompositionMethod(decompositionDict),
 364     mesh_(mesh)
 365 {}
 366
 367
 368 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 369
 370 Foam::labelList Foam::parMetisDecomp::decompose
 371 (
 372     const pointField& cc,
 373     const scalarField& cWeights
 374 )
 375 {
 376     if (cc.size() != mesh_.nCells())
 377     {
 378         FatalErrorIn
 379         (
 380             "parMetisDecomp::decompose"
 381             "(const pointField&, const scalarField&)"
 382         )   << "Can use this decomposition method only for the whole mesh"
 383             << endl
 384             << "and supply one coordinate (cellCentre) for every cell." << endl
 385             << "The number of coordinates " << cc.size() << endl
 386             << "The number of cells in the mesh " << mesh_.nCells()
 387             << exit(FatalError);
 388     }
 389
 390     // For running sequential ...
 391     if (Pstream::nProcs() <= 1)
 392     {
 393         return metisDecomp(decompositionDict_, mesh_).decompose
 394         (
 395             cc,
 396             cWeights
 397         );
 398     }
 399
 400
 401     // Connections
 402     Field<int> adjncy;
 403     // Offsets into adjncy
 404     Field<int> xadj;
 405     calcDistributedCSR
 406     (
 407         mesh_,
 408         adjncy,
 409         xadj
 410     );
 411
 412
 413     // decomposition options. 0 = use defaults
 414     List<int> options(3, 0);
 415     //options[0] = 1;     // don't use defaults but use values below
 416     //options[1] = -1;    // full debug info
 417     //options[2] = 15;    // random number seed
 418
 419     // cell weights (so on the vertices of the dual)
 420     Field<int> cellWeights;
 421
 422     // face weights (so on the edges of the dual)
 423     Field<int> faceWeights;
 424
 425
 426     // Check for externally provided cellweights and if so initialise weights
 427     scalar minWeights = gMin(cWeights);
 428     if (cWeights.size() > 0)
 429     {
 430         if (minWeights <= 0)
 431         {
 432             WarningIn
 433             (
 434                 "metisDecomp::decompose"
 435                 "(const pointField&, const scalarField&)"
 436             )   << "Illegal minimum weight " << minWeights
 437                 << endl;
 438         }
 439
 440         if (cWeights.size() != mesh_.nCells())
 441         {
 442             FatalErrorIn
 443             (
 444                 "parMetisDecomp::decompose"
 445                 "(const pointField&, const scalarField&)"
 446             )   << "Number of cell weights " << cWeights.size()
 447                 << " does not equal number of cells " << mesh_.nCells()
 448                 << exit(FatalError);
 449         }
 450
 451         // Convert to integers.
 452         cellWeights.setSize(cWeights.size());
 453         forAll(cellWeights, i)
 454         {
 455             cellWeights[i] = int(cWeights[i]/minWeights);
 456         }
 457     }
 458
 459
 460     // Check for user supplied weights and decomp options
 461     if (decompositionDict_.found("metisCoeffs"))
 462     {
 463         const dictionary& metisCoeffs =
 464             decompositionDict_.subDict("metisCoeffs");
 465         word weightsFile;
 466
 467         if (metisCoeffs.readIfPresent("cellWeightsFile", weightsFile))
 468         {
 469             Info<< "parMetisDecomp : Using cell-based weights read from "
 470                 << weightsFile << endl;
 471
 472             labelIOField cellIOWeights
 473             (
 474                 IOobject
 475                 (
 476                     weightsFile,
 477                     mesh_.time().timeName(),
 478                     mesh_,
 479                     IOobject::MUST_READ,
 480                     IOobject::AUTO_WRITE
 481                 )
 482             );
 483             cellWeights.transfer(cellIOWeights);
 484
 485             if (cellWeights.size() != mesh_.nCells())
 486             {
 487                 FatalErrorIn
 488                 (
 489                     "parMetisDecomp::decompose"
 490                     "(const pointField&, const scalarField&)"
 491                 )   << "Number of cell weights " << cellWeights.size()
 492                     << " read from " << cellIOWeights.objectPath()
 493                     << " does not equal number of cells " << mesh_.nCells()
 494                     << exit(FatalError);
 495             }
 496         }
 497
 498         if (metisCoeffs.readIfPresent("faceWeightsFile", weightsFile))
 499         {
 500             Info<< "parMetisDecomp : Using face-based weights read from "
 501                 << weightsFile << endl;
 502
 503             labelIOField weights
 504             (
 505                 IOobject
 506                 (
 507                     weightsFile,
 508                     mesh_.time().timeName(),
 509                     mesh_,
 510                     IOobject::MUST_READ,
 511                     IOobject::AUTO_WRITE
 512                 )
 513             );
 514
 515             if (weights.size() != mesh_.nFaces())
 516             {
 517                 FatalErrorIn
 518                 (
 519                     "parMetisDecomp::decompose"
 520                     "(const pointField&, const scalarField&)"
 521                 )   << "Number of face weights " << weights.size()
 522                     << " does not equal number of internal and boundary faces "
 523                     << mesh_.nFaces()
 524                     << exit(FatalError);
 525             }
 526
 527             faceWeights.setSize(adjncy.size());
 528
 529             // Assume symmetric weights. Keep same ordering as adjncy.
 530             List<int> nFacesPerCell(mesh_.nCells(), 0);
 531
 532             // Handle internal faces
 533             for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 534             {
 535                 label w = weights[faceI];
 536
 537                 label own = mesh_.faceOwner()[faceI];
 538                 label nei = mesh_.faceNeighbour()[faceI];
 539
 540                 faceWeights[xadj[own] + nFacesPerCell[own]++] = w;
 541                 faceWeights[xadj[nei] + nFacesPerCell[nei]++] = w;
 542             }
 543             // Coupled boundary faces
 544             const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 545
 546             forAll(patches, patchI)
 547             {
 548                 const polyPatch& pp = patches[patchI];
 549
 550                 if (pp.coupled())
 551                 {
 552                     label faceI = pp.start();
 553
 554                     forAll(pp, i)
 555                     {
 556                         label w = weights[faceI];
 557                         label own = mesh_.faceOwner()[faceI];
 558                         faceWeights[xadj[own] + nFacesPerCell[own]++] = w;
 559                         faceI++;
 560                     }
 561                 }
 562             }
 563         }
 564
 565         if (metisCoeffs.readIfPresent("options", options))
 566         {
 567             Info<< "Using Metis options     " << options
 568                 << nl << endl;
 569
 570             if (options.size() != 3)
 571             {
 572                 FatalErrorIn
 573                 (
 574                     "parMetisDecomp::decompose"
 575                     "(const pointField&, const scalarField&)"
 576                 )   << "Number of options " << options.size()
 577                     << " should be three." << exit(FatalError);
 578             }
 579         }
 580     }
 581
 582
 583     // Do actual decomposition
 584     List<int> finalDecomp;
 585     decompose
 586     (
 587         xadj,
 588         adjncy,
 589         cc,
 590         cellWeights,
 591         faceWeights,
 592         options,
 593         finalDecomp
 594     );
 595
 596     // Copy back to labelList
 597     labelList decomp(finalDecomp.size());
 598     forAll(decomp, i)
 599     {
 600         decomp[i] = finalDecomp[i];
 601     }
 602     return decomp;
 603 }
 604
 605
 606 Foam::labelList Foam::parMetisDecomp::decompose
 607 (
 608     const labelList& cellToRegion,
 609     const pointField& regionPoints,
 610     const scalarField& regionWeights
 611 )
 612 {
 613     const labelList& faceOwner = mesh_.faceOwner();
 614     const labelList& faceNeighbour = mesh_.faceNeighbour();
 615     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 616
 617     if (cellToRegion.size() != mesh_.nCells())
 618     {
 619         FatalErrorIn
 620         (
 621             "parMetisDecomp::decompose(const labelList&, const pointField&)"
 622         )   << "Size of cell-to-coarse map " << cellToRegion.size()
 623             << " differs from number of cells in mesh " << mesh_.nCells()
 624             << exit(FatalError);
 625     }
 626
 627
 628     // Global region numbering engine
 629     globalIndex globalRegions(regionPoints.size());
 630
 631
 632     // Get renumbered owner region on other side of coupled faces
 633     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 634
 635     List<int> globalNeighbour(mesh_.nFaces()-mesh_.nInternalFaces());
 636
 637     forAll(patches, patchI)
 638     {
 639         const polyPatch& pp = patches[patchI];
 640
 641         if (pp.coupled())
 642         {
 643             label faceI = pp.start();
 644             label bFaceI = pp.start() - mesh_.nInternalFaces();
 645
 646             forAll(pp, i)
 647             {
 648                 label ownRegion = cellToRegion[faceOwner[faceI]];
 649                 globalNeighbour[bFaceI++] = globalRegions.toGlobal(ownRegion);
 650                 faceI++;
 651             }
 652         }
 653     }
 654
 655     // Get the cell on the other side of coupled patches
 656     syncTools::swapBoundaryFaceList(mesh_, globalNeighbour, false);
 657
 658
 659     // Get globalCellCells on coarse mesh
 660     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 661
 662     labelListList globalRegionRegions;
 663     {
 664         List<DynamicList<label> > dynRegionRegions(regionPoints.size());
 665
 666         // Internal faces first
 667         forAll(faceNeighbour, faceI)
 668         {
 669             label ownRegion = cellToRegion[faceOwner[faceI]];
 670             label neiRegion = cellToRegion[faceNeighbour[faceI]];
 671
 672             if (ownRegion != neiRegion)
 673             {
 674                 label globalOwn = globalRegions.toGlobal(ownRegion);
 675                 label globalNei = globalRegions.toGlobal(neiRegion);
 676
 677                 if (findIndex(dynRegionRegions[ownRegion], globalNei) == -1)
 678                 {
 679                     dynRegionRegions[ownRegion].append(globalNei);
 680                 }
 681                 if (findIndex(dynRegionRegions[neiRegion], globalOwn) == -1)
 682                 {
 683                     dynRegionRegions[neiRegion].append(globalOwn);
 684                 }
 685             }
 686         }
 687
 688         // Coupled boundary faces
 689         forAll(patches, patchI)
 690         {
 691             const polyPatch& pp = patches[patchI];
 692
 693             if (pp.coupled())
 694             {
 695                 label faceI = pp.start();
 696                 label bFaceI = pp.start() - mesh_.nInternalFaces();
 697
 698                 forAll(pp, i)
 699                 {
 700                     label ownRegion = cellToRegion[faceOwner[faceI]];
 701                     label globalNei = globalNeighbour[bFaceI++];
 702                     faceI++;
 703
 704                     if
 705                     (
 706                         findIndex(dynRegionRegions[ownRegion], globalNei)
 707                      == -1
 708                     )
 709                     {
 710                         dynRegionRegions[ownRegion].append(globalNei);
 711                     }
 712                 }
 713             }
 714         }
 715
 716         globalRegionRegions.setSize(dynRegionRegions.size());
 717         forAll(dynRegionRegions, i)
 718         {
 719             globalRegionRegions[i].transfer(dynRegionRegions[i]);
 720         }
 721     }
 722
 723     labelList regionDecomp
 724     (
 725         decompose
 726         (
 727             globalRegionRegions,
 728             regionPoints,
 729             regionWeights
 730         )
 731     );
 732
 733     // Rework back into decomposition for original mesh_
 734     labelList cellDistribution(cellToRegion.size());
 735
 736     forAll(cellDistribution, cellI)
 737     {
 738         cellDistribution[cellI] = regionDecomp[cellToRegion[cellI]];
 739     }
 740
 741     return cellDistribution;
 742 }
 743
 744
 745 Foam::labelList Foam::parMetisDecomp::decompose
 746 (
 747     const labelListList& globalCellCells,
 748     const pointField& cellCentres,
 749     const scalarField& cWeights
 750 )
 751 {
 752     if (cellCentres.size() != globalCellCells.size())
 753     {
 754         FatalErrorIn
 755         (
 756             "parMetisDecomp::decompose(const labelListList&"
 757             ", const pointField&, const scalarField&)"
 758         )   << "Inconsistent number of cells (" << globalCellCells.size()
 759             << ") and number of cell centres (" << cellCentres.size()
 760             << ") or weights (" << cWeights.size() << ")." << exit(FatalError);
 761     }
 762
 763     // For running sequential ...
 764     if (Pstream::nProcs() <= 1)
 765     {
 766         return metisDecomp(decompositionDict_, mesh_).decompose
 767         (
 768             globalCellCells,
 769             cellCentres,
 770             cWeights
 771         );
 772     }
 773
 774
 775     // Make Metis Distributed CSR (Compressed Storage Format) storage
 776
 777     // Connections
 778     Field<int> adjncy;
 779
 780     // Offsets into adjncy
 781     Field<int> xadj;
 782
 783     calcCSR(globalCellCells, adjncy, xadj);
 784
 785     // decomposition options. 0 = use defaults
 786     List<int> options(3, 0);
 787     //options[0] = 1;     // don't use defaults but use values below
 788     //options[1] = -1;    // full debug info
 789     //options[2] = 15;    // random number seed
 790
 791     // cell weights (so on the vertices of the dual)
 792     Field<int> cellWeights;
 793
 794     // face weights (so on the edges of the dual)
 795     Field<int> faceWeights;
 796
 797
 798     // Check for externally provided cellweights and if so initialise weights
 799     scalar minWeights = gMin(cWeights);
 800     if (cWeights.size() > 0)
 801     {
 802         if (minWeights <= 0)
 803         {
 804             WarningIn
 805             (
 806                 "parMetisDecomp::decompose(const labelListList&"
 807                 ", const pointField&, const scalarField&)"
 808             )   << "Illegal minimum weight " << minWeights
 809                 << endl;
 810         }
 811
 812         if (cWeights.size() != globalCellCells.size())
 813         {
 814             FatalErrorIn
 815             (
 816                 "parMetisDecomp::decompose(const labelListList&"
 817                 ", const pointField&, const scalarField&)"
 818             )   << "Number of cell weights " << cWeights.size()
 819                 << " does not equal number of cells " << globalCellCells.size()
 820                 << exit(FatalError);
 821         }
 822
 823         // Convert to integers.
 824         cellWeights.setSize(cWeights.size());
 825         forAll(cellWeights, i)
 826         {
 827             cellWeights[i] = int(cWeights[i]/minWeights);
 828         }
 829     }
 830
 831
 832     // Check for user supplied weights and decomp options
 833     if (decompositionDict_.found("metisCoeffs"))
 834     {
 835         const dictionary& metisCoeffs =
 836             decompositionDict_.subDict("metisCoeffs");
 837
 838         if (metisCoeffs.readIfPresent("options", options))
 839         {
 840             Info<< "Using Metis options     " << options
 841                 << nl << endl;
 842
 843             if (options.size() != 3)
 844             {
 845                 FatalErrorIn
 846                 (
 847                     "parMetisDecomp::decompose(const labelListList&"
 848                     ", const pointField&, const scalarField&)"
 849                 )   << "Number of options " << options.size()
 850                     << " should be three." << exit(FatalError);
 851             }
 852         }
 853     }
 854
 855
 856     // Do actual decomposition
 857     List<int> finalDecomp;
 858     decompose
 859     (
 860         xadj,
 861         adjncy,
 862         cellCentres,
 863         cellWeights,
 864         faceWeights,
 865         options,
 866         finalDecomp
 867     );
 868
 869     // Copy back to labelList
 870     labelList decomp(finalDecomp.size());
 871     forAll(decomp, i)
 872     {
 873         decomp[i] = finalDecomp[i];
 874     }
 875     return decomp;
 876 }
 877
 878
 879 // ************************************************************************* //