src/foam/algorithms/MeshWave/FaceCellWave.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 "FaceCellWave.H"
  27 #include "polyMesh.H"
  28 #include "processorPolyPatch.H"
  29 #include "cyclicPolyPatch.H"
  30 #include "OPstream.H"
  31 #include "IPstream.H"
  32 #include "PstreamReduceOps.H"
  33 #include "debug.H"
  34 #include "typeInfo.H"
  35
  36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  37
  38 template <class Type>
  39 Foam::debug::tolerancesSwitch Foam::FaceCellWave<Type>::geomTol_
  40 (
  41     "FaceCellWaveGeomTol",
  42     1e-6
  43 );
  44
  45 template <class Type>
  46 Foam::debug::tolerancesSwitch Foam::FaceCellWave<Type>::propagationTol_
  47 (
  48     "FaceCellWavePropagationTol",
  49     0.01
  50 );
  51
  52 // Write to ostream
  53 template <class Type>
  54 Foam::Ostream& Foam::FaceCellWave<Type>::writeFaces
  55 (
  56     const label nFaces,
  57     const labelList& faceLabels,
  58     const List<Type>& faceInfo,
  59     Ostream& os
  60 )
  61 {
  62     // Write list contents depending on data format
  63     if (os.format() == IOstream::ASCII)
  64     {
  65         os << nFaces;
  66
  67         for(label i = 0; i < nFaces; i++)
  68         {
  69             os << ' ' << faceLabels[i];
  70         }
  71         for(label i = 0; i < nFaces; i++)
  72         {
  73             os << ' ' << faceInfo[i];
  74         }
  75     }
  76     else
  77     {
  78         os << nFaces;
  79
  80         for(label i = 0; i < nFaces; i++)
  81         {
  82             os << faceLabels[i];
  83         }
  84         for(label i = 0; i < nFaces; i++)
  85         {
  86             os << faceInfo[i];
  87         }
  88     }
  89     return os;
  90 }
  91
  92
  93 // Read from istream
  94 template <class Type>
  95 Foam::Istream& Foam::FaceCellWave<Type>::readFaces
  96 (
  97     label& nFaces,
  98     labelList& faceLabels,
  99     List<Type>& faceInfo,
 100     Istream& is
 101 )
 102 {
 103     is >> nFaces;
 104
 105     for(label i = 0; i < nFaces; i++)
 106     {
 107         is >> faceLabels[i];
 108     }
 109     for(label i = 0; i < nFaces; i++)
 110     {
 111         is >> faceInfo[i];
 112     }
 113     return is;
 114 }
 115
 116
 117 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 118
 119 // Update info for cellI, at position pt, with information from
 120 // neighbouring face/cell.
 121 // Updates:
 122 //      - changedCell_, changedCells_, nChangedCells_,
 123 //      - statistics: nEvals_, nUnvisitedCells_
 124 template <class Type>
 125 bool Foam::FaceCellWave<Type>::updateCell
 126 (
 127     const label cellI,
 128     const label neighbourFaceI,
 129     const Type& neighbourInfo,
 130     const scalar tol,
 131     Type& cellInfo
 132 )
 133 {
 134     nEvals_++;
 135
 136     bool wasValid = cellInfo.valid();
 137
 138     bool propagate =
 139         cellInfo.updateCell
 140         (
 141             mesh_,
 142             cellI,
 143             neighbourFaceI,
 144             neighbourInfo,
 145             tol
 146         );
 147
 148     if (propagate)
 149     {
 150         if (!changedCell_[cellI])
 151         {
 152             changedCell_[cellI] = true;
 153             changedCells_[nChangedCells_++] = cellI;
 154         }
 155     }
 156
 157     if (!wasValid && cellInfo.valid())
 158     {
 159         --nUnvisitedCells_;
 160     }
 161
 162     return propagate;
 163 }
 164
 165
 166 // Update info for faceI, at position pt, with information from
 167 // neighbouring face/cell.
 168 // Updates:
 169 //      - changedFace_, changedFaces_, nChangedFaces_,
 170 //      - statistics: nEvals_, nUnvisitedFaces_
 171 template <class Type>
 172 bool Foam::FaceCellWave<Type>::updateFace
 173 (
 174     const label faceI,
 175     const label neighbourCellI,
 176     const Type& neighbourInfo,
 177     const scalar tol,
 178     Type& faceInfo
 179 )
 180 {
 181     nEvals_++;
 182
 183     bool wasValid = faceInfo.valid();
 184
 185     bool propagate =
 186         faceInfo.updateFace
 187         (
 188             mesh_,
 189             faceI,
 190             neighbourCellI,
 191             neighbourInfo,
 192             tol
 193         );
 194
 195     if (propagate)
 196     {
 197         if (!changedFace_[faceI])
 198         {
 199             changedFace_[faceI] = true;
 200             changedFaces_[nChangedFaces_++] = faceI;
 201         }
 202     }
 203
 204     if (!wasValid && faceInfo.valid())
 205     {
 206         --nUnvisitedFaces_;
 207     }
 208
 209     return propagate;
 210 }
 211
 212
 213 // Update info for faceI, at position pt, with information from
 214 // same face.
 215 // Updates:
 216 //      - changedFace_, changedFaces_, nChangedFaces_,
 217 //      - statistics: nEvals_, nUnvisitedFaces_
 218 template <class Type>
 219 bool Foam::FaceCellWave<Type>::updateFace
 220 (
 221     const label faceI,
 222     const Type& neighbourInfo,
 223     const scalar tol,
 224     Type& faceInfo
 225 )
 226 {
 227     nEvals_++;
 228
 229     bool wasValid = faceInfo.valid();
 230
 231     bool propagate =
 232         faceInfo.updateFace
 233         (
 234             mesh_,
 235             faceI,
 236             neighbourInfo,
 237             tol
 238         );
 239
 240     if (propagate)
 241     {
 242         if (!changedFace_[faceI])
 243         {
 244             changedFace_[faceI] = true;
 245             changedFaces_[nChangedFaces_++] = faceI;
 246         }
 247     }
 248
 249     if (!wasValid && faceInfo.valid())
 250     {
 251         --nUnvisitedFaces_;
 252     }
 253
 254     return propagate;
 255 }
 256
 257
 258 // For debugging: check status on both sides of cyclic
 259 template <class Type>
 260 void Foam::FaceCellWave<Type>::checkCyclic(const polyPatch& patch) const
 261 {
 262     label cycOffset = patch.size()/2;
 263
 264     for(label patchFaceI = 0; patchFaceI < cycOffset; patchFaceI++)
 265     {
 266         label i1 = patch.start() + patchFaceI;
 267         label i2 = i1 + cycOffset;
 268
 269         if (!allFaceInfo_[i1].sameGeometry(mesh_, allFaceInfo_[i2], geomTol_()))
 270         {
 271             FatalErrorIn("FaceCellWave<Type>::checkCyclic(const polyPatch&)")
 272                 << "problem: i:" << i1 << "  otheri:" << i2
 273                 << "   faceInfo:" << allFaceInfo_[i1]
 274                 << "   otherfaceInfo:" << allFaceInfo_[i2]
 275                 << abort(FatalError);
 276         }
 277
 278         if (changedFace_[i1] != changedFace_[i2])
 279         {
 280             FatalErrorIn("FaceCellWave<Type>::checkCyclic(const polyPatch&)")
 281                 << " problem: i:" << i1 << "  otheri:" << i2
 282                 << "   faceInfo:" << allFaceInfo_[i1]
 283                 << "   otherfaceInfo:" << allFaceInfo_[i2]
 284                 << "   changedFace:" << changedFace_[i1]
 285                 << "   otherchangedFace:" << changedFace_[i2]
 286                 << abort(FatalError);
 287         }
 288     }
 289 }
 290
 291
 292 // Check if patches of given type name are present
 293 template <class Type>
 294 bool Foam::FaceCellWave<Type>::hasPatchType(const word& nameOfType)
 295 {
 296     forAll(mesh_.boundaryMesh(), patchI)
 297     {
 298         if (mesh_.boundaryMesh()[patchI].type() == nameOfType)
 299         {
 300             return true;
 301         }
 302     }
 303     return false;
 304 }
 305
 306
 307 // Copy face information into member data
 308 template <class Type>
 309 void Foam::FaceCellWave<Type>::setFaceInfo
 310 (
 311     const labelList& changedFaces,
 312     const List<Type>& changedFacesInfo
 313 )
 314 {
 315     forAll(changedFaces, changedFaceI)
 316     {
 317         label faceI = changedFaces[changedFaceI];
 318
 319         bool wasValid = allFaceInfo_[faceI].valid();
 320
 321         // Copy info for faceI
 322         allFaceInfo_[faceI] = changedFacesInfo[changedFaceI];
 323
 324         // Maintain count of unset faces
 325         if (!wasValid && allFaceInfo_[faceI].valid())
 326         {
 327             --nUnvisitedFaces_;
 328         }
 329
 330         // Mark faceI as changed, both on list and on face itself.
 331
 332         changedFace_[faceI] = true;
 333         changedFaces_[nChangedFaces_++] = faceI;
 334     }
 335 }
 336
 337
 338 // Merge face information into member data
 339 template <class Type>
 340 void Foam::FaceCellWave<Type>::mergeFaceInfo
 341 (
 342     const polyPatch& patch,
 343     const label nFaces,
 344     const labelList& changedFaces,
 345     const List<Type>& changedFacesInfo,
 346     const bool
 347 )
 348 {
 349     for(label changedFaceI = 0; changedFaceI < nFaces; changedFaceI++)
 350     {
 351         const Type& neighbourWallInfo = changedFacesInfo[changedFaceI];
 352         label patchFaceI = changedFaces[changedFaceI];
 353
 354         label meshFaceI = patch.start() + patchFaceI;
 355
 356         Type& currentWallInfo = allFaceInfo_[meshFaceI];
 357
 358         if (currentWallInfo != neighbourWallInfo)
 359         {
 360             updateFace
 361             (
 362                 meshFaceI,
 363                 neighbourWallInfo,
 364                 propagationTol_(),
 365                 currentWallInfo
 366             );
 367         }
 368     }
 369 }
 370
 371
 372 // Construct compact patchFace change arrays for a (slice of a) single patch.
 373 // changedPatchFaces in local patch numbering.
 374 // Return length of arrays.
 375 template <class Type>
 376 Foam::label Foam::FaceCellWave<Type>::getChangedPatchFaces
 377 (
 378     const polyPatch& patch,
 379     const label startFaceI,
 380     const label nFaces,
 381     labelList& changedPatchFaces,
 382     List<Type>& changedPatchFacesInfo
 383 ) const
 384 {
 385     label nChangedPatchFaces = 0;
 386
 387     for(label i = 0; i < nFaces; i++)
 388     {
 389         label patchFaceI = i + startFaceI;
 390
 391         label meshFaceI = patch.start() + patchFaceI;
 392
 393         if (changedFace_[meshFaceI])
 394         {
 395             changedPatchFaces[nChangedPatchFaces] = patchFaceI;
 396             changedPatchFacesInfo[nChangedPatchFaces] =
 397                 allFaceInfo_[meshFaceI];
 398             nChangedPatchFaces++;
 399         }
 400     }
 401     return nChangedPatchFaces;
 402 }
 403
 404
 405 // Handle leaving domain. Implementation referred to Type
 406 template <class Type>
 407 void Foam::FaceCellWave<Type>::leaveDomain
 408 (
 409     const polyPatch& patch,
 410     const label nFaces,
 411     const labelList& faceLabels,
 412     List<Type>& faceInfo
 413 ) const
 414 {
 415     const vectorField& fc = mesh_.faceCentres();
 416
 417     for(label i = 0; i < nFaces; i++)
 418     {
 419         label patchFaceI = faceLabels[i];
 420
 421         label meshFaceI = patch.start() + patchFaceI;
 422         faceInfo[i].leaveDomain(mesh_, patch, patchFaceI, fc[meshFaceI]);
 423     }
 424 }
 425
 426
 427 // Handle entering domain. Implementation referred to Type
 428 template <class Type>
 429 void Foam::FaceCellWave<Type>::enterDomain
 430 (
 431     const polyPatch& patch,
 432     const label nFaces,
 433     const labelList& faceLabels,
 434     List<Type>& faceInfo
 435 ) const
 436 {
 437     const vectorField& fc = mesh_.faceCentres();
 438
 439     for(label i = 0; i < nFaces; i++)
 440     {
 441         label patchFaceI = faceLabels[i];
 442
 443         label meshFaceI = patch.start() + patchFaceI;
 444         faceInfo[i].enterDomain(mesh_, patch, patchFaceI, fc[meshFaceI]);
 445     }
 446 }
 447
 448
 449 // Transform. Implementation referred to Type
 450 template <class Type>
 451 void Foam::FaceCellWave<Type>::transform
 452 (
 453     const tensorField& rotTensor,
 454     const label nFaces,
 455     List<Type>& faceInfo
 456 )
 457 {
 458     if (rotTensor.size() == 1)
 459     {
 460         const tensor& T = rotTensor[0];
 461
 462         for(label faceI = 0; faceI < nFaces; faceI++)
 463         {
 464             faceInfo[faceI].transform(mesh_, T);
 465         }
 466     }
 467     else
 468     {
 469         for(label faceI = 0; faceI < nFaces; faceI++)
 470         {
 471             faceInfo[faceI].transform(mesh_, rotTensor[faceI]);
 472         }
 473     }
 474 }
 475
 476
 477 // Send face info to neighbour.
 478 template <class Type>
 479 void Foam::FaceCellWave<Type>::sendPatchInfo
 480 (
 481     const label neighbour,
 482     const label nFaces,
 483     const labelList& faceLabels,
 484     const List<Type>& faceInfo
 485 ) const
 486 {
 487     OPstream toNeighbour(Pstream::blocking, neighbour);
 488
 489     writeFaces(nFaces, faceLabels, faceInfo, toNeighbour);
 490 }
 491
 492
 493 // Receive face info from neighbour
 494 template <class Type>
 495 Foam::label Foam::FaceCellWave<Type>::receivePatchInfo
 496 (
 497     const label neighbour,
 498     labelList& faceLabels,
 499     List<Type>& faceInfo
 500 ) const
 501 {
 502     IPstream fromNeighbour(Pstream::blocking, neighbour);
 503
 504     label nFaces = 0;
 505     readFaces(nFaces, faceLabels, faceInfo, fromNeighbour);
 506
 507     return nFaces;
 508 }
 509
 510
 511 // Offset mesh face. Used for transferring from one cyclic half to the other.
 512 template <class Type>
 513 void Foam::FaceCellWave<Type>::offset
 514 (
 515     const polyPatch&,
 516     const label cycOffset,
 517     const label nFaces,
 518     labelList& faces
 519 )
 520 {
 521     for(label faceI = 0; faceI < nFaces; faceI++)
 522     {
 523         faces[faceI] += cycOffset;
 524     }
 525 }
 526
 527
 528 // Tranfer all the information to/from neighbouring processors
 529 template <class Type>
 530 void Foam::FaceCellWave<Type>::handleProcPatches()
 531 {
 532     // Send all
 533
 534     forAll(mesh_.boundaryMesh(), patchI)
 535     {
 536         const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 537
 538         if (isA<processorPolyPatch>(patch))
 539         {
 540             // Allocate buffers
 541             label nSendFaces;
 542             labelList sendFaces(patch.size());
 543             List<Type> sendFacesInfo(patch.size());
 544
 545             // Determine which faces changed on current patch
 546             nSendFaces = getChangedPatchFaces
 547             (
 548                 patch,
 549                 0,
 550                 patch.size(),
 551                 sendFaces,
 552                 sendFacesInfo
 553             );
 554
 555             // Adapt wallInfo for leaving domain
 556             leaveDomain
 557             (
 558                 patch,
 559                 nSendFaces,
 560                 sendFaces,
 561                 sendFacesInfo
 562             );
 563
 564             const processorPolyPatch& procPatch =
 565                 refCast<const processorPolyPatch>(patch);
 566
 567             if (debug)
 568             {
 569                 Pout<< " Processor patch " << patchI << ' ' << patch.name()
 570                     << " communicating with " << procPatch.neighbProcNo()
 571                     << "  Sending:" << nSendFaces
 572                     << endl;
 573             }
 574
 575             sendPatchInfo
 576             (
 577                 procPatch.neighbProcNo(),
 578                 nSendFaces,
 579                 sendFaces,
 580                 sendFacesInfo
 581             );
 582         }
 583     }
 584
 585     // Receive all
 586
 587     forAll(mesh_.boundaryMesh(), patchI)
 588     {
 589         const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 590
 591         if (isA<processorPolyPatch>(patch))
 592         {
 593             const processorPolyPatch& procPatch =
 594                 refCast<const processorPolyPatch>(patch);
 595
 596             // Allocate buffers
 597             label nReceiveFaces;
 598             labelList receiveFaces(patch.size());
 599             List<Type> receiveFacesInfo(patch.size());
 600
 601             nReceiveFaces = receivePatchInfo
 602             (
 603                 procPatch.neighbProcNo(),
 604                 receiveFaces,
 605                 receiveFacesInfo
 606             );
 607
 608             if (debug)
 609             {
 610                 Pout<< " Processor patch " << patchI << ' ' << patch.name()
 611                     << " communicating with " << procPatch.neighbProcNo()
 612                     << "  Receiving:" << nReceiveFaces
 613                     << endl;
 614             }
 615
 616             // Apply transform to received data for non-parallel planes
 617             if (!procPatch.parallel())
 618             {
 619                 transform
 620                 (
 621                     procPatch.forwardT(),
 622                     nReceiveFaces,
 623                     receiveFacesInfo
 624                 );
 625             }
 626
 627             // Adapt wallInfo for entering domain
 628             enterDomain
 629             (
 630                 patch,
 631                 nReceiveFaces,
 632                 receiveFaces,
 633                 receiveFacesInfo
 634             );
 635
 636             // Merge received info
 637             mergeFaceInfo
 638             (
 639                 patch,
 640                 nReceiveFaces,
 641                 receiveFaces,
 642                 receiveFacesInfo,
 643                 procPatch.parallel()
 644             );
 645         }
 646     }
 647 }
 648
 649
 650 // Transfer information across cyclic halves.
 651 // Note: Cyclic is two patches in one: one side from 0..size/2 and other
 652 // side from size/2 .. size.
 653 template <class Type>
 654 void Foam::FaceCellWave<Type>::handleCyclicPatches()
 655 {
 656     forAll(mesh_.boundaryMesh(), patchI)
 657     {
 658         const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 659
 660         if (isA<cyclicPolyPatch>(patch))
 661         {
 662             label halfSize = patch.size()/2;
 663
 664             // Allocate buffers
 665             label nSendFaces;
 666             labelList sendFaces(halfSize);
 667             List<Type> sendFacesInfo(halfSize);
 668
 669             label nReceiveFaces;
 670             labelList receiveFaces(halfSize);
 671             List<Type> receiveFacesInfo(halfSize);
 672
 673             // Half1: Determine which faces changed. Use sendFaces for storage
 674             nSendFaces = getChangedPatchFaces
 675             (
 676                 patch,
 677                 0,
 678                 halfSize,
 679                 sendFaces,
 680                 sendFacesInfo
 681             );
 682
 683             // Half2: Determine which faces changed. Use receiveFaces_  ,,
 684             nReceiveFaces = getChangedPatchFaces
 685             (
 686                 patch,
 687                 halfSize,
 688                 halfSize,
 689                 receiveFaces,
 690                 receiveFacesInfo
 691             );
 692
 693             //Info<< "Half1:" << endl;
 694             //writeFaces(nSendFaces, sendFaces, sendFacesInfo, Info);
 695             //Info<< endl;
 696             //
 697             //Info<< "Half2:" << endl;
 698             //writeFaces(nReceiveFaces, receiveFaces, receiveFacesInfo, Info);
 699             //Info<< endl;
 700
 701
 702             // Half1: Adapt wallInfo for leaving domain
 703             leaveDomain
 704             (
 705                 patch,
 706                 nSendFaces,
 707                 sendFaces,
 708                 sendFacesInfo
 709             );
 710             // Half2: Adapt wallInfo for leaving domain
 711             leaveDomain
 712             (
 713                 patch,
 714                 nReceiveFaces,
 715                 receiveFaces,
 716                 receiveFacesInfo
 717             );
 718
 719             // Half1: 'transfer' to other side by offsetting patchFaceI
 720             offset(patch, halfSize, nSendFaces, sendFaces);
 721
 722             // Half2: 'transfer' to other side
 723             offset(patch, -halfSize, nReceiveFaces, receiveFaces);
 724
 725             // Apply rotation for non-parallel planes
 726             const cyclicPolyPatch& cycPatch =
 727                 refCast<const cyclicPolyPatch>(patch);
 728
 729             if (!cycPatch.parallel())
 730             {
 731                 // sendFaces = received data from half1
 732                 transform
 733                 (
 734                     cycPatch.forwardT(),
 735                     nSendFaces,
 736                     sendFacesInfo
 737                 );
 738
 739                 // receiveFaces = received data from half2
 740                 transform
 741                 (
 742                     cycPatch.forwardT(),
 743                     nReceiveFaces,
 744                     receiveFacesInfo
 745                 );
 746             }
 747
 748             if (debug)
 749             {
 750                 Pout<< " Cyclic patch " << patchI << ' ' << patch.name()
 751                     << "  Changed on first half : " << nSendFaces
 752                     << "  Changed on second half : " << nReceiveFaces
 753                     << endl;
 754             }
 755
 756             // Half1: Adapt wallInfo for entering domain
 757             enterDomain
 758             (
 759                 patch,
 760                 nSendFaces,
 761                 sendFaces,
 762                 sendFacesInfo
 763             );
 764
 765             // Half2: Adapt wallInfo for entering domain
 766             enterDomain
 767             (
 768                 patch,
 769                 nReceiveFaces,
 770                 receiveFaces,
 771                 receiveFacesInfo
 772             );
 773
 774             // Merge into global storage
 775             mergeFaceInfo
 776             (
 777                 patch,
 778                 nSendFaces,
 779                 sendFaces,
 780                 sendFacesInfo,
 781                 cycPatch.parallel()
 782             );
 783             // Merge into global storage
 784             mergeFaceInfo
 785             (
 786                 patch,
 787                 nReceiveFaces,
 788                 receiveFaces,
 789                 receiveFacesInfo,
 790                 cycPatch.parallel()
 791             );
 792
 793             if (debug)
 794             {
 795                 checkCyclic(patch);
 796             }
 797         }
 798     }
 799 }
 800
 801
 802 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 803
 804 // Set up only. Use setFaceInfo and iterate() to do actual calculation.
 805 template <class Type>
 806 Foam::FaceCellWave<Type>::FaceCellWave
 807 (
 808     const polyMesh& mesh,
 809     UList<Type>& allFaceInfo,
 810     UList<Type>& allCellInfo
 811 )
 812 :
 813     mesh_(mesh),
 814     allFaceInfo_(allFaceInfo),
 815     allCellInfo_(allCellInfo),
 816     changedFace_(mesh_.nFaces(), false),
 817     changedFaces_(mesh_.nFaces()),
 818     nChangedFaces_(0),
 819     changedCell_(mesh_.nCells(), false),
 820     changedCells_(mesh_.nCells()),
 821     nChangedCells_(0),
 822     hasCyclicPatches_(hasPatchType(cyclicPolyPatch::typeName)),
 823     nEvals_(0),
 824     nUnvisitedCells_(mesh_.nCells()),
 825     nUnvisitedFaces_(mesh_.nFaces()),
 826     iter_(0)
 827 {}
 828
 829
 830 // Iterate, propagating changedFacesInfo across mesh, until no change (or
 831 // maxIter reached). Initial cell values specified.
 832 template <class Type>
 833 Foam::FaceCellWave<Type>::FaceCellWave
 834 (
 835     const polyMesh& mesh,
 836     const labelList& changedFaces,
 837     const List<Type>& changedFacesInfo,
 838     UList<Type>& allFaceInfo,
 839     UList<Type>& allCellInfo,
 840     const label maxIter
 841 )
 842 :
 843     mesh_(mesh),
 844     allFaceInfo_(allFaceInfo),
 845     allCellInfo_(allCellInfo),
 846     changedFace_(mesh_.nFaces(), false),
 847     changedFaces_(mesh_.nFaces()),
 848     nChangedFaces_(0),
 849     changedCell_(mesh_.nCells(), false),
 850     changedCells_(mesh_.nCells()),
 851     nChangedCells_(0),
 852     hasCyclicPatches_(hasPatchType(cyclicPolyPatch::typeName)),
 853     nEvals_(0),
 854     nUnvisitedCells_(mesh_.nCells()),
 855     nUnvisitedFaces_(mesh_.nFaces()),
 856     iter_(0)
 857 {
 858     // Copy initial changed faces data
 859     setFaceInfo(changedFaces, changedFacesInfo);
 860
 861     // Iterate until nothing changes
 862     iterate(maxIter);
 863
 864     if ((maxIter > 0) && (iter_ >= maxIter))
 865     {
 866         FatalErrorIn
 867         (
 868             "FaceCellWave<Type>::FaceCellWave"
 869             "(const polyMesh&, const labelList&, const List<Type>,"
 870             " UList<Type>&, UList<Type>&, const label maxIter)"
 871         )
 872             << "Maximum number of iterations reached. Increase maxIter." << nl
 873             << "    maxIter:" << maxIter << endl
 874             << "    nChangedCells:" << nChangedCells_ << endl
 875             << "    nChangedFaces:" << nChangedFaces_ << endl
 876             << exit(FatalError);
 877     }
 878 }
 879
 880
 881 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 882
 883
 884 template <class Type>
 885 Foam::label Foam::FaceCellWave<Type>::getUnsetCells() const
 886 {
 887     return nUnvisitedCells_;
 888 }
 889
 890
 891 template <class Type>
 892 Foam::label Foam::FaceCellWave<Type>::getUnsetFaces() const
 893 {
 894     return nUnvisitedFaces_;
 895 }
 896
 897
 898
 899 // Propagate cell to face
 900 template <class Type>
 901 Foam::label Foam::FaceCellWave<Type>::faceToCell()
 902 {
 903     const labelList& owner = mesh_.faceOwner();
 904     const labelList& neighbour = mesh_.faceNeighbour();
 905     label nInternalFaces = mesh_.nInternalFaces();
 906
 907     for
 908     (
 909         label changedFaceI = 0;
 910         changedFaceI < nChangedFaces_;
 911         changedFaceI++
 912     )
 913     {
 914         label faceI = changedFaces_[changedFaceI];
 915         if (!changedFace_[faceI])
 916         {
 917             FatalErrorIn("FaceCellWave<Type>::faceToCell()")
 918                 << "Face " << faceI
 919                 << " not marked as having been changed"
 920                 << abort(FatalError);
 921         }
 922
 923
 924         const Type& neighbourWallInfo = allFaceInfo_[faceI];
 925
 926         // Evaluate all connected cells
 927
 928         // Owner
 929         label cellI = owner[faceI];
 930         Type& currentWallInfo = allCellInfo_[cellI];
 931
 932         if (currentWallInfo != neighbourWallInfo)
 933         {
 934             updateCell
 935             (
 936                 cellI,
 937                 faceI,
 938                 neighbourWallInfo,
 939                 propagationTol_(),
 940                 currentWallInfo
 941             );
 942         }
 943
 944         // Neighbour. Hack for check if face has neighbour.
 945         if (faceI < nInternalFaces)
 946         {
 947             cellI = neighbour[faceI];
 948             Type& currentWallInfo2 = allCellInfo_[cellI];
 949
 950             if (currentWallInfo2 != neighbourWallInfo)
 951             {
 952                 updateCell
 953                 (
 954                     cellI,
 955                     faceI,
 956                     neighbourWallInfo,
 957                     propagationTol_(),
 958                     currentWallInfo2
 959                 );
 960             }
 961         }
 962
 963         // Reset status of face
 964         changedFace_[faceI] = false;
 965     }
 966
 967     // Handled all changed faces by now
 968     nChangedFaces_ = 0;
 969
 970     if (debug)
 971     {
 972         Pout<< " Changed cells            : " << nChangedCells_ << endl;
 973     }
 974
 975     // Sum nChangedCells over all procs
 976     label totNChanged = nChangedCells_;
 977
 978     reduce(totNChanged, sumOp<label>());
 979
 980     return totNChanged;
 981 }
 982
 983
 984 // Propagate cell to face
 985 template <class Type>
 986 Foam::label Foam::FaceCellWave<Type>::cellToFace()
 987 {
 988     const cellList& cells = mesh_.cells();
 989
 990     for
 991     (
 992         label changedCellI = 0;
 993         changedCellI < nChangedCells_;
 994         changedCellI++
 995     )
 996     {
 997         label cellI = changedCells_[changedCellI];
 998         if (!changedCell_[cellI])
 999         {
1000             FatalErrorIn("FaceCellWave<Type>::cellToFace()")
1001                 << "Cell " << cellI
1002                 << " not marked as having been changed"
1003                 << abort(FatalError);
1004         }
1005
1006         const Type& neighbourWallInfo = allCellInfo_[cellI];
1007
1008         // Evaluate all connected faces
1009
1010         const labelList& faceLabels = cells[cellI];
1011         forAll(faceLabels, faceLabelI)
1012         {
1013             label faceI = faceLabels[faceLabelI];
1014             Type& currentWallInfo = allFaceInfo_[faceI];
1015
1016             if (currentWallInfo != neighbourWallInfo)
1017             {
1018                 updateFace
1019                 (
1020                     faceI,
1021                     cellI,
1022                     neighbourWallInfo,
1023                     propagationTol_(),
1024                     currentWallInfo
1025                 );
1026             }
1027         }
1028
1029         // Reset status of cell
1030         changedCell_[cellI] = false;
1031     }
1032
1033     // Handled all changed cells by now
1034     nChangedCells_ = 0;
1035
1036     if (hasCyclicPatches_)
1037     {
1038         // Transfer changed faces across cyclic halves
1039         handleCyclicPatches();
1040     }
1041     if (Pstream::parRun())
1042     {
1043         // Transfer changed faces from neighbouring processors.
1044         handleProcPatches();
1045     }
1046
1047     if (debug)
1048     {
1049         Pout<< " Changed faces            : " << nChangedFaces_ << endl;
1050     }
1051
1052     // Sum nChangedFaces over all procs
1053     label totNChanged = nChangedFaces_;
1054
1055     reduce(totNChanged, sumOp<label>());
1056
1057     return totNChanged;
1058 }
1059
1060
1061 // Iterate
1062 template <class Type>
1063 Foam::label Foam::FaceCellWave<Type>::iterate(const label maxIter)
1064 {
1065     if (hasCyclicPatches_)
1066     {
1067         // Transfer changed faces across cyclic halves
1068         handleCyclicPatches();
1069     }
1070     if (Pstream::parRun())
1071     {
1072         // Transfer changed faces from neighbouring processors.
1073         handleProcPatches();
1074     }
1075
1076     while(iter_ < maxIter)
1077     {
1078         if (debug)
1079         {
1080             Pout<< " Iteration " << iter_ << endl;
1081         }
1082
1083         nEvals_ = 0;
1084
1085         label nCells = faceToCell();
1086
1087         if (debug)
1088         {
1089             Pout<< " Total changed cells      : " << nCells << endl;
1090         }
1091
1092         if (nCells == 0)
1093         {
1094             break;
1095         }
1096
1097         label nFaces = cellToFace();
1098
1099         if (debug)
1100         {
1101             Pout<< " Total changed faces      : " << nFaces << nl
1102                 << " Total evaluations        : " << nEvals_ << nl
1103                 << " Remaining unvisited cells: " << nUnvisitedCells_ << nl
1104                 << " Remaining unvisited faces: " << nUnvisitedFaces_ << endl;
1105         }
1106
1107         if (nFaces == 0)
1108         {
1109             break;
1110         }
1111
1112         ++iter_;
1113     }
1114
1115     return nUnvisitedCells_;
1116 }
1117
1118 // ************************************************************************* //