src/meshTools/regionSplit/regionSplit.C

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