src/decompositionAgglomeration/parMetisDecomp/parMetisDecomp.C

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