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