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