src/meshTools/regionSplit/regionSplit.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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 "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     DynamicList<label>& 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         DynamicList<label> 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         DynamicList<label> 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
 213             (
 214                 isA<cyclicPolyPatch>(pp)
 215              && refCast<const cyclicPolyPatch>(pp).owner()
 216             )
 217             {
 218                 // Transfer from neighbourPatch to here or vice versa.
 219
 220                 const cyclicPolyPatch& cycPatch =
 221                     refCast<const cyclicPolyPatch>(pp);
 222
 223                 label faceI = cycPatch.start();
 224
 225                 forAll(cycPatch, i)
 226                 {
 227                     label otherFaceI = cycPatch.transformGlobalFace(faceI);
 228
 229                     transferCoupledFaceRegion
 230                     (
 231                         faceI,
 232                         otherFaceI,
 233                         faceRegion,
 234                         newChangedFaces
 235                     );
 236
 237                     faceI++;
 238                 }
 239             }
 240         }
 241         forAll(explicitConnections, i)
 242         {
 243             transferCoupledFaceRegion
 244             (
 245                 explicitConnections[i][0],
 246                 explicitConnections[i][1],
 247                 faceRegion,
 248                 newChangedFaces
 249             );
 250         }
 251
 252         //if (debug)
 253         //{
 254         //    Pout<< "regionSplit::fillSeedMask : changedFaces after sync:"
 255         //        << newChangedFaces.size() << endl;
 256         //}
 257
 258         changedFaces.transfer(newChangedFaces);
 259     }
 260 }
 261
 262
 263 Foam::label Foam::regionSplit::calcRegionSplit
 264 (
 265     const boolList& blockedFace,
 266     const List<labelPair>& explicitConnections,
 267
 268     labelList& cellRegion
 269 ) const
 270 {
 271     if (debug)
 272     {
 273         if (blockedFace.size())
 274         {
 275             // Check that blockedFace is synced.
 276             boolList syncBlockedFace(blockedFace);
 277             syncTools::swapFaceList(mesh_, syncBlockedFace);
 278
 279             forAll(syncBlockedFace, faceI)
 280             {
 281                 if (syncBlockedFace[faceI] != blockedFace[faceI])
 282                 {
 283                     FatalErrorIn
 284                     (
 285                         "regionSplit::calcRegionSplit(..)"
 286                     )   << "Face " << faceI << " not synchronised. My value:"
 287                         << blockedFace[faceI] << "  coupled value:"
 288                         << syncBlockedFace[faceI]
 289                         << abort(FatalError);
 290                 }
 291             }
 292         }
 293     }
 294
 295     // Region per face.
 296     // -1 unassigned
 297     // -2 blocked
 298     labelList faceRegion(mesh_.nFaces(), -1);
 299
 300     if (blockedFace.size())
 301     {
 302         forAll(blockedFace, faceI)
 303         {
 304             if (blockedFace[faceI])
 305             {
 306                 faceRegion[faceI] = -2;
 307             }
 308         }
 309     }
 310
 311
 312     // Assign local regions
 313     // ~~~~~~~~~~~~~~~~~~~~
 314
 315     // Start with region 0
 316     label nRegions = 0;
 317
 318     label unsetCellI = 0;
 319
 320     do
 321     {
 322         // Find first unset cell
 323
 324         for (; unsetCellI < mesh_.nCells(); unsetCellI++)
 325         {
 326             if (cellRegion[unsetCellI] == -1)
 327             {
 328                 break;
 329             }
 330         }
 331
 332         if (unsetCellI >= mesh_.nCells())
 333         {
 334             break;
 335         }
 336
 337         fillSeedMask
 338         (
 339             explicitConnections,
 340             cellRegion,
 341             faceRegion,
 342             unsetCellI,
 343             nRegions
 344         );
 345
 346         // Current unsetCell has now been handled. Go to next region.
 347         nRegions++;
 348         unsetCellI++;
 349     }
 350     while (true);
 351
 352
 353     if (debug)
 354     {
 355         forAll(cellRegion, cellI)
 356         {
 357             if (cellRegion[cellI] < 0)
 358             {
 359                 FatalErrorIn("regionSplit::calcRegionSplit(..)")
 360                     << "cell:" << cellI << " region:" << cellRegion[cellI]
 361                     << abort(FatalError);
 362             }
 363         }
 364
 365         forAll(faceRegion, faceI)
 366         {
 367             if (faceRegion[faceI] == -1)
 368             {
 369                 FatalErrorIn("regionSplit::calcRegionSplit(..)")
 370                     << "face:" << faceI << " region:" << faceRegion[faceI]
 371                     << abort(FatalError);
 372             }
 373         }
 374     }
 375
 376
 377
 378     // Assign global regions
 379     // ~~~~~~~~~~~~~~~~~~~~~
 380     // Offset local regions to create unique global regions.
 381
 382     globalIndex globalRegions(nRegions);
 383
 384
 385     // Merge global regions
 386     // ~~~~~~~~~~~~~~~~~~~~
 387     // Regions across non-blocked proc patches get merged.
 388     // This will set merged global regions to be the min of both.
 389     // (this will create gaps in the global region list so they will get
 390     // merged later on)
 391
 392     // Map from global to merged global
 393     labelList mergedGlobal(identity(globalRegions.size()));
 394
 395
 396     // See if any regions get merged. Only nessecary for parallel
 397     while (Pstream::parRun())
 398     {
 399         if (debug)
 400         {
 401             Pout<< nl << "-- Starting Iteration --" << endl;
 402         }
 403
 404         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 405
 406         // Send global regions across (or -2 if blocked face)
 407         forAll(patches, patchI)
 408         {
 409             const polyPatch& pp = patches[patchI];
 410
 411             if (isA<processorPolyPatch>(pp))
 412             {
 413                 labelList myGlobalRegions(pp.size());
 414
 415                 label faceI = pp.start();
 416
 417                 forAll(pp, i)
 418                 {
 419                     if (faceRegion[faceI] < 0)
 420                     {
 421                         myGlobalRegions[i] = faceRegion[faceI];
 422                     }
 423                     else
 424                     {
 425                         myGlobalRegions[i] = mergedGlobal
 426                         [globalRegions.toGlobal(faceRegion[faceI])];
 427                     }
 428
 429                     faceI++;
 430                 }
 431
 432                 OPstream toProcNbr
 433                 (
 434                     Pstream::blocking,
 435                     refCast<const processorPolyPatch>(pp).neighbProcNo()
 436                 );
 437
 438                 toProcNbr << myGlobalRegions;
 439             }
 440         }
 441
 442
 443         // Receive global regions
 444
 445         label nMerged = 0;
 446
 447         forAll(patches, patchI)
 448         {
 449             const polyPatch& pp = patches[patchI];
 450
 451             if (isA<processorPolyPatch>(pp))
 452             {
 453                 const processorPolyPatch& procPp =
 454                     refCast<const processorPolyPatch>(pp);
 455
 456                 IPstream fromProcNbr(Pstream::blocking, procPp.neighbProcNo());
 457
 458                 labelList nbrRegions(fromProcNbr);
 459
 460
 461                 // Compare with my regions to see which get merged.
 462
 463                 label faceI = pp.start();
 464
 465                 forAll(pp, i)
 466                 {
 467                     if
 468                     (
 469                         faceRegion[faceI] < 0
 470                      || nbrRegions[i] < 0
 471                     )
 472                     {
 473                         if (faceRegion[faceI] != nbrRegions[i])
 474                         {
 475                             FatalErrorIn("regionSplit::calcRegionSplit(..)")
 476                                 << "On patch:" << pp.name()
 477                                 << " face:" << faceI
 478                                 << " my local region:" << faceRegion[faceI]
 479                                 << " neighbouring region:"
 480                                 << nbrRegions[i] << nl
 481                                 << "Maybe your blockedFaces are not"
 482                                 << " synchronized across coupled faces?"
 483                                 << abort(FatalError);
 484                         }
 485                     }
 486                     else
 487                     {
 488                         label uncompactGlobal =
 489                             globalRegions.toGlobal(faceRegion[faceI]);
 490
 491                         label myGlobal = mergedGlobal[uncompactGlobal];
 492
 493                         if (myGlobal != nbrRegions[i])
 494                         {
 495                             label minRegion = min(myGlobal, nbrRegions[i]);
 496
 497                             if (debug)
 498                             {
 499                                 Pout<< "Merging region " << myGlobal
 500                                     << " (on proc " << Pstream::myProcNo()
 501                                     << ") and region " << nbrRegions[i]
 502                                     << " (on proc " << procPp.neighbProcNo()
 503                                     << ") into region " << minRegion << endl;
 504                             }
 505
 506                             mergedGlobal[uncompactGlobal] = minRegion;
 507                             mergedGlobal[myGlobal] = minRegion;
 508                             mergedGlobal[nbrRegions[i]] = minRegion;
 509
 510                             nMerged++;
 511                         }
 512                     }
 513
 514                     faceI++;
 515                 }
 516             }
 517         }
 518
 519
 520         reduce(nMerged, sumOp<label>());
 521
 522         if (debug)
 523         {
 524             Pout<< "nMerged:" << nMerged << endl;
 525         }
 526
 527         if (nMerged == 0)
 528         {
 529             break;
 530         }
 531
 532         // Merge the compacted regions.
 533         Pstream::listCombineGather(mergedGlobal, minEqOp<label>());
 534         Pstream::listCombineScatter(mergedGlobal);
 535     }
 536
 537
 538     // Compact global regions
 539     // ~~~~~~~~~~~~~~~~~~~~~~
 540
 541     // All procs will have the same global mergedGlobal region.
 542     // There might be gaps in it however so compact.
 543
 544     labelList mergedToCompacted(globalRegions.size(), -1);
 545
 546     label compactI = 0;
 547
 548     forAll(mergedGlobal, i)
 549     {
 550         label merged = mergedGlobal[i];
 551
 552         if (mergedToCompacted[merged] == -1)
 553         {
 554             mergedToCompacted[merged] = compactI++;
 555         }
 556     }
 557
 558     if (debug)
 559     {
 560         Pout<< "Compacted down to " << compactI << " regions." << endl;
 561     }
 562
 563     // Renumber cellRegion to be global regions
 564     forAll(cellRegion, cellI)
 565     {
 566         label region = cellRegion[cellI];
 567
 568         if (region >= 0)
 569         {
 570             label merged = mergedGlobal[globalRegions.toGlobal(region)];
 571
 572             cellRegion[cellI] = mergedToCompacted[merged];
 573         }
 574     }
 575
 576     return compactI;
 577 }
 578
 579
 580 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 581
 582 Foam::regionSplit::regionSplit(const polyMesh& mesh)
 583 :
 584     labelList(mesh.nCells(), -1),
 585     mesh_(mesh),
 586     nRegions_(calcRegionSplit(boolList(0, false), List<labelPair>(0), *this))
 587 {}
 588
 589
 590 Foam::regionSplit::regionSplit
 591 (
 592     const polyMesh& mesh,
 593     const boolList& blockedFace
 594 )
 595 :
 596     labelList(mesh.nCells(), -1),
 597     mesh_(mesh),
 598     nRegions_(calcRegionSplit(blockedFace, List<labelPair>(0), *this))
 599 {}
 600
 601
 602 Foam::regionSplit::regionSplit
 603 (
 604     const polyMesh& mesh,
 605     const boolList& blockedFace,
 606     const List<labelPair>& explicitConnections
 607 )
 608 :
 609     labelList(mesh.nCells(), -1),
 610     mesh_(mesh),
 611     nRegions_(calcRegionSplit(blockedFace, explicitConnections, *this))
 612 {}
 613
 614
 615 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 616
 617
 618 // ************************************************************************* //