applications/utilities/parallelProcessing/decomposePar/decomposeMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 InClass
  25     domainDecomposition
  26
  27 Description
  28     Private member of domainDecomposition.
  29     Decomposes the mesh into bits
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "domainDecomposition.H"
  34 #include "IOstreams.H"
  35 #include "SLPtrList.H"
  36 #include "boolList.H"
  37 #include "cellList.H"
  38 #include "primitiveMesh.H"
  39 #include "cyclicPolyPatch.H"
  40
  41 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  42
  43 void domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
  44 {
  45     // Decide which cell goes to which processor
  46     distributeCells();
  47
  48     // Distribute the cells according to the given processor label
  49
  50     // calculate the addressing information for the original mesh
  51     Info<< "\nCalculating original mesh data" << endl;
  52
  53     // set references to the original mesh
  54     const polyBoundaryMesh& patches = boundaryMesh();
  55
  56     // Access all faces to grab the zones
  57     const faceList& fcs = allFaces();
  58     const labelList& owner = faceOwner();
  59     const labelList& neighbour = faceNeighbour();
  60
  61     // loop through the list of processor labels for the cell and add the
  62     // cell shape to the list of cells for the appropriate processor
  63
  64     Info<< "\nDistributing cells to processors" << endl;
  65
  66     // Memory management
  67     {
  68         List<SLList<label> > procCellList(nProcs_);
  69
  70         forAll (cellToProc_, celli)
  71         {
  72             if (cellToProc_[celli] >= nProcs_)
  73             {
  74                 FatalErrorIn("domainDecomposition::decomposeMesh()")
  75                     << "Impossible processor label " << cellToProc_[celli]
  76                     << "for cell " << celli
  77                     << abort(FatalError);
  78             }
  79             else
  80             {
  81                 procCellList[cellToProc_[celli]].append(celli);
  82             }
  83         }
  84
  85         // Convert linked lists into normal lists
  86         forAll (procCellList, procI)
  87         {
  88             procCellAddressing_[procI] = procCellList[procI];
  89         }
  90     }
  91
  92     Info << "\nDistributing faces to processors" << endl;
  93
  94     // Loop through internal faces and decide which processor they belong to
  95     // First visit all internal faces. If cells at both sides belong to the
  96     // same processor, the face is an internal face. If they are different,
  97     // it belongs to both processors.
  98
  99     // Memory management
 100     {
 101         List<SLList<label> > procFaceList(nProcs_);
 102
 103         forAll (neighbour, facei)
 104         {
 105             if (cellToProc_[owner[facei]] == cellToProc_[neighbour[facei]])
 106             {
 107                 // Face internal to processor
 108                 procFaceList[cellToProc_[owner[facei]]].append(facei);
 109             }
 110         }
 111
 112         // Record number of internal faces on each processor
 113         forAll (procFaceList, procI)
 114         {
 115             nInternalProcFaces_[procI] = procFaceList[procI].size();
 116         }
 117
 118         // Detect inter-processor boundaries
 119
 120         // Neighbour processor for each subdomain
 121         List<SLList<label> > interProcBoundaries(nProcs_);
 122
 123         // Face labels belonging to each inter-processor boundary
 124         List<SLList<SLList<label> > > interProcBFaces(nProcs_);
 125
 126         List<SLList<label> > procPatchIndex(nProcs_);
 127
 128         forAll (neighbour, facei)
 129         {
 130             if (cellToProc_[owner[facei]] != cellToProc_[neighbour[facei]])
 131             {
 132                 // inter - processor patch face found. Go through the list of
 133                 // inside boundaries for the owner processor and try to find
 134                 // this inter-processor patch.
 135
 136                 label ownerProc = cellToProc_[owner[facei]];
 137                 label neighbourProc = cellToProc_[neighbour[facei]];
 138
 139                 SLList<label>::iterator curInterProcBdrsOwnIter =
 140                     interProcBoundaries[ownerProc].begin();
 141
 142                 SLList<SLList<label> >::iterator curInterProcBFacesOwnIter =
 143                     interProcBFaces[ownerProc].begin();
 144
 145                 bool interProcBouFound = false;
 146
 147                 // WARNING: Synchronous SLList iterators
 148
 149                 for
 150                 (
 151                     ;
 152                     curInterProcBdrsOwnIter
 153                  != interProcBoundaries[ownerProc].end()
 154                  && curInterProcBFacesOwnIter
 155                  != interProcBFaces[ownerProc].end();
 156                     ++curInterProcBdrsOwnIter, ++curInterProcBFacesOwnIter
 157                 )
 158                 {
 159                     if (curInterProcBdrsOwnIter() == neighbourProc)
 160                     {
 161                         // the inter - processor boundary exists. Add the face
 162                         interProcBouFound = true;
 163
 164                         curInterProcBFacesOwnIter().append(facei);
 165
 166                         SLList<label>::iterator curInterProcBdrsNeiIter =
 167                             interProcBoundaries[neighbourProc].begin();
 168
 169                         SLList<SLList<label> >::iterator
 170                             curInterProcBFacesNeiIter =
 171                             interProcBFaces[neighbourProc].begin();
 172
 173                         bool neighbourFound = false;
 174
 175                         // WARNING: Synchronous SLList iterators
 176
 177                         for
 178                         (
 179                             ;
 180                             curInterProcBdrsNeiIter !=
 181                             interProcBoundaries[neighbourProc].end()
 182                          && curInterProcBFacesNeiIter !=
 183                             interProcBFaces[neighbourProc].end();
 184                             ++curInterProcBdrsNeiIter,
 185                             ++curInterProcBFacesNeiIter
 186                         )
 187                         {
 188                             if (curInterProcBdrsNeiIter() == ownerProc)
 189                             {
 190                                 // boundary found. Add the face
 191                                 neighbourFound = true;
 192
 193                                 curInterProcBFacesNeiIter().append(facei);
 194                             }
 195
 196                             if (neighbourFound) break;
 197                         }
 198
 199                         if (interProcBouFound && !neighbourFound)
 200                         {
 201                             FatalErrorIn
 202                             (
 203                                 "domainDecomposition::decomposeMesh()"
 204                             )   << "Inconsistency in inter - "
 205                                 << "processor boundary lists for processors "
 206                                 << ownerProc << " and " << neighbourProc
 207                                 << abort(FatalError);
 208                         }
 209                     }
 210
 211                     if (interProcBouFound) break;
 212                 }
 213
 214                 if (!interProcBouFound)
 215                 {
 216                     // inter - processor boundaries do not exist and need to
 217                     // be created
 218
 219                     // set the new addressing information
 220
 221                     // owner
 222                     interProcBoundaries[ownerProc].append(neighbourProc);
 223                     interProcBFaces[ownerProc].append(SLList<label>(facei));
 224
 225                     // neighbour
 226                     interProcBoundaries[neighbourProc].append(ownerProc);
 227                     interProcBFaces[neighbourProc].append
 228                     (
 229                         SLList<label>(facei)
 230                     );
 231                 }
 232             }
 233         }
 234
 235         // Loop through patches. For cyclic boundaries detect inter-processor
 236         // faces; for all other, add faces to the face list and remember start
 237         // and size of all patches.
 238
 239         // for all processors, set the size of start index and patch size
 240         // lists to the number of patches in the mesh
 241         forAll (procPatchSize_, procI)
 242         {
 243             procPatchSize_[procI].setSize(patches.size());
 244             procPatchStartIndex_[procI].setSize(patches.size());
 245         }
 246
 247         forAll (patches, patchi)
 248         {
 249             // Reset size and start index for all processors
 250             forAll (procPatchSize_, procI)
 251             {
 252                 procPatchSize_[procI][patchi] = 0;
 253                 procPatchStartIndex_[procI][patchi] =
 254                     procFaceList[procI].size();
 255             }
 256
 257             const label patchStart = patches[patchi].start();
 258
 259             if (!isA<cyclicPolyPatch>(patches[patchi]))
 260             {
 261                 // Normal patch. Add faces to processor where the cell
 262                 // next to the face lives
 263
 264                 const unallocLabelList& patchFaceCells =
 265                     patches[patchi].faceCells();
 266
 267                 forAll (patchFaceCells, facei)
 268                 {
 269                     const label curProc = cellToProc_[patchFaceCells[facei]];
 270
 271                     // add the face
 272                     procFaceList[curProc].append(patchStart + facei);
 273
 274                     // increment the number of faces for this patch
 275                     procPatchSize_[curProc][patchi]++;
 276                 }
 277             }
 278             else
 279             {
 280                 // Cyclic patch special treatment
 281
 282                 const polyPatch& cPatch = patches[patchi];
 283
 284                 const label cycOffset = cPatch.size()/2;
 285
 286                 // Set reference to faceCells for both patches
 287                 const labelList::subList firstFaceCells
 288                 (
 289                     cPatch.faceCells(),
 290                     cycOffset
 291                 );
 292
 293                 const labelList::subList secondFaceCells
 294                 (
 295                     cPatch.faceCells(),
 296                     cycOffset,
 297                     cycOffset
 298                 );
 299
 300                 forAll (firstFaceCells, facei)
 301                 {
 302                     if
 303                     (
 304                         cellToProc_[firstFaceCells[facei]]
 305                      != cellToProc_[secondFaceCells[facei]]
 306                     )
 307                     {
 308                         // This face becomes an inter-processor boundary face
 309                         // inter - processor patch face found. Go through
 310                         // the list of inside boundaries for the owner
 311                         // processor and try to find this inter-processor
 312                         // patch.
 313
 314                         cyclicParallel_ = true;
 315
 316                         label ownerProc = cellToProc_[firstFaceCells[facei]];
 317                         label neighbourProc =
 318                             cellToProc_[secondFaceCells[facei]];
 319
 320                         SLList<label>::iterator curInterProcBdrsOwnIter =
 321                             interProcBoundaries[ownerProc].begin();
 322
 323                         SLList<SLList<label> >::iterator
 324                             curInterProcBFacesOwnIter =
 325                             interProcBFaces[ownerProc].begin();
 326
 327                         bool interProcBouFound = false;
 328
 329                         // WARNING: Synchronous SLList iterators
 330
 331                         for
 332                         (
 333                             ;
 334                             curInterProcBdrsOwnIter !=
 335                             interProcBoundaries[ownerProc].end()
 336                          && curInterProcBFacesOwnIter !=
 337                             interProcBFaces[ownerProc].end();
 338                             ++curInterProcBdrsOwnIter,
 339                             ++curInterProcBFacesOwnIter
 340                         )
 341                         {
 342                             if (curInterProcBdrsOwnIter() == neighbourProc)
 343                             {
 344                                 // the inter - processor boundary exists.
 345                                 // Add the face
 346                                 interProcBouFound = true;
 347
 348                                 curInterProcBFacesOwnIter().append
 349                                     (patchStart + facei);
 350
 351                                 SLList<label>::iterator curInterProcBdrsNeiIter
 352                                    = interProcBoundaries[neighbourProc].begin();
 353
 354                                 SLList<SLList<label> >::iterator
 355                                     curInterProcBFacesNeiIter =
 356                                     interProcBFaces[neighbourProc].begin();
 357
 358                                 bool neighbourFound = false;
 359
 360                                 // WARNING: Synchronous SLList iterators
 361
 362                                 for
 363                                 (
 364                                     ;
 365                                     curInterProcBdrsNeiIter
 366                                    != interProcBoundaries[neighbourProc].end()
 367                                  && curInterProcBFacesNeiIter
 368                                    != interProcBFaces[neighbourProc].end();
 369                                     ++curInterProcBdrsNeiIter,
 370                                     ++curInterProcBFacesNeiIter
 371                                 )
 372                                 {
 373                                     if (curInterProcBdrsNeiIter() == ownerProc)
 374                                     {
 375                                         // boundary found. Add the face
 376                                         neighbourFound = true;
 377
 378                                         curInterProcBFacesNeiIter()
 379                                             .append
 380                                             (
 381                                                 patchStart
 382                                               + cycOffset
 383                                               + facei
 384                                             );
 385                                     }
 386
 387                                     if (neighbourFound) break;
 388                                 }
 389
 390                                 if (interProcBouFound && !neighbourFound)
 391                                 {
 392                                     FatalErrorIn
 393                                     (
 394                                         "domainDecomposition::decomposeMesh()"
 395                                     )   << "Inconsistency in inter-processor "
 396                                         << "boundary lists for processors "
 397                                         << ownerProc << " and "
 398                                         << neighbourProc
 399                                         << " in cyclic boundary matching"
 400                                         << abort(FatalError);
 401                                 }
 402                             }
 403
 404                             if (interProcBouFound) break;
 405                         }
 406
 407                         if (!interProcBouFound)
 408                         {
 409                             // inter - processor boundaries do not exist
 410                             // and need to be created
 411
 412                             // set the new addressing information
 413
 414                             // owner
 415                             interProcBoundaries[ownerProc]
 416                                 .append(neighbourProc);
 417                             interProcBFaces[ownerProc]
 418                                 .append(SLList<label>(patchStart + facei));
 419
 420                             // neighbour
 421                             interProcBoundaries[neighbourProc]
 422                                 .append(ownerProc);
 423                             interProcBFaces[neighbourProc]
 424                                 .append
 425                                 (
 426                                     SLList<label>
 427                                     (
 428                                         patchStart
 429                                       + cycOffset
 430                                       + facei
 431                                     )
 432                                 );
 433                         }
 434                     }
 435                     else
 436                     {
 437                         // This cyclic face remains on the processor
 438                         label ownerProc = cellToProc_[firstFaceCells[facei]];
 439
 440                         // add the face
 441                         procFaceList[ownerProc].append(patchStart + facei);
 442
 443                         // increment the number of faces for this patch
 444                         procPatchSize_[ownerProc][patchi]++;
 445
 446                         // Note: I cannot add the other side of the cyclic
 447                         // boundary here because this would violate the order.
 448                         // They will be added in a separate loop below
 449                         // HJ, 15/Jan/2001
 450                     }
 451                 }
 452
 453                 // Ordering in cyclic boundaries is important.
 454                 // Add the other half of cyclic faces for cyclic boundaries
 455                 // that remain on the processor
 456                 forAll (secondFaceCells, facei)
 457                 {
 458                     if
 459                     (
 460                         cellToProc_[firstFaceCells[facei]]
 461                      == cellToProc_[secondFaceCells[facei]]
 462                     )
 463                     {
 464                         // This cyclic face remains on the processor
 465                         label ownerProc = cellToProc_[firstFaceCells[facei]];
 466
 467                         // add the second face
 468                         procFaceList[ownerProc].append
 469                             (patchStart + cycOffset + facei);
 470
 471                         // increment the number of faces for this patch
 472                         procPatchSize_[ownerProc][patchi]++;
 473                     }
 474                 }
 475             }
 476         }
 477
 478         // Face zone treatment.  HJ, 27/Mar/2009
 479         // Face zones identified as global will be present on all CPUs
 480         List<SLList<label> > procZoneFaceList(nProcs_);
 481
 482         if (decompositionDict_.found("globalFaceZones"))
 483         {
 484             wordList fzNames(decompositionDict_.lookup("globalFaceZones"));
 485
 486             const faceZoneMesh& fz = faceZones();
 487
 488             forAll (fzNames, nameI)
 489             {
 490                 const label zoneID = fz.findZoneID(fzNames[nameI]);
 491
 492                 if (zoneID == -1)
 493                 {
 494                     FatalErrorIn("domainDecomposition::decomposeMesh()")
 495                         << "Unknown global face zone " << fzNames[nameI]
 496                         << nl << "Valid face zones are" << fz.names()
 497                         << exit(FatalError);
 498                 }
 499
 500                 Info<< "Preserving global face zone " << fzNames[nameI]
 501                     << endl;
 502
 503                 const faceZone& curFz = fz[zoneID];
 504
 505                 // Go through all the faces in the zone.  If the owner of the
 506                 // face equals to current processor, it has already been added;
 507                 // otherwise, add the face to all processor face lists
 508                 forAll (curFz, faceI)
 509                 {
 510                     const label curFaceID = curFz[faceI];
 511
 512                     if (curFaceID < owner.size())
 513                     {
 514                         // This is an active mesh face, and it already belongs
 515                         // to one CPU.  Find out which and add it to the others
 516
 517                         const label curProc = cellToProc_[owner[curFaceID]];
 518
 519                         forAll (procZoneFaceList, procI)
 520                         {
 521                             if (procI != curProc)
 522                             {
 523                                 procZoneFaceList[procI].append(curFaceID);
 524                             }
 525                         }
 526                     }
 527                     else
 528                     {
 529                         // This is a stand-alone face, add it to all processors
 530                         forAll (procFaceList, procI)
 531                         {
 532                             procZoneFaceList[procI].append(curFaceID);
 533                         }
 534                     }
 535                 }
 536             }
 537         }
 538
 539
 540         // Convert linked lists into normal lists
 541         // Add inter-processor boundaries and remember start indices
 542
 543         forAll (procFaceList, procI)
 544         {
 545             // Get internal and regular boundary processor faces
 546             const SLList<label>& curProcFaces = procFaceList[procI];
 547
 548             // Get reference to processor face addressing
 549             labelList& curProcFaceAddressing = procFaceAddressing_[procI];
 550
 551             labelList& curProcNeighbourProcessors =
 552                 procNeighbourProcessors_[procI];
 553
 554             labelList& curProcProcessorPatchSize =
 555                 procProcessorPatchSize_[procI];
 556
 557             labelList& curProcProcessorPatchStartIndex =
 558                 procProcessorPatchStartIndex_[procI];
 559
 560             // calculate the size
 561             label nFacesOnProcessor = curProcFaces.size();
 562
 563             for
 564             (
 565                 SLList<SLList<label> >::const_iterator curInterProcBFacesIter =
 566                     interProcBFaces[procI].begin();
 567                 curInterProcBFacesIter != interProcBFaces[procI].end();
 568                 ++curInterProcBFacesIter
 569             )
 570             {
 571                 nFacesOnProcessor += curInterProcBFacesIter().size();
 572             }
 573
 574             // Add stand-alone global zone faces
 575             nFacesOnProcessor += procZoneFaceList[procI].size();
 576
 577             curProcFaceAddressing.setSize(nFacesOnProcessor);
 578
 579             // Fill in the list. Calculate turning index.
 580             // Turning index will be -1 only for some faces on processor
 581             // boundaries, i.e. the ones where the current processor ID
 582             // is in the cell which is a face neighbour.
 583             // Turning index is stored as the sign of the face addressing list
 584
 585             label nFaces = 0;
 586
 587             // Add internal and boundary faces
 588             // Remember to increment the index by one such that the
 589             // turning index works properly.  HJ, 5/Dec/2001
 590             for
 591             (
 592                 SLList<label>::const_iterator curProcFacesIter =
 593                     curProcFaces.begin();
 594                 curProcFacesIter != curProcFaces.end();
 595                 ++curProcFacesIter
 596             )
 597             {
 598                 curProcFaceAddressing[nFaces] = curProcFacesIter() + 1;
 599                 nFaces++;
 600             }
 601
 602             // Add inter-processor boundary faces. At the beginning of each
 603             // patch, grab the patch start index and size
 604
 605             curProcNeighbourProcessors.setSize
 606             (
 607                 interProcBoundaries[procI].size()
 608             );
 609
 610             curProcProcessorPatchSize.setSize
 611             (
 612                 interProcBoundaries[procI].size()
 613             );
 614
 615             curProcProcessorPatchStartIndex.setSize
 616             (
 617                 interProcBoundaries[procI].size()
 618             );
 619
 620             label nProcPatches = 0;
 621
 622             SLList<label>::iterator curInterProcBdrsIter =
 623                 interProcBoundaries[procI].begin();
 624
 625             SLList<SLList<label> >::iterator curInterProcBFacesIter =
 626                 interProcBFaces[procI].begin();
 627
 628             for
 629             (
 630                 ;
 631                 curInterProcBdrsIter != interProcBoundaries[procI].end()
 632              && curInterProcBFacesIter != interProcBFaces[procI].end();
 633                 ++curInterProcBdrsIter, ++curInterProcBFacesIter
 634             )
 635             {
 636                 curProcNeighbourProcessors[nProcPatches] =
 637                     curInterProcBdrsIter();
 638
 639                 // Get start index for processor patch
 640                 curProcProcessorPatchStartIndex[nProcPatches] = nFaces;
 641
 642                 label& curSize =
 643                     curProcProcessorPatchSize[nProcPatches];
 644
 645                 curSize = 0;
 646
 647                 // Add faces for this processor boundary
 648
 649                 for
 650                 (
 651                     SLList<label>::iterator curFacesIter =
 652                         curInterProcBFacesIter().begin();
 653                     curFacesIter != curInterProcBFacesIter().end();
 654                     ++curFacesIter
 655                 )
 656                 {
 657                     // Add the face
 658
 659                     // Remember to increment the index by one such that the
 660                     // turning index works properly.  HJ, 5/Dec/2001
 661                     if (cellToProc_[owner[curFacesIter()]] == procI)
 662                     {
 663                         curProcFaceAddressing[nFaces] = curFacesIter() + 1;
 664                     }
 665                     else
 666                     {
 667                         // Turning face
 668                         curProcFaceAddressing[nFaces] = -(curFacesIter() + 1);
 669                     }
 670
 671                     // increment the size
 672                     curSize++;
 673
 674                     nFaces++;
 675                 }
 676
 677                 nProcPatches++;
 678             }
 679
 680             // Record number of live faces
 681             nLiveProcFaces_[procI] = nFaces;
 682
 683             // Add stand-alone face zone faces
 684             const SLList<label>& curProcZoneFaces = procZoneFaceList[procI];
 685
 686             for
 687             (
 688                 SLList<label>::const_iterator curProcZoneFacesIter =
 689                     curProcZoneFaces.begin();
 690                 curProcZoneFacesIter != curProcZoneFaces.end();
 691                 ++curProcZoneFacesIter
 692             )
 693             {
 694                 curProcFaceAddressing[nFaces] = curProcZoneFacesIter() + 1;
 695                 nFaces++;
 696             }
 697         } // End for all processors
 698     } // End of memory management
 699
 700     Info << "\nCalculating processor boundary addressing" << endl;
 701     // For every patch of processor boundary, find the index of the original
 702     // patch. Mis-alignment is caused by the fact that patches with zero size
 703     // are omitted. For processor patches, set index to -1.
 704     // At the same time, filter the procPatchSize_ and procPatchStartIndex_
 705     // lists to exclude zero-size patches
 706     forAll (procPatchSize_, procI)
 707     {
 708         // Make a local copy of old lists
 709         const labelList oldPatchSizes = procPatchSize_[procI];
 710
 711         const labelList oldPatchStarts = procPatchStartIndex_[procI];
 712
 713         labelList& curPatchSizes = procPatchSize_[procI];
 714
 715         labelList& curPatchStarts = procPatchStartIndex_[procI];
 716
 717         const labelList& curProcessorPatchSizes =
 718             procProcessorPatchSize_[procI];
 719
 720         labelList& curBoundaryAddressing = procBoundaryAddressing_[procI];
 721
 722         curBoundaryAddressing.setSize
 723         (
 724             oldPatchSizes.size()
 725           + curProcessorPatchSizes.size()
 726         );
 727
 728         label nPatches = 0;
 729
 730         forAll (oldPatchSizes, patchi)
 731         {
 732             if (!filterEmptyPatches || oldPatchSizes[patchi] > 0)
 733             {
 734                 curBoundaryAddressing[nPatches] = patchi;
 735
 736                 curPatchSizes[nPatches] = oldPatchSizes[patchi];
 737
 738                 curPatchStarts[nPatches] = oldPatchStarts[patchi];
 739
 740                 nPatches++;
 741             }
 742         }
 743
 744         // reset to the size of live patches
 745         curPatchSizes.setSize(nPatches);
 746         curPatchStarts.setSize(nPatches);
 747
 748         forAll (curProcessorPatchSizes, procPatchI)
 749         {
 750             curBoundaryAddressing[nPatches] = -1;
 751
 752             nPatches++;
 753         }
 754
 755         curBoundaryAddressing.setSize(nPatches);
 756     }
 757
 758     Info << "\nDistributing points to processors" << endl;
 759     // For every processor, loop through the list of faces for the processor.
 760     // For every face, loop through the list of points and mark the point as
 761     // used for the processor. Collect the list of used points for the
 762     // processor.
 763
 764     // Record number of live points on each processor
 765     labelList nLivePoints(nProcs_, 0);
 766
 767     forAll (procPointAddressing_, procI)
 768     {
 769         // Dimension list to all points in the mesh.  HJ, 27/Mar/2009
 770         boolList pointLabels(allPoints().size(), false);
 771
 772         // Get reference to list of used faces
 773         const labelList& procFaceLabels = procFaceAddressing_[procI];
 774
 775         // Collect the used points
 776         labelList& procPointLabels = procPointAddressing_[procI];
 777
 778         procPointLabels.setSize(pointLabels.size());
 779
 780         // Two-pass algorithm:
 781         // First loop through live faces and record them in the points list
 782         // Second, visit all inactive zone faces and record the points
 783
 784         label nUsedPoints = 0;
 785
 786         // First pass: live faces
 787
 788         for (label faceI = 0; faceI < nLiveProcFaces_[procI]; faceI++)
 789         {
 790             // Because of the turning index, some labels may be negative
 791             const labelList& facePoints = fcs[mag(procFaceLabels[faceI]) - 1];
 792
 793             forAll (facePoints, pointI)
 794             {
 795                 // Mark the point as used
 796                 pointLabels[facePoints[pointI]] = true;
 797             }
 798         }
 799
 800         forAll (pointLabels, pointI)
 801         {
 802             if (pointLabels[pointI])
 803             {
 804                 procPointLabels[nUsedPoints] = pointI;
 805
 806                 nUsedPoints++;
 807             }
 808         }
 809
 810         // Record number of live points
 811         nLivePoints[procI] = nUsedPoints;
 812
 813         // Second pass: zone faces
 814
 815         // Reset point usage list
 816         boolList pointLabelsSecondPass(allPoints().size(), false);
 817
 818         for
 819         (
 820             label faceI = nLiveProcFaces_[procI];
 821             faceI < procFaceLabels.size();
 822             faceI++
 823         )
 824         {
 825             // Because of the turning index, some labels may be negative
 826             const labelList& facePoints = fcs[mag(procFaceLabels[faceI]) - 1];
 827
 828             forAll (facePoints, pointI)
 829             {
 830                 // Mark the point as used
 831                 if (!pointLabels[facePoints[pointI]])
 832                 {
 833                     pointLabelsSecondPass[facePoints[pointI]] = true;
 834                 }
 835             }
 836         }
 837
 838         forAll (pointLabelsSecondPass, pointI)
 839         {
 840             if (pointLabelsSecondPass[pointI])
 841             {
 842                 procPointLabels[nUsedPoints] = pointI;
 843
 844                 nUsedPoints++;
 845             }
 846         }
 847
 848         // Reset the size of used points
 849         procPointLabels.setSize(nUsedPoints);
 850     }
 851
 852     // Gather data about globally shared points
 853
 854     // Memory management
 855     {
 856         // Dimension list to all points in the mesh.  HJ, 27/Mar/2009
 857         labelList pointsUsage(allPoints().size(), 0);
 858
 859         // Globally shared points are the ones used by more than 2 processors
 860         // Size the list approximately and gather the points
 861         labelHashSet gSharedPoints
 862         (
 863             min(100, nPoints()/1000)
 864         );
 865
 866         // Loop through all the processors and mark up points used by
 867         // processor boundaries.  When a point is used twice, it is a
 868         // globally shared point
 869
 870         for (label procI = 0; procI < nProcs_; procI++)
 871         {
 872             // Get list of face labels
 873             const labelList& curFaceLabels = procFaceAddressing_[procI];
 874
 875             // Get start of processor faces
 876             const labelList& curProcessorPatchStarts =
 877                 procProcessorPatchStartIndex_[procI];
 878
 879             const labelList& curProcessorPatchSizes =
 880                 procProcessorPatchSize_[procI];
 881
 882             // Reset the lookup list
 883             pointsUsage = 0;
 884
 885             forAll (curProcessorPatchStarts, patchi)
 886             {
 887                 const label curStart = curProcessorPatchStarts[patchi];
 888                 const label curEnd = curStart + curProcessorPatchSizes[patchi];
 889
 890                 for
 891                 (
 892                     label faceI = curStart;
 893                     faceI < curEnd;
 894                     faceI++
 895                 )
 896                 {
 897                     // Mark the original face as used
 898                     // Remember to decrement the index by one (turning index)
 899                     // HJ, 5/Dec/2001
 900                     const label curF = mag(curFaceLabels[faceI]) - 1;
 901
 902                     const face& f = fcs[curF];
 903
 904                     forAll (f, pointI)
 905                     {
 906                         if (pointsUsage[f[pointI]] == 0)
 907                         {
 908                             // Point not previously used
 909                             pointsUsage[f[pointI]] = patchi + 1;
 910                         }
 911                         else if (pointsUsage[f[pointI]] != patchi + 1)
 912                         {
 913                             // Point used by some other patch = global point!
 914                             gSharedPoints.insert(f[pointI]);
 915                         }
 916                     }
 917                 }
 918             }
 919         }
 920
 921         // Grab the result from the hash list
 922         globallySharedPoints_ = gSharedPoints.toc();
 923         sort(globallySharedPoints_);
 924     }
 925 }