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