applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 "domainDecomposition.H"
  28 #include "Time.H"
  29 #include "dictionary.H"
  30 #include "labelIOList.H"
  31 #include "processorPolyPatch.H"
  32 #include "fvMesh.H"
  33 #include "OSspecific.H"
  34 #include "Map.H"
  35 #include "globalMeshData.H"
  36 #include "DynamicList.H"
  37
  38 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  39
  40 domainDecomposition::domainDecomposition(const IOobject& io)
  41 :
  42     fvMesh(io),
  43     decompositionDict_
  44     (
  45         IOobject
  46         (
  47             "decomposeParDict",
  48             time().system(),
  49             *this,
  50             IOobject::MUST_READ,
  51             IOobject::NO_WRITE
  52         )
  53     ),
  54     nProcs_(readInt(decompositionDict_.lookup("numberOfSubdomains"))),
  55     distributed_(false),
  56     cellToProc_(nCells()),
  57     procPointAddressing_(nProcs_),
  58     procFaceAddressing_(nProcs_),
  59     nInternalProcFaces_(nProcs_),
  60     nLiveProcFaces_(nProcs_),
  61     procCellAddressing_(nProcs_),
  62     procBoundaryAddressing_(nProcs_),
  63     procPatchSize_(nProcs_),
  64     procPatchStartIndex_(nProcs_),
  65     procNeighbourProcessors_(nProcs_),
  66     procProcessorPatchSize_(nProcs_),
  67     procProcessorPatchStartIndex_(nProcs_),
  68     globallySharedPoints_(0),
  69     cyclicParallel_(false)
  70 {
  71     if (decompositionDict_.found("distributed"))
  72     {
  73         distributed_ = Switch(decompositionDict_.lookup("distributed"));
  74     }
  75 }
  76
  77
  78 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  79
  80 domainDecomposition::~domainDecomposition()
  81 {}
  82
  83
  84 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  85
  86 bool domainDecomposition::writeDecomposition()
  87 {
  88     Info<< "\nConstructing processor meshes" << endl;
  89
  90     // Make a lookup map for globally shared points
  91     Map<label> sharedPointLookup(2*globallySharedPoints_.size());
  92
  93     forAll (globallySharedPoints_, pointi)
  94     {
  95         sharedPointLookup.insert(globallySharedPoints_[pointi], pointi);
  96     }
  97
  98
  99     // Mark point/faces/cells that are in zones.  Bad coding - removed
 100     // HJ, 31/Mar/2009
 101
 102     label totProcFaces = 0;
 103     label maxProcPatches = 0;
 104     label maxProcFaces = 0;
 105
 106
 107     // Write out the meshes
 108     for (label procI = 0; procI < nProcs_; procI++)
 109     {
 110         // Create processor points
 111         const labelList& curPointLabels = procPointAddressing_[procI];
 112
 113         // Access list of all points in the mesh.  HJ, 27/Mar/2009
 114         const pointField& meshPoints = allPoints();
 115
 116         labelList pointLookup(meshPoints.size(), -1);
 117
 118         pointField procPoints(curPointLabels.size());
 119
 120         forAll (curPointLabels, pointi)
 121         {
 122             procPoints[pointi] = meshPoints[curPointLabels[pointi]];
 123
 124             pointLookup[curPointLabels[pointi]] = pointi;
 125         }
 126
 127         // Create processor faces
 128         const labelList& curFaceLabels = procFaceAddressing_[procI];
 129
 130         // Access list of all faces in the mesh.  HJ, 27/Mar/2009
 131         const faceList& meshFaces = allFaces();
 132
 133         labelList faceLookup(meshFaces.size(), -1);
 134
 135         faceList procFaces(curFaceLabels.size());
 136
 137         forAll (curFaceLabels, facei)
 138         {
 139             // Mark the original face as used
 140             // Remember to decrement the index by one (turning index)
 141             // HJ, 5/Dec/2001
 142             label curF = mag(curFaceLabels[facei]) - 1;
 143
 144             faceLookup[curF] = facei;
 145
 146             // get the original face
 147             labelList origFaceLabels;
 148
 149             if (curFaceLabels[facei] >= 0)
 150             {
 151                 // face not turned
 152                 origFaceLabels = meshFaces[curF];
 153             }
 154             else
 155             {
 156                 origFaceLabels = meshFaces[curF].reverseFace();
 157             }
 158
 159             // translate face labels into local point list
 160             face& procFaceLabels = procFaces[facei];
 161
 162             procFaceLabels.setSize(origFaceLabels.size());
 163
 164             forAll (origFaceLabels, pointi)
 165             {
 166                 procFaceLabels[pointi] = pointLookup[origFaceLabels[pointi]];
 167             }
 168         }
 169
 170         // Create cell lookup
 171         labelList cellLookup(nCells(), -1);
 172         const labelList& curCellLabels = procCellAddressing_[procI];
 173
 174         forAll (curCellLabels, cellI)
 175         {
 176             cellLookup[curCellLabels[cellI]] = cellI;
 177         }
 178
 179         // Get complete owner-neighour addressing in the mesh
 180         const labelList& own = faceOwner();
 181         const labelList& nei = faceNeighbour();
 182
 183         // Calculate owner and neighbour list
 184         // Owner list is sized to number of live faces
 185         // Neighbour list is sized to number of internal faces
 186
 187         labelList procOwner(nLiveProcFaces_[procI]);
 188
 189         // Note: loop over owner, not all faces: sizes are different
 190         forAll (procOwner, faceI)
 191         {
 192             // Remember to decrement the index by one (turning index)
 193             // HJ, 28/Mar/2009
 194             label curF = mag(curFaceLabels[faceI]) - 1;
 195
 196             if (curFaceLabels[faceI] >= 0)
 197             {
 198                 procOwner[faceI] = cellLookup[own[curF]];
 199             }
 200             else
 201             {
 202                 procOwner[faceI] = cellLookup[nei[curF]];
 203             }
 204         }
 205
 206         labelList procNeighbour(nInternalProcFaces_[procI]);
 207
 208         // Note: loop over neighbour, not all faces: sizes are different
 209         forAll (procNeighbour, faceI)
 210         {
 211             // Remember to decrement the index by one (turning index)
 212             // HJ, 28/Mar/2009
 213             label curF = mag(curFaceLabels[faceI]) - 1;
 214
 215             if (curFaceLabels[faceI] >= 0)
 216             {
 217                 procNeighbour[faceI] = cellLookup[nei[curF]];
 218             }
 219             else
 220             {
 221                 procNeighbour[faceI] = cellLookup[own[curF]];
 222             }
 223         }
 224
 225         // Create processor cells.  No longer needed: using owner and neighbour
 226         // HJ, 28/Mar/2009
 227 //         const cellList& meshCells = cells();
 228
 229 //         cellList procCells(curCellLabels.size());
 230
 231 //         forAll (curCellLabels, cellI)
 232 //         {
 233 //             const labelList& origCellLabels = meshCells[curCellLabels[cellI]];
 234
 235 //             cell& curCell = procCells[cellI];
 236
 237 //             curCell.setSize(origCellLabels.size());
 238
 239 //             forAll (origCellLabels, cellFaceI)
 240 //             {
 241 //                 curCell[cellFaceI] = faceLookup[origCellLabels[cellFaceI]];
 242 //             }
 243 //         }
 244
 245         // Create processor mesh without a boundary
 246
 247         fileName processorCasePath
 248         (
 249             time().caseName()/fileName(word("processor") + Foam::name(procI))
 250         );
 251
 252         // make the processor directory
 253         mkDir(time().rootPath()/processorCasePath);
 254
 255         // create a database
 256         Time processorDb
 257         (
 258             Time::controlDictName,
 259             time().rootPath(),
 260             processorCasePath,
 261             "system",
 262             "constant",
 263             true
 264         );
 265
 266         // Create the mesh
 267         polyMesh procMesh
 268         (
 269             IOobject
 270             (
 271                 this->polyMesh::name(),  // region name of undecomposed mesh
 272                 pointsInstance(),
 273                 processorDb
 274             ),
 275             xferMove(procPoints),
 276             xferMove(procFaces),
 277             xferMove(procOwner),
 278             xferMove(procNeighbour),
 279             false          // Do not sync par
 280 //   xferMove(procCells)   // Old-fashioned mesh creation using cells.
 281                            // Deprecated: using face owner/neighbour
 282                            // HJ, 30/Mar/2009
 283         );
 284
 285         // Create processor boundary patches
 286         const labelList& curBoundaryAddressing = procBoundaryAddressing_[procI];
 287
 288         const labelList& curPatchSizes = procPatchSize_[procI];
 289
 290         const labelList& curPatchStarts = procPatchStartIndex_[procI];
 291
 292         const labelList& curNeighbourProcessors =
 293             procNeighbourProcessors_[procI];
 294
 295         const labelList& curProcessorPatchSizes =
 296             procProcessorPatchSize_[procI];
 297
 298         const labelList& curProcessorPatchStarts =
 299             procProcessorPatchStartIndex_[procI];
 300
 301         const polyPatchList& meshPatches = boundaryMesh();
 302
 303         List<polyPatch*> procPatches
 304         (
 305             curPatchSizes.size()
 306           + curProcessorPatchSizes.size(),
 307             reinterpret_cast<polyPatch*>(0)
 308         );
 309
 310         label nPatches = 0;
 311
 312         forAll (curPatchSizes, patchi)
 313         {
 314             procPatches[nPatches] =
 315                 meshPatches[curBoundaryAddressing[patchi]].clone
 316                 (
 317                     procMesh.boundaryMesh(),
 318                     nPatches,
 319                     curPatchSizes[patchi],
 320                     curPatchStarts[patchi]
 321                 ).ptr();
 322
 323             nPatches++;
 324         }
 325
 326         forAll (curProcessorPatchSizes, procPatchI)
 327         {
 328             procPatches[nPatches] =
 329                 new processorPolyPatch
 330                 (
 331                     word("procBoundary") + Foam::name(procI)
 332                   + word("to")
 333                   + Foam::name(curNeighbourProcessors[procPatchI]),
 334                     curProcessorPatchSizes[procPatchI],
 335                     curProcessorPatchStarts[procPatchI],
 336                     nPatches,
 337                     procMesh.boundaryMesh(),
 338                     procI,
 339                     curNeighbourProcessors[procPatchI]
 340             );
 341
 342             nPatches++;
 343         }
 344
 345         // Add boundary patches
 346         procMesh.addPatches(procPatches);
 347
 348         // Create and add zones
 349
 350         // Note:
 351         // This coding was all wrong, as each point/face/cell may only
 352         //  belong to a single zone.
 353         // Additionally, ordering of points/faces/cells in the processor mesh
 354         // needs to match the ordering in global mesh zones.  Full rewrite.
 355         // HJ, 30/Mar/2009
 356
 357         // Create zones if needed
 358         if
 359         (
 360             pointZones().size() > 0
 361          || faceZones().size() > 0
 362          || cellZones().size() > 0
 363         )
 364         {
 365             // Point zones
 366             List<pointZone*> procPz(pointZones().size());
 367
 368             {
 369                 const pointZoneMesh& pz = pointZones();
 370
 371                 // Go through all the zoned points and find out if they
 372                 // belong to a processor.  If so, add it to the zone as
 373                 // necessary
 374                 forAll (pz, zoneI)
 375                 {
 376                     const labelList& zonePoints = pz[zoneI];
 377
 378                     labelList procZonePoints(zonePoints.size());
 379                     label nZonePoints = 0;
 380
 381                     forAll (zonePoints, pointI)
 382                     {
 383                         const label localIndex =
 384                             pointLookup[zonePoints[pointI]];
 385
 386                         if (localIndex >= 0)
 387                         {
 388                             // Point live on processor: add to zone
 389                             procZonePoints[nZonePoints] = localIndex;
 390                             nZonePoints++;
 391                         }
 392                     }
 393
 394                     // Add the zone
 395                     procZonePoints.setSize(nZonePoints);
 396
 397                     procPz[zoneI] = new pointZone
 398                     (
 399                         pz[zoneI].name(),
 400                         procZonePoints,
 401                         zoneI,
 402                         procMesh.pointZones()
 403                     );
 404                 }
 405             }
 406
 407
 408             // Face zones
 409             List<faceZone*> procFz(faceZones().size());
 410
 411             {
 412                 const faceZoneMesh& fz = faceZones();
 413
 414                 forAll (fz, zoneI)
 415                 {
 416                     const labelList& zoneFaces = fz[zoneI];
 417                     const boolList& flipMap = fz[zoneI].flipMap();
 418
 419                     // Go through all the zoned faces and find out if they
 420                     // belong to a processor.  If so, add it to the zone as
 421                     // necessary
 422
 423                     labelList procZoneFaces(zoneFaces.size());
 424                     boolList procZoneFaceFlips(zoneFaces.size());
 425                     label nZoneFaces = 0;
 426
 427                     forAll (zoneFaces, faceI)
 428                     {
 429                         const label localIndex = faceLookup[zoneFaces[faceI]];
 430
 431                         if (localIndex >= 0)
 432                         {
 433                             // Face is present on the processor
 434
 435                             // Add the face to the zone
 436                             procZoneFaces[nZoneFaces] = localIndex;
 437
 438                             // Grab the flip
 439                             bool flip = flipMap[faceI];
 440
 441                             if (curFaceLabels[localIndex] < 0)
 442                             {
 443                                 flip = !flip;
 444                             }
 445
 446                             procZoneFaceFlips[nZoneFaces] = flip;
 447                             nZoneFaces++;
 448                         }
 449                     }
 450
 451                     // Add the zone
 452                     procZoneFaces.setSize(nZoneFaces);
 453                     procZoneFaceFlips.setSize(nZoneFaces);
 454
 455                     procFz[zoneI] = new faceZone
 456                     (
 457                         fz[zoneI].name(),
 458                         procZoneFaces,
 459                         procZoneFaceFlips,
 460                         zoneI,
 461                         procMesh.faceZones()
 462                     );
 463                 }
 464             }
 465
 466             // Cell zones
 467             List<cellZone*> procCz(cellZones().size());
 468
 469             {
 470                 const cellZoneMesh& cz = cellZones();
 471
 472                 // Go through all the zoned cells and find out if they
 473                 // belong to a processor.  If so, add it to the zone as
 474                 // necessary
 475
 476                 forAll (cz, zoneI)
 477                 {
 478                     const labelList& zoneCells = cz[zoneI];
 479
 480                     labelList procZoneCells(zoneCells.size());
 481                     label nZoneCells = 0;
 482
 483                     forAll (zoneCells, cellI)
 484                     {
 485                         const label localIndex = cellLookup[zoneCells[cellI]];
 486
 487                         if (localIndex >= 0)
 488                         {
 489                             procZoneCells[nZoneCells] = localIndex;
 490                             nZoneCells++;
 491                         }
 492                     }
 493
 494                     // Add the zone
 495                     procZoneCells.setSize(nZoneCells);
 496
 497                     procCz[zoneI] = new cellZone
 498                     (
 499                         cz[zoneI].name(),
 500                         procZoneCells,
 501                         zoneI,
 502                         procMesh.cellZones()
 503                     );
 504                 }
 505             }
 506
 507             // Add zones
 508             procMesh.addZones(procPz, procFz, procCz);
 509         }
 510
 511         // Set the precision of the points data to 10
 512         IOstream::defaultPrecision(10);
 513
 514         procMesh.write();
 515
 516         Info<< endl
 517             << "Processor " << procI << nl
 518             << "    Number of cells = " << procMesh.nCells()
 519             << endl;
 520
 521         label nBoundaryFaces = 0;
 522         label nProcPatches = 0;
 523         label nProcFaces = 0;
 524
 525         forAll (procMesh.boundaryMesh(), patchi)
 526         {
 527             if
 528             (
 529                 procMesh.boundaryMesh()[patchi].type()
 530              == processorPolyPatch::typeName
 531             )
 532             {
 533                 const processorPolyPatch& ppp =
 534                 refCast<const processorPolyPatch>
 535                 (
 536                     procMesh.boundaryMesh()[patchi]
 537                 );
 538
 539                 Info<< "    Number of faces shared with processor "
 540                     << ppp.neighbProcNo() << " = " << ppp.size() << endl;
 541
 542                 nProcPatches++;
 543                 nProcFaces += ppp.size();
 544             }
 545             else
 546             {
 547                 nBoundaryFaces += procMesh.boundaryMesh()[patchi].size();
 548             }
 549         }
 550
 551         Info<< "    Number of processor patches = " << nProcPatches << nl
 552             << "    Number of processor faces = " << nProcFaces << nl
 553             << "    Number of boundary faces = " << nBoundaryFaces << endl;
 554
 555         totProcFaces += nProcFaces;
 556         maxProcPatches = max(maxProcPatches, nProcPatches);
 557         maxProcFaces = max(maxProcFaces, nProcFaces);
 558
 559         // create and write the addressing information
 560         labelIOList pointProcAddressing
 561         (
 562             IOobject
 563             (
 564                 "pointProcAddressing",
 565                 procMesh.facesInstance(),
 566                 procMesh.meshSubDir,
 567                 procMesh,
 568                 IOobject::NO_READ,
 569                 IOobject::NO_WRITE
 570             ),
 571             procPointAddressing_[procI]
 572         );
 573         pointProcAddressing.write();
 574
 575         labelIOList faceProcAddressing
 576         (
 577             IOobject
 578             (
 579                 "faceProcAddressing",
 580                 procMesh.facesInstance(),
 581                 procMesh.meshSubDir,
 582                 procMesh,
 583                 IOobject::NO_READ,
 584                 IOobject::NO_WRITE
 585             ),
 586             procFaceAddressing_[procI]
 587         );
 588         faceProcAddressing.write();
 589
 590         labelIOList cellProcAddressing
 591         (
 592             IOobject
 593             (
 594                 "cellProcAddressing",
 595                 procMesh.facesInstance(),
 596                 procMesh.meshSubDir,
 597                 procMesh,
 598                 IOobject::NO_READ,
 599                 IOobject::NO_WRITE
 600             ),
 601             procCellAddressing_[procI]
 602         );
 603         cellProcAddressing.write();
 604
 605         labelIOList boundaryProcAddressing
 606         (
 607             IOobject
 608             (
 609                 "boundaryProcAddressing",
 610                 procMesh.facesInstance(),
 611                 procMesh.meshSubDir,
 612                 procMesh,
 613                 IOobject::NO_READ,
 614                 IOobject::NO_WRITE
 615             ),
 616             procBoundaryAddressing_[procI]
 617         );
 618         boundaryProcAddressing.write();
 619     }
 620
 621     Info<< nl
 622         << "Number of processor faces = " << totProcFaces/2 << nl
 623         << "Max number of processor patches = " << maxProcPatches << nl
 624         << "Max number of faces between processors = " << maxProcFaces
 625         << endl;
 626
 627     return true;
 628 }
 629
 630
 631 // ************************************************************************* //