src/decompositionMethods/parMetisDecomp/parMetisDecomp.C

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