applications/utilities/parallelProcessing/decomposePar/domainDecomposition.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 "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 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  39
  40 void domainDecomposition::mark
  41 (
  42     const labelList& zoneElems,
  43     const label zoneI,
  44     labelList& elementToZone
  45 )
  46 {
  47     forAll(zoneElems, i)
  48     {
  49         label pointi = zoneElems[i];
  50
  51         if (elementToZone[pointi] == -1)
  52         {
  53             // First occurrence
  54             elementToZone[pointi] = zoneI;
  55         }
  56         else if (elementToZone[pointi] >= 0)
  57         {
  58             // Multiple zones
  59             elementToZone[pointi] = -2;
  60         }
  61     }
  62 }
  63
  64
  65 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  66
  67 // from components
  68 domainDecomposition::domainDecomposition(const IOobject& io)
  69 :
  70     fvMesh(io),
  71     decompositionDict_
  72     (
  73         IOobject
  74         (
  75             "decomposeParDict",
  76             time().system(),
  77             *this,
  78             IOobject::MUST_READ,
  79             IOobject::NO_WRITE
  80         )
  81     ),
  82     nProcs_(readInt(decompositionDict_.lookup("numberOfSubdomains"))),
  83     distributed_(false),
  84     cellToProc_(nCells()),
  85     procPointAddressing_(nProcs_),
  86     procFaceAddressing_(nProcs_),
  87     procCellAddressing_(nProcs_),
  88     procBoundaryAddressing_(nProcs_),
  89     procPatchSize_(nProcs_),
  90     procPatchStartIndex_(nProcs_),
  91     procNeighbourProcessors_(nProcs_),
  92     procProcessorPatchSize_(nProcs_),
  93     procProcessorPatchStartIndex_(nProcs_),
  94     globallySharedPoints_(0),
  95     cyclicParallel_(false)
  96 {
  97     if (decompositionDict_.found("distributed"))
  98     {
  99         Switch distributed(decompositionDict_.lookup("distributed"));
 100         distributed_ = distributed;
 101     }
 102 }
 103
 104
 105 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 106
 107 domainDecomposition::~domainDecomposition()
 108 {}
 109
 110
 111 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 112
 113 bool domainDecomposition::writeDecomposition()
 114 {
 115     Info<< "\nConstructing processor meshes" << endl;
 116
 117     // Make a lookup map for globally shared points
 118     Map<label> sharedPointLookup(2*globallySharedPoints_.size());
 119
 120     forAll (globallySharedPoints_, pointi)
 121     {
 122         sharedPointLookup.insert(globallySharedPoints_[pointi], pointi);
 123     }
 124
 125
 126     // Mark point/faces/cells that are in zones.
 127     // -1   : not in zone
 128     // -2   : in multiple zones
 129     // >= 0 : in single given zone
 130     // This will give direct lookup of elements that are in a single zone
 131     // and we'll only have to revert back to searching through all zones
 132     // for the duplicate elements
 133
 134     // Point zones
 135     labelList pointToZone(points().size(), -1);
 136
 137     forAll(pointZones(), zoneI)
 138     {
 139         mark(pointZones()[zoneI], zoneI, pointToZone);
 140     }
 141
 142     // Face zones
 143     labelList faceToZone(faces().size(), -1);
 144
 145     forAll(faceZones(), zoneI)
 146     {
 147         mark(faceZones()[zoneI], zoneI, faceToZone);
 148     }
 149
 150     // Cell zones
 151     labelList cellToZone(nCells(), -1);
 152
 153     forAll(cellZones(), zoneI)
 154     {
 155         mark(cellZones()[zoneI], zoneI, cellToZone);
 156     }
 157
 158
 159     label totProcFaces = 0;
 160     label maxProcPatches = 0;
 161     label maxProcFaces = 0;
 162
 163
 164     // Write out the meshes
 165     for (label procI = 0; procI < nProcs_; procI++)
 166     {
 167         // Create processor points
 168         const labelList& curPointLabels = procPointAddressing_[procI];
 169
 170         const pointField& meshPoints = points();
 171
 172         labelList pointLookup(nPoints(), -1);
 173
 174         pointField procPoints(curPointLabels.size());
 175
 176         forAll (curPointLabels, pointi)
 177         {
 178             procPoints[pointi] = meshPoints[curPointLabels[pointi]];
 179
 180             pointLookup[curPointLabels[pointi]] = pointi;
 181         }
 182
 183         // Create processor faces
 184         const labelList& curFaceLabels = procFaceAddressing_[procI];
 185
 186         const faceList& meshFaces = faces();
 187
 188         labelList faceLookup(nFaces(), -1);
 189
 190         faceList procFaces(curFaceLabels.size());
 191
 192         forAll (curFaceLabels, facei)
 193         {
 194             // Mark the original face as used
 195             // Remember to decrement the index by one (turning index)
 196             //
 197             label curF = mag(curFaceLabels[facei]) - 1;
 198
 199             faceLookup[curF] = facei;
 200
 201             // get the original face
 202             labelList origFaceLabels;
 203
 204             if (curFaceLabels[facei] >= 0)
 205             {
 206                 // face not turned
 207                 origFaceLabels = meshFaces[curF];
 208             }
 209             else
 210             {
 211                 origFaceLabels = meshFaces[curF].reverseFace();
 212             }
 213
 214             // translate face labels into local point list
 215             face& procFaceLabels = procFaces[facei];
 216
 217             procFaceLabels.setSize(origFaceLabels.size());
 218
 219             forAll (origFaceLabels, pointi)
 220             {
 221                 procFaceLabels[pointi] = pointLookup[origFaceLabels[pointi]];
 222             }
 223         }
 224
 225         // Create processor cells
 226         const labelList& curCellLabels = procCellAddressing_[procI];
 227
 228         const cellList& meshCells = cells();
 229
 230         cellList procCells(curCellLabels.size());
 231
 232         forAll (curCellLabels, celli)
 233         {
 234             const labelList& origCellLabels = meshCells[curCellLabels[celli]];
 235
 236             cell& curCell = procCells[celli];
 237
 238             curCell.setSize(origCellLabels.size());
 239
 240             forAll (origCellLabels, cellFaceI)
 241             {
 242                 curCell[cellFaceI] = faceLookup[origCellLabels[cellFaceI]];
 243             }
 244         }
 245
 246         // Create processor mesh without a boundary
 247
 248         fileName processorCasePath
 249         (
 250             time().caseName()/fileName(word("processor") + Foam::name(procI))
 251         );
 252
 253         // make the processor directory
 254         mkDir(time().rootPath()/processorCasePath);
 255
 256         // create a database
 257         Time processorDb
 258         (
 259             Time::controlDictName,
 260             time().rootPath(),
 261             processorCasePath,
 262             "system",
 263             "constant"
 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(procCells)
 278         );
 279
 280         // Create processor boundary patches
 281         const labelList& curBoundaryAddressing = procBoundaryAddressing_[procI];
 282
 283         const labelList& curPatchSizes = procPatchSize_[procI];
 284
 285         const labelList& curPatchStarts = procPatchStartIndex_[procI];
 286
 287         const labelList& curNeighbourProcessors =
 288             procNeighbourProcessors_[procI];
 289
 290         const labelList& curProcessorPatchSizes =
 291             procProcessorPatchSize_[procI];
 292
 293         const labelList& curProcessorPatchStarts =
 294             procProcessorPatchStartIndex_[procI];
 295
 296         const polyPatchList& meshPatches = boundaryMesh();
 297
 298         List<polyPatch*> procPatches
 299         (
 300             curPatchSizes.size()
 301           + curProcessorPatchSizes.size(),
 302             reinterpret_cast<polyPatch*>(0)
 303         );
 304
 305         label nPatches = 0;
 306
 307         forAll (curPatchSizes, patchi)
 308         {
 309             procPatches[nPatches] =
 310                 meshPatches[curBoundaryAddressing[patchi]].clone
 311                 (
 312                     procMesh.boundaryMesh(),
 313                     nPatches,
 314                     curPatchSizes[patchi],
 315                     curPatchStarts[patchi]
 316                 ).ptr();
 317
 318             nPatches++;
 319         }
 320
 321         forAll (curProcessorPatchSizes, procPatchI)
 322         {
 323             procPatches[nPatches] =
 324                 new processorPolyPatch
 325                 (
 326                     word("procBoundary") + Foam::name(procI)
 327                   + word("to")
 328                   + Foam::name(curNeighbourProcessors[procPatchI]),
 329                     curProcessorPatchSizes[procPatchI],
 330                     curProcessorPatchStarts[procPatchI],
 331                     nPatches,
 332                     procMesh.boundaryMesh(),
 333                     procI,
 334                     curNeighbourProcessors[procPatchI]
 335             );
 336
 337             nPatches++;
 338         }
 339
 340         // Add boundary patches
 341         procMesh.addPatches(procPatches);
 342
 343         // Create and add zones
 344
 345         // Point zones
 346         {
 347             const pointZoneMesh& pz = pointZones();
 348
 349             // Go through all the zoned points and find out if they
 350             // belong to a zone.  If so, add it to the zone as
 351             // necessary
 352             List<DynamicList<label> > zonePoints(pz.size());
 353
 354             // Estimate size
 355             forAll(zonePoints, zoneI)
 356             {
 357                 zonePoints[zoneI].setCapacity(pz[zoneI].size() / nProcs_);
 358             }
 359
 360             // Use the pointToZone map to find out the single zone (if any),
 361             // use slow search only for shared points.
 362             forAll (curPointLabels, pointi)
 363             {
 364                 label curPoint = curPointLabels[pointi];
 365
 366                 label zoneI = pointToZone[curPoint];
 367
 368                 if (zoneI >= 0)
 369                 {
 370                     // Single zone.
 371                     zonePoints[zoneI].append(pointi);
 372                 }
 373                 else if (zoneI == -2)
 374                 {
 375                     // Multiple zones. Lookup.
 376                     forAll(pz, zoneI)
 377                     {
 378                         label index = pz[zoneI].whichPoint(curPoint);
 379
 380                         if (index != -1)
 381                         {
 382                             zonePoints[zoneI].append(pointi);
 383                         }
 384                     }
 385                 }
 386             }
 387
 388             procMesh.pointZones().clearAddressing();
 389             procMesh.pointZones().setSize(zonePoints.size());
 390             forAll(zonePoints, zoneI)
 391             {
 392                 procMesh.pointZones().set
 393                 (
 394                     zoneI,
 395                     pz[zoneI].clone
 396                     (
 397                         procMesh.pointZones(),
 398                         zoneI,
 399                         zonePoints[zoneI].shrink()
 400                     )
 401                 );
 402             }
 403
 404             if (pz.size())
 405             {
 406                 // Force writing on all processors
 407                 procMesh.pointZones().writeOpt() = IOobject::AUTO_WRITE;
 408             }
 409         }
 410
 411         // Face zones
 412         {
 413             const faceZoneMesh& fz = faceZones();
 414
 415             // Go through all the zoned face and find out if they
 416             // belong to a zone.  If so, add it to the zone as
 417             // necessary
 418             List<DynamicList<label> > zoneFaces(fz.size());
 419             List<DynamicList<bool> > zoneFaceFlips(fz.size());
 420
 421             // Estimate size
 422             forAll(zoneFaces, zoneI)
 423             {
 424                 label procSize = fz[zoneI].size() / nProcs_;
 425
 426                 zoneFaces[zoneI].setCapacity(procSize);
 427                 zoneFaceFlips[zoneI].setCapacity(procSize);
 428             }
 429
 430             // Go through all the zoned faces and find out if they
 431             // belong to a zone.  If so, add it to the zone as
 432             // necessary
 433             forAll (curFaceLabels, facei)
 434             {
 435                 // Remember to decrement the index by one (turning index)
 436                 //
 437                 label curF = mag(curFaceLabels[facei]) - 1;
 438
 439                 label zoneI = faceToZone[curF];
 440
 441                 if (zoneI >= 0)
 442                 {
 443                     // Single zone. Add the face
 444                     zoneFaces[zoneI].append(facei);
 445
 446                     label index = fz[zoneI].whichFace(curF);
 447
 448                     bool flip = fz[zoneI].flipMap()[index];
 449
 450                     if (curFaceLabels[facei] < 0)
 451                     {
 452                         flip = !flip;
 453                     }
 454
 455                     zoneFaceFlips[zoneI].append(flip);
 456                 }
 457                 else if (zoneI == -2)
 458                 {
 459                     // Multiple zones. Lookup.
 460                     forAll(fz, zoneI)
 461                     {
 462                         label index = fz[zoneI].whichFace(curF);
 463
 464                         if (index != -1)
 465                         {
 466                             zoneFaces[zoneI].append(facei);
 467
 468                             bool flip = fz[zoneI].flipMap()[index];
 469
 470                             if (curFaceLabels[facei] < 0)
 471                             {
 472                                 flip = !flip;
 473                             }
 474
 475                             zoneFaceFlips[zoneI].append(flip);
 476                         }
 477                     }
 478                 }
 479             }
 480
 481             procMesh.faceZones().clearAddressing();
 482             procMesh.faceZones().setSize(zoneFaces.size());
 483             forAll(zoneFaces, zoneI)
 484             {
 485                 procMesh.faceZones().set
 486                 (
 487                     zoneI,
 488                     fz[zoneI].clone
 489                     (
 490                         zoneFaces[zoneI].shrink(),          // addressing
 491                         zoneFaceFlips[zoneI].shrink(),      // flipmap
 492                         zoneI,
 493                         procMesh.faceZones()
 494                     )
 495                 );
 496             }
 497
 498             if (fz.size())
 499             {
 500                 // Force writing on all processors
 501                 procMesh.faceZones().writeOpt() = IOobject::AUTO_WRITE;
 502             }
 503         }
 504
 505         // Cell zones
 506         {
 507             const cellZoneMesh& cz = cellZones();
 508
 509             // Go through all the zoned cells and find out if they
 510             // belong to a zone.  If so, add it to the zone as
 511             // necessary
 512             List<DynamicList<label> > zoneCells(cz.size());
 513
 514             // Estimate size
 515             forAll(zoneCells, zoneI)
 516             {
 517                 zoneCells[zoneI].setCapacity(cz[zoneI].size() / nProcs_);
 518             }
 519
 520             forAll (curCellLabels, celli)
 521             {
 522                 label curCellI = curCellLabels[celli];
 523
 524                 label zoneI = cellToZone[curCellI];
 525
 526                 if (zoneI >= 0)
 527                 {
 528                     // Single zone.
 529                     zoneCells[zoneI].append(celli);
 530                 }
 531                 else if (zoneI == -2)
 532                 {
 533                     // Multiple zones. Lookup.
 534                     forAll(cz, zoneI)
 535                     {
 536                         label index = cz[zoneI].whichCell(curCellI);
 537
 538                         if (index != -1)
 539                         {
 540                             zoneCells[zoneI].append(celli);
 541                         }
 542                     }
 543                 }
 544             }
 545
 546             procMesh.cellZones().clearAddressing();
 547             procMesh.cellZones().setSize(zoneCells.size());
 548             forAll(zoneCells, zoneI)
 549             {
 550                 procMesh.cellZones().set
 551                 (
 552                     zoneI,
 553                     cz[zoneI].clone
 554                     (
 555                         zoneCells[zoneI].shrink(),
 556                         zoneI,
 557                         procMesh.cellZones()
 558                     )
 559                 );
 560             }
 561
 562             if (cz.size())
 563             {
 564                 // Force writing on all processors
 565                 procMesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
 566             }
 567         }
 568
 569         // Set the precision of the points data to 10
 570         IOstream::defaultPrecision(10);
 571
 572         procMesh.write();
 573
 574         Info<< endl
 575             << "Processor " << procI << nl
 576             << "    Number of cells = " << procMesh.nCells()
 577             << endl;
 578
 579         label nBoundaryFaces = 0;
 580         label nProcPatches = 0;
 581         label nProcFaces = 0;
 582
 583         forAll (procMesh.boundaryMesh(), patchi)
 584         {
 585             if
 586             (
 587                 procMesh.boundaryMesh()[patchi].type()
 588              == processorPolyPatch::typeName
 589             )
 590             {
 591                 const processorPolyPatch& ppp =
 592                 refCast<const processorPolyPatch>
 593                 (
 594                     procMesh.boundaryMesh()[patchi]
 595                 );
 596
 597                 Info<< "    Number of faces shared with processor "
 598                     << ppp.neighbProcNo() << " = " << ppp.size() << endl;
 599
 600                 nProcPatches++;
 601                 nProcFaces += ppp.size();
 602             }
 603             else
 604             {
 605                 nBoundaryFaces += procMesh.boundaryMesh()[patchi].size();
 606             }
 607         }
 608
 609         Info<< "    Number of processor patches = " << nProcPatches << nl
 610             << "    Number of processor faces = " << nProcFaces << nl
 611             << "    Number of boundary faces = " << nBoundaryFaces << endl;
 612
 613         totProcFaces += nProcFaces;
 614         maxProcPatches = max(maxProcPatches, nProcPatches);
 615         maxProcFaces = max(maxProcFaces, nProcFaces);
 616
 617         // create and write the addressing information
 618         labelIOList pointProcAddressing
 619         (
 620             IOobject
 621             (
 622                 "pointProcAddressing",
 623                 procMesh.facesInstance(),
 624                 procMesh.meshSubDir,
 625                 procMesh,
 626                 IOobject::NO_READ,
 627                 IOobject::NO_WRITE
 628             ),
 629             procPointAddressing_[procI]
 630         );
 631         pointProcAddressing.write();
 632
 633         labelIOList faceProcAddressing
 634         (
 635             IOobject
 636             (
 637                 "faceProcAddressing",
 638                 procMesh.facesInstance(),
 639                 procMesh.meshSubDir,
 640                 procMesh,
 641                 IOobject::NO_READ,
 642                 IOobject::NO_WRITE
 643             ),
 644             procFaceAddressing_[procI]
 645         );
 646         faceProcAddressing.write();
 647
 648         labelIOList cellProcAddressing
 649         (
 650             IOobject
 651             (
 652                 "cellProcAddressing",
 653                 procMesh.facesInstance(),
 654                 procMesh.meshSubDir,
 655                 procMesh,
 656                 IOobject::NO_READ,
 657                 IOobject::NO_WRITE
 658             ),
 659             procCellAddressing_[procI]
 660         );
 661         cellProcAddressing.write();
 662
 663         labelIOList boundaryProcAddressing
 664         (
 665             IOobject
 666             (
 667                 "boundaryProcAddressing",
 668                 procMesh.facesInstance(),
 669                 procMesh.meshSubDir,
 670                 procMesh,
 671                 IOobject::NO_READ,
 672                 IOobject::NO_WRITE
 673             ),
 674             procBoundaryAddressing_[procI]
 675         );
 676         boundaryProcAddressing.write();
 677     }
 678
 679     Info<< nl
 680         << "Number of processor faces = " << totProcFaces/2 << nl
 681         << "Max number of processor patches = " << maxProcPatches << nl
 682         << "Max number of faces between processors = " << maxProcFaces
 683         << endl;
 684
 685     return true;
 686 }
 687
 688
 689 // ************************************************************************* //