src/meshTools/regionSplit/regionSplit.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 \*---------------------------------------------------------------------------*/
  25
  26 #include "regionSplit.H"
  27 #include "cyclicPolyPatch.H"
  28 #include "processorPolyPatch.H"
  29 #include "globalIndex.H"
  30 #include "syncTools.H"
  31
  32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36
  37 defineTypeNameAndDebug(regionSplit, 0);
  38
  39 }
  40
  41 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  42
  43 // Handle (non-processor) coupled faces.
  44 void Foam::regionSplit::transferCoupledFaceRegion
  45 (
  46     const label faceI,
  47     const label otherFaceI,
  48
  49     labelList& faceRegion,
  50     dynamicLabelList& newChangedFaces
  51 ) const
  52 {
  53     if (faceRegion[faceI] >= 0)
  54     {
  55         if (faceRegion[otherFaceI] == -1)
  56         {
  57             faceRegion[otherFaceI] = faceRegion[faceI];
  58             newChangedFaces.append(otherFaceI);
  59         }
  60         else if (faceRegion[otherFaceI] == -2)
  61         {
  62             // otherFaceI blocked but faceI is not. Is illegal for coupled
  63             // faces, not for explicit connections.
  64         }
  65         else if (faceRegion[otherFaceI] != faceRegion[faceI])
  66         {
  67             FatalErrorIn
  68             (
  69                   "regionSplit::transferCoupledFaceRegion"
  70                   "(const label, const label, labelList&, labelList&) const"
  71               )   << "Problem : coupled face " << faceI
  72                   << " on patch " << mesh_.boundaryMesh().whichPatch(faceI)
  73                   << " has region " << faceRegion[faceI]
  74                   << " but coupled face " << otherFaceI
  75                   << " has region " << faceRegion[otherFaceI]
  76                   << endl
  77                   << "Is your blocked faces specification"
  78                   << " synchronized across coupled boundaries?"
  79                   << abort(FatalError);
  80         }
  81     }
  82     else if (faceRegion[faceI] == -1)
  83     {
  84         if (faceRegion[otherFaceI] >= 0)
  85         {
  86             faceRegion[faceI] = faceRegion[otherFaceI];
  87             newChangedFaces.append(faceI);
  88         }
  89         else if (faceRegion[otherFaceI] == -2)
  90         {
  91             // otherFaceI blocked but faceI is not. Is illegal for coupled
  92             // faces, not for explicit connections.
  93         }
  94     }
  95 }
  96
  97
  98 void Foam::regionSplit::fillSeedMask
  99 (
 100     const List<labelPair>& explicitConnections,
 101     labelList& cellRegion,
 102     labelList& faceRegion,
 103     const label seedCellID,
 104     const label markValue
 105 ) const
 106 {
 107     // Do seed cell
 108     cellRegion[seedCellID] = markValue;
 109
 110
 111     // Collect faces on seed cell
 112     const cell& cFaces = mesh_.cells()[seedCellID];
 113
 114     label nFaces = 0;
 115
 116     labelList changedFaces(cFaces.size());
 117
 118     forAll(cFaces, i)
 119     {
 120         label faceI = cFaces[i];
 121
 122         if (faceRegion[faceI] == -1)
 123         {
 124             faceRegion[faceI] = markValue;
 125             changedFaces[nFaces++] = faceI;
 126         }
 127     }
 128     changedFaces.setSize(nFaces);
 129
 130
 131     // Loop over changed faces. MeshWave in small.
 132
 133     while (changedFaces.size())
 134     {
 135         //if (debug)
 136         //{
 137         //    Pout<< "regionSplit::fillSeedMask : changedFaces:"
 138         //        << changedFaces.size() << endl;
 139         //}
 140
 141         dynamicLabelList changedCells(changedFaces.size());
 142
 143         forAll(changedFaces, i)
 144         {
 145             label faceI = changedFaces[i];
 146
 147             label own = mesh_.faceOwner()[faceI];
 148
 149             if (cellRegion[own] == -1)
 150             {
 151                 cellRegion[own] = markValue;
 152                 changedCells.append(own);
 153             }
 154
 155             if (mesh_.isInternalFace(faceI))
 156             {
 157                 label nei = mesh_.faceNeighbour()[faceI];
 158
 159                 if (cellRegion[nei] == -1)
 160                 {
 161                     cellRegion[nei] = markValue;
 162                     changedCells.append(nei);
 163                 }
 164             }
 165         }
 166
 167
 168         //if (debug)
 169         //{
 170         //    Pout<< "regionSplit::fillSeedMask : changedCells:"
 171         //        << changedCells.size() << endl;
 172         //}
 173
 174         // Loop over changedCells and collect faces
 175         dynamicLabelList newChangedFaces(changedCells.size());
 176
 177         forAll(changedCells, i)
 178         {
 179             label cellI = changedCells[i];
 180
 181             const cell& cFaces = mesh_.cells()[cellI];
 182
 183             forAll(cFaces, cFaceI)
 184             {
 185                 label faceI = cFaces[cFaceI];
 186
 187                 if (faceRegion[faceI] == -1)
 188                 {
 189                     faceRegion[faceI] = markValue;
 190                     newChangedFaces.append(faceI);
 191                 }
 192             }
 193         }
 194
 195
 196         //if (debug)
 197         //{
 198         //    Pout<< "regionSplit::fillSeedMask : changedFaces before sync:"
 199         //        << changedFaces.size() << endl;
 200         //}
 201
 202
 203         // Check for changes to any locally coupled face.
 204         // Global connections are done later.
 205
 206         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 207
 208         forAll(patches, patchI)
 209         {
 210             const polyPatch& pp = patches[patchI];
 211
 212             if (isA<cyclicPolyPatch>(pp))
 213             {
 214                 label faceI = pp.start();
 215
 216                 label halfSz = pp.size()/2;
 217
 218                 for (label i = 0; i < halfSz; i++)
 219                 {
 220                     label otherFaceI = refCast<const cyclicPolyPatch>(pp)
 221                         .transformGlobalFace(faceI);
 222
 223                     transferCoupledFaceRegion
 224                     (
 225                         faceI,
 226                         otherFaceI,
 227                         faceRegion,
 228                         newChangedFaces
 229                     );
 230
 231                     faceI++;
 232                 }
 233             }
 234         }
 235         forAll(explicitConnections, i)
 236         {
 237             transferCoupledFaceRegion
 238             (
 239                 explicitConnections[i][0],
 240                 explicitConnections[i][1],
 241                 faceRegion,
 242                 newChangedFaces
 243             );
 244         }
 245
 246         //if (debug)
 247         //{
 248         //    Pout<< "regionSplit::fillSeedMask : changedFaces after sync:"
 249         //        << newChangedFaces.size() << endl;
 250         //}
 251
 252         changedFaces.transfer(newChangedFaces);
 253     }
 254 }
 255
 256
 257 Foam::label Foam::regionSplit::calcRegionSplit
 258 (
 259     const boolList& blockedFace,
 260     const List<labelPair>& explicitConnections,
 261
 262     labelList& cellRegion
 263 ) const
 264 {
 265     if (debug)
 266     {
 267         if (blockedFace.size())
 268         {
 269             // Check that blockedFace is synced.
 270             boolList syncBlockedFace(blockedFace);
 271             syncTools::swapFaceList(mesh_, syncBlockedFace, false);
 272
 273             forAll(syncBlockedFace, faceI)
 274             {
 275                 if (syncBlockedFace[faceI] != blockedFace[faceI])
 276                 {
 277                     FatalErrorIn
 278                     (
 279                         "regionSplit::calcRegionSplit(..)"
 280                     )   << "Face " << faceI << " not synchronised. My value:"
 281                         << blockedFace[faceI] << "  coupled value:"
 282                         << syncBlockedFace[faceI]
 283                         << abort(FatalError);
 284                 }
 285             }
 286         }
 287     }
 288
 289     // Region per face.
 290     // -1 unassigned
 291     // -2 blocked
 292     labelList faceRegion(mesh_.nFaces(), -1);
 293
 294     if (blockedFace.size())
 295     {
 296         forAll(blockedFace, faceI)
 297         {
 298             if (blockedFace[faceI])
 299             {
 300                 faceRegion[faceI] = -2;
 301             }
 302         }
 303     }
 304
 305
 306     // Assign local regions
 307     // ~~~~~~~~~~~~~~~~~~~~
 308
 309     // Start with region 0
 310     label nRegions = 0;
 311
 312     label unsetCellI = 0;
 313
 314     do
 315     {
 316         // Find first unset cell
 317
 318         for (; unsetCellI < mesh_.nCells(); unsetCellI++)
 319         {
 320             if (cellRegion[unsetCellI] == -1)
 321             {
 322                 break;
 323             }
 324         }
 325
 326         if (unsetCellI >= mesh_.nCells())
 327         {
 328             break;
 329         }
 330
 331         fillSeedMask
 332         (
 333             explicitConnections,
 334             cellRegion,
 335             faceRegion,
 336             unsetCellI,
 337             nRegions
 338         );
 339
 340         // Current unsetCell has now been handled. Go to next region.
 341         nRegions++;
 342         unsetCellI++;
 343     }
 344     while(true);
 345
 346
 347     if (debug)
 348     {
 349         forAll(cellRegion, cellI)
 350         {
 351             if (cellRegion[cellI] < 0)
 352             {
 353                 FatalErrorIn("regionSplit::calcRegionSplit(..)")
 354                     << "cell:" << cellI << " region:" << cellRegion[cellI]
 355                     << abort(FatalError);
 356             }
 357         }
 358
 359         forAll(faceRegion, faceI)
 360         {
 361             if (faceRegion[faceI] == -1)
 362             {
 363                 FatalErrorIn("regionSplit::calcRegionSplit(..)")
 364                     << "face:" << faceI << " region:" << faceRegion[faceI]
 365                     << abort(FatalError);
 366             }
 367         }
 368     }
 369
 370
 371
 372     // Assign global regions
 373     // ~~~~~~~~~~~~~~~~~~~~~
 374     // Offset local regions to create unique global regions.
 375
 376     globalIndex globalRegions(nRegions);
 377
 378
 379     // Merge global regions
 380     // ~~~~~~~~~~~~~~~~~~~~
 381     // Regions across non-blocked proc patches get merged.
 382     // This will set merged global regions to be the min of both.
 383     // (this will create gaps in the global region list so they will get
 384     // merged later on)
 385
 386     // Map from global to merged global
 387     labelList mergedGlobal(identity(globalRegions.size()));
 388
 389
 390     // See if any regions get merged. Only nessecary for parallel
 391     while (Pstream::parRun())
 392     {
 393         if (debug)
 394         {
 395             Pout<< nl << "-- Starting Iteration --" << endl;
 396         }
 397
 398         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 399
 400         // Send global regions across (or -2 if blocked face)
 401         forAll(patches, patchI)
 402         {
 403             const polyPatch& pp = patches[patchI];
 404
 405             if (isA<processorPolyPatch>(pp))
 406             {
 407                 labelList myGlobalRegions(pp.size());
 408
 409                 label faceI = pp.start();
 410
 411                 forAll(pp, i)
 412                 {
 413                     if (faceRegion[faceI] < 0)
 414                     {
 415                         myGlobalRegions[i] = faceRegion[faceI];
 416                     }
 417                     else
 418                     {
 419                         myGlobalRegions[i] = mergedGlobal
 420                         [globalRegions.toGlobal(faceRegion[faceI])];
 421                     }
 422
 423                     faceI++;
 424                 }
 425
 426                 OPstream toProcNbr
 427                 (
 428                     Pstream::blocking,
 429                     refCast<const processorPolyPatch>(pp).neighbProcNo()
 430                 );
 431
 432                 toProcNbr << myGlobalRegions;
 433             }
 434         }
 435
 436
 437         // Receive global regions
 438
 439         label nMerged = 0;
 440
 441         forAll(patches, patchI)
 442         {
 443             const polyPatch& pp = patches[patchI];
 444
 445             if (isA<processorPolyPatch>(pp))
 446             {
 447                 const processorPolyPatch& procPp =
 448                     refCast<const processorPolyPatch>(pp);
 449
 450                 IPstream fromProcNbr(Pstream::blocking, procPp.neighbProcNo());
 451
 452                 labelList nbrRegions(fromProcNbr);
 453
 454
 455                 // Compare with my regions to see which get merged.
 456
 457                 label faceI = pp.start();
 458
 459                 forAll(pp, i)
 460                 {
 461                     if
 462                     (
 463                         faceRegion[faceI] < 0
 464                      || nbrRegions[i] < 0
 465                     )
 466                     {
 467                         if (faceRegion[faceI] != nbrRegions[i])
 468                         {
 469                             FatalErrorIn("regionSplit::calcRegionSplit(..)")
 470                                 << "On patch:" << pp.name()
 471                                 << " face:" << faceI
 472                                 << " my local region:" << faceRegion[faceI]
 473                                 << " neighbouring region:"
 474                                 << nbrRegions[i] << nl
 475                                 << "Maybe your blockedFaces are not"
 476                                 << " synchronized across coupled faces?"
 477                                 << abort(FatalError);
 478                         }
 479                     }
 480                     else
 481                     {
 482                         label uncompactGlobal =
 483                             globalRegions.toGlobal(faceRegion[faceI]);
 484
 485                         label myGlobal = mergedGlobal[uncompactGlobal];
 486
 487                         if (myGlobal != nbrRegions[i])
 488                         {
 489                             label minRegion = min(myGlobal, nbrRegions[i]);
 490
 491                             if (debug)
 492                             {
 493                                 Pout<< "Merging region " << myGlobal
 494                                     << " (on proc " << Pstream::myProcNo()
 495                                     << ") and region " << nbrRegions[i]
 496                                     << " (on proc " << procPp.neighbProcNo()
 497                                     << ") into region " << minRegion << endl;
 498                             }
 499
 500                             mergedGlobal[uncompactGlobal] = minRegion;
 501                             mergedGlobal[myGlobal] = minRegion;
 502                             mergedGlobal[nbrRegions[i]] = minRegion;
 503
 504                             nMerged++;
 505                         }
 506                     }
 507
 508                     faceI++;
 509                 }
 510             }
 511         }
 512
 513
 514         reduce(nMerged, sumOp<label>());
 515
 516         if (debug)
 517         {
 518             Pout<< "nMerged:" << nMerged << endl;
 519         }
 520
 521         if (nMerged == 0)
 522         {
 523             break;
 524         }
 525
 526         // Merge the compacted regions.
 527         Pstream::listCombineGather(mergedGlobal, minEqOp<label>());
 528         Pstream::listCombineScatter(mergedGlobal);
 529     }
 530
 531
 532     // Compact global regions
 533     // ~~~~~~~~~~~~~~~~~~~~~~
 534
 535     // All procs will have the same global mergedGlobal region.
 536     // There might be gaps in it however so compact.
 537
 538     labelList mergedToCompacted(globalRegions.size(), -1);
 539
 540     label compactI = 0;
 541
 542     forAll(mergedGlobal, i)
 543     {
 544         label merged = mergedGlobal[i];
 545
 546         if (mergedToCompacted[merged] == -1)
 547         {
 548             mergedToCompacted[merged] = compactI++;
 549         }
 550     }
 551
 552     if (debug)
 553     {
 554         Pout<< "Compacted down to " << compactI << " regions." << endl;
 555     }
 556
 557     // Renumber cellRegion to be global regions
 558     forAll(cellRegion, cellI)
 559     {
 560         label region = cellRegion[cellI];
 561
 562         if (region >= 0)
 563         {
 564             label merged = mergedGlobal[globalRegions.toGlobal(region)];
 565
 566             cellRegion[cellI] = mergedToCompacted[merged];
 567         }
 568     }
 569
 570     return compactI;
 571 }
 572
 573
 574 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 575
 576 Foam::regionSplit::regionSplit(const polyMesh& mesh)
 577 :
 578     labelList(mesh.nCells(), -1),
 579     mesh_(mesh),
 580     nRegions_(calcRegionSplit(boolList(0, false), List<labelPair>(0), *this))
 581 {}
 582
 583
 584 Foam::regionSplit::regionSplit
 585 (
 586     const polyMesh& mesh,
 587     const boolList& blockedFace
 588 )
 589 :
 590     labelList(mesh.nCells(), -1),
 591     mesh_(mesh),
 592     nRegions_(calcRegionSplit(blockedFace, List<labelPair>(0), *this))
 593 {}
 594
 595
 596 Foam::regionSplit::regionSplit
 597 (
 598     const polyMesh& mesh,
 599     const boolList& blockedFace,
 600     const List<labelPair>& explicitConnections
 601 )
 602 :
 603     labelList(mesh.nCells(), -1),
 604     mesh_(mesh),
 605     nRegions_(calcRegionSplit(blockedFace, explicitConnections, *this))
 606 {}
 607
 608
 609 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 610
 611
 612 // ************************************************************************* //