src/meshTools/PointEdgeWave/PointEdgeWave.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 "PointEdgeWave.H"
  27 #include "polyMesh.H"
  28 #include "processorPolyPatch.H"
  29 #include "cyclicPolyPatch.H"
  30 #include "ggiPolyPatch.H"
  31 #include "OPstream.H"
  32 #include "IPstream.H"
  33 #include "PstreamCombineReduceOps.H"
  34 #include "debug.H"
  35 #include "typeInfo.H"
  36
  37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  38
  39 template <class Type>
  40 Foam::scalar Foam::PointEdgeWave<Type>::propagationTol_ = 0.01;
  41
  42
  43 // Offset labelList. Used for transferring from one cyclic half to the other.
  44 template <class Type>
  45 void Foam::PointEdgeWave<Type>::offset(const label val, labelList& elems)
  46 {
  47     forAll(elems, i)
  48     {
  49         elems[i] += val;
  50     }
  51 }
  52
  53
  54 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  55
  56 // Gets point-point correspondence.  Is
  57 // - list of halfA points (in cyclic patch points)
  58 // - list of halfB points (can overlap with A!)
  59 // - for every patchPoint its corresponding point
  60 template <class Type>
  61 void Foam::PointEdgeWave<Type>::calcCyclicAddressing()
  62 {
  63     label cycHalf = 0;
  64
  65     forAll(mesh_.boundaryMesh(), patchI)
  66     {
  67         const polyPatch& patch = mesh_.boundaryMesh()[patchI];
  68
  69         if (isA<cyclicPolyPatch>(patch))
  70         {
  71             label halfSize = patch.size()/2;
  72
  73             SubList<face> halfAFaces
  74             (
  75                 mesh_.faces(),
  76                 halfSize,
  77                 patch.start()
  78             );
  79
  80             cycHalves_.set
  81             (
  82                 cycHalf++,
  83                 new primitivePatch(halfAFaces, mesh_.points())
  84             );
  85
  86             SubList<face> halfBFaces
  87             (
  88                 mesh_.faces(),
  89                 halfSize,
  90                 patch.start() + halfSize
  91             );
  92
  93             cycHalves_.set
  94             (
  95                 cycHalf++,
  96                 new primitivePatch(halfBFaces, mesh_.points())
  97             );
  98         }
  99     }
 100 }
 101
 102
 103 // Handle leaving domain. Implementation referred to Type
 104 template <class Type>
 105 void Foam::PointEdgeWave<Type>::leaveDomain
 106 (
 107     const polyPatch& meshPatch,
 108     const primitivePatch& patch,
 109     const labelList& patchPointLabels,
 110     List<Type>& pointInfo
 111 ) const
 112 {
 113     const labelList& meshPoints = patch.meshPoints();
 114
 115     forAll(patchPointLabels, i)
 116     {
 117         label patchPointI = patchPointLabels[i];
 118
 119         const point& pt = patch.points()[meshPoints[patchPointI]];
 120
 121         pointInfo[i].leaveDomain(meshPatch, patchPointI, pt);
 122     }
 123 }
 124
 125
 126 // Handle entering domain. Implementation referred to Type
 127 template <class Type>
 128 void Foam::PointEdgeWave<Type>::enterDomain
 129 (
 130     const polyPatch& meshPatch,
 131     const primitivePatch& patch,
 132     const labelList& patchPointLabels,
 133     List<Type>& pointInfo
 134 ) const
 135 {
 136     const labelList& meshPoints = patch.meshPoints();
 137
 138     forAll(patchPointLabels, i)
 139     {
 140         label patchPointI = patchPointLabels[i];
 141
 142         const point& pt = patch.points()[meshPoints[patchPointI]];
 143
 144         pointInfo[i].enterDomain(meshPatch, patchPointI, pt);
 145     }
 146 }
 147
 148
 149 // Transform. Implementation referred to Type
 150 template <class Type>
 151 void Foam::PointEdgeWave<Type>::transform
 152 (
 153     const tensorField& rotTensor,
 154     List<Type>& pointInfo
 155 ) const
 156 {
 157     if (rotTensor.size() == 1)
 158     {
 159         const tensor& T = rotTensor[0];
 160
 161         forAll(pointInfo, i)
 162         {
 163             pointInfo[i].transform(T);
 164         }
 165     }
 166     else
 167     {
 168         FatalErrorIn
 169         (
 170             "PointEdgeWave<Type>::transform(const tensorField&, List<Type>&)"
 171         )   << "Parallel cyclics not supported"
 172             << abort(FatalError);
 173
 174         forAll(pointInfo, i)
 175         {
 176             pointInfo[i].transform(rotTensor[i]);
 177         }
 178     }
 179 }
 180
 181
 182 // Update info for pointI, at position pt, with information from
 183 // neighbouring edge.
 184 // Updates:
 185 //      - changedPoint_, changedPoints_, nChangedPoints_,
 186 //      - statistics: nEvals_, nUnvisitedPoints_
 187 template <class Type>
 188 bool Foam::PointEdgeWave<Type>::updatePoint
 189 (
 190     const label pointI,
 191     const label neighbourEdgeI,
 192     const Type& neighbourInfo,
 193     const scalar tol,
 194     Type& pointInfo
 195 )
 196 {
 197     nEvals_++;
 198
 199     bool wasValid = pointInfo.valid();
 200
 201     bool propagate =
 202         pointInfo.updatePoint
 203         (
 204             mesh_,
 205             pointI,
 206             neighbourEdgeI,
 207             neighbourInfo,
 208             tol
 209         );
 210
 211     if (propagate)
 212     {
 213         if (!changedPoint_[pointI])
 214         {
 215             changedPoint_[pointI] = true;
 216             changedPoints_[nChangedPoints_++] = pointI;
 217         }
 218     }
 219
 220     if (!wasValid && pointInfo.valid())
 221     {
 222         --nUnvisitedPoints_;
 223     }
 224
 225     return propagate;
 226 }
 227
 228
 229 // Update info for pointI, at position pt, with information from
 230 // same point.
 231 // Updates:
 232 //      - changedPoint_, changedPoints_, nChangedPoints_,
 233 //      - statistics: nEvals_, nUnvisitedPoints_
 234 template <class Type>
 235 bool Foam::PointEdgeWave<Type>::updatePoint
 236 (
 237     const label pointI,
 238     const Type& neighbourInfo,
 239     const scalar tol,
 240     Type& pointInfo
 241 )
 242 {
 243     nEvals_++;
 244
 245     bool wasValid = pointInfo.valid();
 246
 247     bool propagate =
 248         pointInfo.updatePoint
 249         (
 250             mesh_,
 251             pointI,
 252             neighbourInfo,
 253             tol
 254         );
 255
 256     if (propagate)
 257     {
 258         if (!changedPoint_[pointI])
 259         {
 260             changedPoint_[pointI] = true;
 261             changedPoints_[nChangedPoints_++] = pointI;
 262         }
 263     }
 264
 265     if (!wasValid && pointInfo.valid())
 266     {
 267         --nUnvisitedPoints_;
 268     }
 269
 270     return propagate;
 271 }
 272
 273
 274 // Update info for edgeI, at position pt, with information from
 275 // neighbouring point.
 276 // Updates:
 277 //      - changedEdge_, changedEdges_, nChangedEdges_,
 278 //      - statistics: nEvals_, nUnvisitedEdge_
 279 template <class Type>
 280 bool Foam::PointEdgeWave<Type>::updateEdge
 281 (
 282     const label edgeI,
 283     const label neighbourPointI,
 284     const Type& neighbourInfo,
 285     const scalar tol,
 286     Type& edgeInfo
 287 )
 288 {
 289     nEvals_++;
 290
 291     bool wasValid = edgeInfo.valid();
 292
 293     bool propagate =
 294         edgeInfo.updateEdge
 295         (
 296             mesh_,
 297             edgeI,
 298             neighbourPointI,
 299             neighbourInfo,
 300             tol
 301         );
 302
 303     if (propagate)
 304     {
 305         if (!changedEdge_[edgeI])
 306         {
 307             changedEdge_[edgeI] = true;
 308             changedEdges_[nChangedEdges_++] = edgeI;
 309         }
 310     }
 311
 312     if (!wasValid && edgeInfo.valid())
 313     {
 314         --nUnvisitedEdges_;
 315     }
 316
 317     return propagate;
 318 }
 319
 320
 321 // Check if patches of given type name are present
 322 template <class Type>
 323 template <class PatchType>
 324 Foam::label Foam::PointEdgeWave<Type>::countPatchType() const
 325 {
 326     label nPatches = 0;
 327
 328     forAll(mesh_.boundaryMesh(), patchI)
 329     {
 330         if (isA<PatchType>(mesh_.boundaryMesh()[patchI]))
 331         {
 332             nPatches++;
 333         }
 334     }
 335     return nPatches;
 336 }
 337
 338
 339 // Collect changed patch points
 340 template <class Type>
 341 void Foam::PointEdgeWave<Type>::getChangedPatchPoints
 342 (
 343     const primitivePatch& patch,
 344
 345     DynamicList<Type>& patchInfo,
 346     dynamicLabelList& patchPoints,
 347     dynamicLabelList& owner,
 348     dynamicLabelList& ownerIndex
 349 ) const
 350 {
 351     const labelList& meshPoints = patch.meshPoints();
 352     const faceList& localFaces = patch.localFaces();
 353     const labelListList& pointFaces = patch.pointFaces();
 354
 355     forAll(meshPoints, patchPointI)
 356     {
 357         label meshPointI = meshPoints[patchPointI];
 358
 359         if (changedPoint_[meshPointI])
 360         {
 361             patchInfo.append(allPointInfo_[meshPointI]);
 362
 363             //Pout << "Sending " << meshPointI << " " << mesh_.points()[meshPointI] << " o = " << patchInfo[patchInfo.size()-1].origin() << " " << allPointInfo_[meshPointI] << endl;
 364
 365             patchPoints.append(patchPointI);
 366
 367             label patchFaceI = pointFaces[patchPointI][0];
 368
 369             const face& f = localFaces[patchFaceI];
 370
 371             label index = findIndex(f, patchPointI);
 372
 373             owner.append(patchFaceI);
 374             ownerIndex.append(index);
 375         }
 376     }
 377
 378     patchInfo.shrink();
 379     patchPoints.shrink();
 380     owner.shrink();
 381     ownerIndex.shrink();
 382 }
 383
 384
 385 // Update overall for changed patch points
 386 template <class Type>
 387 void Foam::PointEdgeWave<Type>::updateFromPatchInfo
 388 (
 389     const polyPatch& meshPatch,
 390     const primitivePatch& patch,
 391     const labelList& owner,
 392     const labelList& ownerIndex,
 393     List<Type>& patchInfo
 394 )
 395 {
 396     const faceList& localFaces = patch.localFaces();
 397     const labelList& meshPoints = patch.meshPoints();
 398
 399     // Get patch and mesh points.
 400     labelList changedPatchPoints(patchInfo.size());
 401     labelList changedMeshPoints(patchInfo.size());
 402
 403     forAll(owner, i)
 404     {
 405         label faceI = owner[i];
 406
 407         const face& f = localFaces[faceI];
 408
 409         label index = (f.size() - ownerIndex[i]) % f.size();
 410
 411         changedPatchPoints[i] = f[index];
 412         changedMeshPoints[i] = meshPoints[f[index]];
 413     }
 414
 415     // Adapt for entering domain
 416     enterDomain(meshPatch, patch, changedPatchPoints, patchInfo);
 417
 418     // Merge received info
 419     forAll(patchInfo, i)
 420     {
 421         updatePoint
 422         (
 423             changedMeshPoints[i],
 424             patchInfo[i],
 425             propagationTol_,
 426             allPointInfo_[changedMeshPoints[i]]
 427         );
 428     }
 429 }
 430
 431
 432
 433 // Transfer all the information to/from neighbouring processors
 434 template <class Type>
 435 void Foam::PointEdgeWave<Type>::handleProcPatches()
 436 {
 437     // 1. Send all point info on processor patches. Send as
 438     // face label + offset in face.
 439
 440     forAll(mesh_.boundaryMesh(), patchI)
 441     {
 442         const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 443
 444         if (isA<processorPolyPatch>(patch))
 445         {
 446             // Get all changed points in relative addressing
 447
 448             DynamicList<Type> patchInfo(patch.nPoints());
 449             dynamicLabelList patchPoints(patch.nPoints());
 450             dynamicLabelList owner(patch.nPoints());
 451             dynamicLabelList ownerIndex(patch.nPoints());
 452
 453             getChangedPatchPoints
 454             (
 455                 patch,
 456                 patchInfo,
 457                 patchPoints,
 458                 owner,
 459                 ownerIndex
 460             );
 461
 462             // Adapt for leaving domain
 463             leaveDomain(patch, patch, patchPoints, patchInfo);
 464
 465             const processorPolyPatch& procPatch =
 466                 refCast<const processorPolyPatch>(patch);
 467
 468             if (debug)
 469             {
 470                 Pout<< "Processor patch " << patchI << ' ' << patch.name()
 471                     << " communicating with " << procPatch.neighbProcNo()
 472                     << "  Sending:" << patchInfo.size() << endl;
 473             }
 474
 475             {
 476                 OPstream toNeighbour
 477                 (
 478                     Pstream::blocking,
 479                     procPatch.neighbProcNo()
 480                 );
 481
 482                 toNeighbour << owner << ownerIndex << patchInfo;
 483             }
 484         }
 485     }
 486
 487
 488     //
 489     // 2. Receive all point info on processor patches.
 490     //
 491
 492     forAll(mesh_.boundaryMesh(), patchI)
 493     {
 494         const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 495
 496         if (isA<processorPolyPatch>(patch))
 497         {
 498             const processorPolyPatch& procPatch =
 499                 refCast<const processorPolyPatch>(patch);
 500
 501             List<Type> patchInfo;
 502             labelList owner;
 503             labelList ownerIndex;
 504             {
 505                 IPstream fromNeighbour
 506                 (
 507                     Pstream::blocking,
 508                     procPatch.neighbProcNo()
 509                 );
 510
 511                 fromNeighbour >> owner >> ownerIndex >> patchInfo;
 512             }
 513
 514             if (debug)
 515             {
 516                 Pout<< "Processor patch " << patchI << ' ' << patch.name()
 517                     << " communicating with " << procPatch.neighbProcNo()
 518                     << "  Received:" << patchInfo.size() << endl;
 519             }
 520
 521             // Apply transform to received data for non-parallel planes
 522             if (!procPatch.parallel())
 523             {
 524                 transform(procPatch.forwardT(), patchInfo);
 525             }
 526
 527             updateFromPatchInfo
 528             (
 529                 patch,
 530                 patch,
 531                 owner,
 532                 ownerIndex,
 533                 patchInfo
 534             );
 535         }
 536     }
 537
 538
 539
 540     //
 541     // 3. Handle all shared points
 542     //    (Note:irrespective if changed or not for now)
 543     //
 544
 545     const globalMeshData& pd = mesh_.globalData();
 546
 547     List<Type> sharedData(pd.nGlobalPoints());
 548
 549     forAll(pd.sharedPointLabels(), i)
 550     {
 551         label meshPointI = pd.sharedPointLabels()[i];
 552
 553         // Fill my entries in the shared points
 554         sharedData[pd.sharedPointAddr()[i]] = allPointInfo_[meshPointI];
 555     }
 556
 557     // Combine on master. Reduce operator has to handle a list and call
 558     // Type.updatePoint for all elements
 559     combineReduce(sharedData, listUpdateOp<Type>());
 560
 561     forAll(pd.sharedPointLabels(), i)
 562     {
 563         label meshPointI = pd.sharedPointLabels()[i];
 564
 565         // Retrieve my entries from the shared points
 566         updatePoint
 567         (
 568             meshPointI,
 569             sharedData[pd.sharedPointAddr()[i]],
 570             propagationTol_,
 571             allPointInfo_[meshPointI]
 572         );
 573     }
 574 }
 575
 576
 577 template <class Type>
 578 void Foam::PointEdgeWave<Type>::handleCyclicPatches()
 579 {
 580     // 1. Send all point info on cyclic patches. Send as
 581     // face label + offset in face.
 582
 583     label cycHalf = 0;
 584
 585     forAll(mesh_.boundaryMesh(), patchI)
 586     {
 587         const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 588
 589         if (isA<cyclicPolyPatch>(patch))
 590         {
 591             const primitivePatch& halfA = cycHalves_[cycHalf++];
 592             const primitivePatch& halfB = cycHalves_[cycHalf++];
 593
 594             // HalfA : get all changed points in relative addressing
 595
 596             DynamicList<Type> halfAInfo(halfA.nPoints());
 597             dynamicLabelList halfAPoints(halfA.nPoints());
 598             dynamicLabelList halfAOwner(halfA.nPoints());
 599             dynamicLabelList halfAIndex(halfA.nPoints());
 600
 601             getChangedPatchPoints
 602             (
 603                 halfA,
 604                 halfAInfo,
 605                 halfAPoints,
 606                 halfAOwner,
 607                 halfAIndex
 608             );
 609
 610             // HalfB : get all changed points in relative addressing
 611
 612             DynamicList<Type> halfBInfo(halfB.nPoints());
 613             dynamicLabelList halfBPoints(halfB.nPoints());
 614             dynamicLabelList halfBOwner(halfB.nPoints());
 615             dynamicLabelList halfBIndex(halfB.nPoints());
 616
 617             getChangedPatchPoints
 618             (
 619                 halfB,
 620                 halfBInfo,
 621                 halfBPoints,
 622                 halfBOwner,
 623                 halfBIndex
 624             );
 625
 626
 627             // HalfA : adapt for leaving domain
 628             leaveDomain(patch, halfA, halfAPoints, halfAInfo);
 629
 630             // HalfB : adapt for leaving domain
 631             leaveDomain(patch, halfB, halfBPoints, halfBInfo);
 632
 633
 634             // Apply rotation for non-parallel planes
 635             const cyclicPolyPatch& cycPatch =
 636                 refCast<const cyclicPolyPatch>(patch);
 637
 638             if (!cycPatch.parallel())
 639             {
 640                 // received data from half1
 641                 transform(cycPatch.forwardT(), halfAInfo);
 642
 643                 // received data from half2
 644                 transform(cycPatch.forwardT(), halfBInfo);
 645             }
 646
 647             if (debug)
 648             {
 649                 Pout<< "Cyclic patch " << patchI << ' ' << patch.name()
 650                     << "  Changed on first half : " << halfAInfo.size()
 651                     << "  Changed on second half : " << halfBInfo.size()
 652                     << endl;
 653             }
 654
 655             // Half1: update with data from halfB
 656             updateFromPatchInfo
 657             (
 658                 patch,
 659                 halfA,
 660                 halfBOwner,
 661                 halfBIndex,
 662                 halfBInfo
 663             );
 664
 665             // Half2: update with data from halfA
 666             updateFromPatchInfo
 667             (
 668                 patch,
 669                 halfB,
 670                 halfAOwner,
 671                 halfAIndex,
 672                 halfAInfo
 673             );
 674
 675             if (debug)
 676             {
 677                 //checkCyclic(patch);
 678             }
 679         }
 680     }
 681 }
 682
 683
 684 // Update overall for changed patch points
 685 template <class Type>
 686 void Foam::PointEdgeWave<Type>::updateFromPatchInfo
 687 (
 688     const ggiPolyPatch& to,
 689     const labelListList& shadowAddr,
 690     const labelList& owner,
 691     const labelList& ownerIndex,
 692     List<Type>& patchInfo
 693 )
 694 {
 695     //const labelList& meshPoints = to.meshPoints();
 696     const pointField& points = to.points();
 697     const List<face>& allFaces = mesh_.allFaces();
 698
 699     forAll(patchInfo, i)
 700     {
 701         label fID = to.shadow().zone()[owner[i]];
 702         label pID = allFaces[fID][ownerIndex[i]];
 703         point p = points[pID];
 704
 705         // Update in sending zone without propagation
 706         if(fID >= mesh_.nFaces())
 707         {
 708             allPointInfo_[pID].updatePoint
 709             (
 710                 mesh_,
 711                 pID,
 712                 patchInfo[i],
 713                 propagationTol_
 714             );
 715         }
 716
 717         const labelList& addr = shadowAddr[owner[i]];
 718
 719         if(addr.size() > 0)
 720         {
 721             label nearestPoint = -1;
 722             label nearestFace = -1;
 723             scalar dist = GREAT;
 724
 725             forAll(addr, saI)
 726             {
 727                 label fID2 = to.zone()[addr[saI]];
 728
 729                 const face& f = allFaces[fID2];
 730                 forAll(f, pI)
 731                 {
 732                     label pID2 = f[pI];
 733
 734                     scalar d = magSqr(points[pID2] - p);
 735
 736                     if(d < dist)
 737                     {
 738                         nearestFace = fID2;
 739                         nearestPoint = pID2;
 740                         dist = d;
 741                     }
 742                     else if(nearestPoint == pID2 && fID2 < mesh_.nFaces())
 743                     {
 744                         // Choose face in patch over face in zone
 745                         nearestFace = fID2;
 746                     }
 747                 }
 748             }
 749
 750             patchInfo[i].enterDomain(to, nearestPoint, points[nearestPoint]);
 751
 752             if(nearestFace < mesh_.nFaces())
 753             {
 754                 // Update in receiving patch with propagation
 755                 updatePoint
 756                 (
 757                     nearestPoint,
 758                     patchInfo[i],
 759                     propagationTol_,
 760                     allPointInfo_[nearestPoint]
 761                 );
 762             }
 763             else
 764             {
 765                 // Update in receiving zone without propagation
 766                 allPointInfo_[nearestPoint].updatePoint
 767                 (
 768                     mesh_,
 769                     nearestPoint,
 770                     patchInfo[i],
 771                     propagationTol_
 772                 );
 773             }
 774         }
 775     }
 776 }
 777
 778
 779 // Transfer all the information to/from neighbouring processors
 780 template <class Type>
 781 void Foam::PointEdgeWave<Type>::handleGgiPatches()
 782 {
 783     forAll(mesh_.boundaryMesh(), patchI)
 784     {
 785         const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 786
 787         if (isA<ggiPolyPatch>(patch))
 788         {
 789             const ggiPolyPatch& master = refCast<const ggiPolyPatch>(patch);
 790             const ggiPolyPatch& slave = master.shadow();
 791
 792             if(master.master() && (master.localParallel() || master.size()))
 793             {
 794                 // 1. Collect all point info on master side
 795                 DynamicList<Type> masterInfo(master.nPoints());
 796                 dynamicLabelList masterOwner(master.nPoints());
 797                 dynamicLabelList masterOwnerIndex(master.nPoints());
 798
 799                 {
 800                     dynamicLabelList patchPoints(master.nPoints());
 801
 802                     // Get all changed points in relative addressing
 803                     getChangedPatchPoints
 804                     (
 805                         master,
 806                         masterInfo,
 807                         patchPoints,
 808                         masterOwner,
 809                         masterOwnerIndex
 810                     );
 811
 812                     forAll(masterOwner, i)
 813                     {
 814                         masterOwner[i] =
 815                             master.zoneAddressing()[masterOwner[i]];
 816                     }
 817
 818                     // Adapt for leaving domain
 819                     leaveDomain(master, master, patchPoints, masterInfo);
 820
 821                     if (debug)
 822                     {
 823                         Pout<< "Ggi patch " << master.index() << ' ' << master.name()
 824                             << "  Sending:" << masterInfo.size() << endl;
 825                     }
 826                 }
 827
 828                 // 2. Collect all point info on slave side
 829                 DynamicList<Type> slaveInfo(slave.nPoints());
 830                 dynamicLabelList slaveOwner(slave.nPoints());
 831                 dynamicLabelList slaveOwnerIndex(slave.nPoints());
 832
 833                 {
 834                     dynamicLabelList patchPoints(slave.nPoints());
 835
 836                     // Get all changed points in relative addressing
 837                     getChangedPatchPoints
 838                     (
 839                         slave,
 840                         slaveInfo,
 841                         patchPoints,
 842                         slaveOwner,
 843                         slaveOwnerIndex
 844                     );
 845
 846                     forAll(slaveOwner, i)
 847                     {
 848                         slaveOwner[i] =
 849                             slave.zoneAddressing()[slaveOwner[i]];
 850                     }
 851
 852                     // Adapt for leaving domain
 853                     leaveDomain(slave, slave, patchPoints, slaveInfo);
 854
 855                     if (debug)
 856                     {
 857                         Pout<< "Ggi patch " << slave.index() << ' ' << slave.name()
 858                             << "  Sending:" << slaveInfo.size() << endl;
 859                     }
 860
 861                 }
 862
 863                 if(!master.localParallel())
 864                 {
 865                     combineReduce(masterInfo, listAppendOp<Type>());
 866                     combineReduce(slaveInfo, listAppendOp<Type>());
 867                     combineReduce(masterOwner, listAppendOp<label>());
 868                     combineReduce(slaveOwner, listAppendOp<label>());
 869                     combineReduce(masterOwnerIndex, listAppendOp<label>());
 870                     combineReduce(slaveOwnerIndex, listAppendOp<label>());
 871                 }
 872
 873                 // 3. Apply point info on master & slave side
 874
 875                 if (debug)
 876                 {
 877                     Pout<< "Ggi patch " << slave.index() << ' ' << slave.name()
 878                         << "  Received:" << masterInfo.size() << endl;
 879                     Pout<< "Ggi patch " << master.index() << ' ' << master.name()
 880                         << "  Received:" << slaveInfo.size() << endl;
 881                 }
 882
 883                 // Apply transform to received data for non-parallel planes
 884                 if (!master.parallel())
 885                 {
 886                     transform(master.forwardT(), masterInfo);
 887                     transform(slave.forwardT(), slaveInfo);
 888                 }
 889
 890                 updateFromPatchInfo
 891                 (
 892                     slave,
 893                     master.patchToPatch().masterAddr(),
 894                     masterOwner,
 895                     masterOwnerIndex,
 896                     masterInfo
 897                 );
 898
 899                 updateFromPatchInfo
 900                 (
 901                     master,
 902                     master.patchToPatch().slaveAddr(),
 903                     slaveOwner,
 904                     slaveOwnerIndex,
 905                     slaveInfo
 906                 );
 907             }
 908         }
 909     }
 910 }
 911
 912
 913 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 914
 915 // Iterate, propagating changedPointsInfo across mesh, until no change (or
 916 // maxIter reached). Initial point values specified.
 917 template <class Type>
 918 Foam::PointEdgeWave<Type>::PointEdgeWave
 919 (
 920     const polyMesh& mesh,
 921     const labelList& changedPoints,
 922     const List<Type>& changedPointsInfo,
 923
 924     List<Type>& allPointInfo,
 925     List<Type>& allEdgeInfo,
 926
 927     const label maxIter
 928 )
 929 :
 930     mesh_(mesh),
 931     allPointInfo_(allPointInfo),
 932     allEdgeInfo_(allEdgeInfo),
 933     changedPoint_(mesh_.nPoints(), false),
 934     changedPoints_(mesh_.nPoints()),
 935     nChangedPoints_(0),
 936     changedEdge_(mesh_.nEdges(), false),
 937     changedEdges_(mesh_.nEdges()),
 938     nChangedEdges_(0),
 939     nCyclicPatches_(countPatchType<cyclicPolyPatch>()),
 940     nGgiPatches_(countPatchType<ggiPolyPatch>()),
 941     cycHalves_(2*nCyclicPatches_),
 942     nEvals_(0),
 943     nUnvisitedPoints_(mesh_.nPoints()),
 944     nUnvisitedEdges_(mesh_.nEdges())
 945 {
 946     if
 947     (
 948         allPointInfo_.size() != mesh_.nPoints()
 949         && allPointInfo_.size() != mesh_.allPoints().size()
 950     )
 951     {
 952         FatalErrorIn
 953         (
 954             "PointEdgeWave<Type>::PointEdgeWave"
 955             "(const polyMesh&, const labelList&, const List<Type>,"
 956             " List<Type>&, List<Type>&, const label maxIter)"
 957         )   << "size of pointInfo work array is not equal to the number"
 958             << " of points in the mesh" << endl
 959             << "    pointInfo   :" << allPointInfo_.size() << endl
 960             << "    mesh.nPoints:" << mesh_.nPoints()
 961             << exit(FatalError);
 962     }
 963     if (allEdgeInfo_.size() != mesh_.nEdges())
 964     {
 965         FatalErrorIn
 966         (
 967             "PointEdgeWave<Type>::PointEdgeWave"
 968             "(const polyMesh&, const labelList&, const List<Type>,"
 969             " List<Type>&, List<Type>&, const label maxIter)"
 970         )   << "size of edgeInfo work array is not equal to the number"
 971             << " of edges in the mesh" << endl
 972             << "    edgeInfo   :" << allEdgeInfo_.size() << endl
 973             << "    mesh.nEdges:" << mesh_.nEdges()
 974             << exit(FatalError);
 975     }
 976
 977
 978     // Calculate cyclic halves addressing.
 979     if (nCyclicPatches_ > 0)
 980     {
 981         calcCyclicAddressing();
 982     }
 983
 984
 985     // Set from initial changed points data
 986     setPointInfo(changedPoints, changedPointsInfo);
 987
 988     if (debug)
 989     {
 990         Pout<< "Seed points               : " << nChangedPoints_ << endl;
 991     }
 992
 993     // Iterate until nothing changes
 994     label iter = iterate(maxIter);
 995
 996     if ((maxIter > 0) && (iter >= maxIter))
 997     {
 998         FatalErrorIn
 999         (
1000             "PointEdgeWave<Type>::PointEdgeWave"
1001             "(const polyMesh&, const labelList&, const List<Type>,"
1002             " List<Type>&, List<Type>&, const label maxIter)"
1003         )   << "Maximum number of iterations reached. Increase maxIter." << nl
1004             << "    maxIter:" << maxIter << endl
1005             << "    nChangedPoints:" << nChangedPoints_ << endl
1006             << "    nChangedEdges:" << nChangedEdges_ << endl
1007             << exit(FatalError);
1008     }
1009 }
1010
1011
1012 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
1013
1014 template <class Type>
1015 Foam::PointEdgeWave<Type>::~PointEdgeWave()
1016 {}
1017
1018
1019 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1020
1021
1022 template <class Type>
1023 Foam::label Foam::PointEdgeWave<Type>::getUnsetPoints() const
1024 {
1025     return nUnvisitedPoints_;
1026 }
1027
1028
1029 template <class Type>
1030 Foam::label Foam::PointEdgeWave<Type>::getUnsetEdges() const
1031 {
1032     return nUnvisitedEdges_;
1033 }
1034
1035
1036 // Copy point information into member data
1037 template <class Type>
1038 void Foam::PointEdgeWave<Type>::setPointInfo
1039 (
1040     const labelList& changedPoints,
1041     const List<Type>& changedPointsInfo
1042 )
1043 {
1044     forAll(changedPoints, changedPointI)
1045     {
1046         label pointI = changedPoints[changedPointI];
1047
1048         bool wasValid = allPointInfo_[pointI].valid();
1049
1050         // Copy info for pointI
1051         allPointInfo_[pointI] = changedPointsInfo[changedPointI];
1052
1053         // Maintain count of unset points
1054         if (!wasValid && allPointInfo_[pointI].valid())
1055         {
1056             --nUnvisitedPoints_;
1057         }
1058
1059         // Mark pointI as changed, both on list and on point itself.
1060
1061         if (!changedPoint_[pointI])
1062         {
1063             changedPoint_[pointI] = true;
1064             changedPoints_[nChangedPoints_++] = pointI;
1065         }
1066     }
1067 }
1068
1069
1070 // Propagate information from edge to point. Return number of points changed.
1071 template <class Type>
1072 Foam::label Foam::PointEdgeWave<Type>::edgeToPoint()
1073 {
1074     for
1075     (
1076         label changedEdgeI = 0;
1077         changedEdgeI < nChangedEdges_;
1078         changedEdgeI++
1079     )
1080     {
1081         label edgeI = changedEdges_[changedEdgeI];
1082
1083         if (!changedEdge_[edgeI])
1084         {
1085             FatalErrorIn("PointEdgeWave<Type>::edgeToPoint()")
1086                 << "edge " << edgeI
1087                 << " not marked as having been changed" << nl
1088                 << "This might be caused by multiple occurences of the same"
1089                 << " seed point." << abort(FatalError);
1090         }
1091
1092
1093         const Type& neighbourWallInfo = allEdgeInfo_[edgeI];
1094
1095         // Evaluate all connected points (= edge endpoints)
1096         const edge& e = mesh_.edges()[edgeI];
1097
1098         forAll(e, eI)
1099         {
1100             Type& currentWallInfo = allPointInfo_[e[eI]];
1101
1102             if (currentWallInfo != neighbourWallInfo)
1103             {
1104                 updatePoint
1105                 (
1106                     e[eI],
1107                     edgeI,
1108                     neighbourWallInfo,
1109                     propagationTol_,
1110                     currentWallInfo
1111                 );
1112             }
1113         }
1114
1115         // Reset status of edge
1116         changedEdge_[edgeI] = false;
1117     }
1118
1119     // Handled all changed edges by now
1120     nChangedEdges_ = 0;
1121
1122     if (nCyclicPatches_ > 0)
1123     {
1124         // Transfer changed points across cyclic halves
1125         handleCyclicPatches();
1126     }
1127     if (nGgiPatches_ > 0)
1128     {
1129         // Transfer changed points across cyclic halves
1130         handleGgiPatches();
1131     }
1132     if (Pstream::parRun())
1133     {
1134         // Transfer changed points from neighbouring processors.
1135         handleProcPatches();
1136     }
1137
1138     if (debug)
1139     {
1140         Pout<< "Changed points            : " << nChangedPoints_ << endl;
1141     }
1142
1143     // Sum nChangedPoints over all procs
1144     label totNChanged = nChangedPoints_;
1145
1146     reduce(totNChanged, sumOp<label>());
1147
1148     return totNChanged;
1149 }
1150
1151
1152 // Propagate information from point to edge. Return number of edges changed.
1153 template <class Type>
1154 Foam::label Foam::PointEdgeWave<Type>::pointToEdge()
1155 {
1156     const labelListList& pointEdges = mesh_.pointEdges();
1157
1158     for
1159     (
1160         label changedPointI = 0;
1161         changedPointI < nChangedPoints_;
1162         changedPointI++
1163     )
1164     {
1165         label pointI = changedPoints_[changedPointI];
1166
1167         if (!changedPoint_[pointI])
1168         {
1169             FatalErrorIn("PointEdgeWave<Type>::pointToEdge()")
1170                 << "Point " << pointI
1171                 << " not marked as having been changed" << nl
1172                 << "This might be caused by multiple occurences of the same"
1173                 << " seed point." << abort(FatalError);
1174         }
1175
1176         const Type& neighbourWallInfo = allPointInfo_[pointI];
1177
1178         // Evaluate all connected edges
1179
1180         const labelList& edgeLabels = pointEdges[pointI];
1181         forAll(edgeLabels, edgeLabelI)
1182         {
1183             label edgeI = edgeLabels[edgeLabelI];
1184
1185             Type& currentWallInfo = allEdgeInfo_[edgeI];
1186
1187             if (currentWallInfo != neighbourWallInfo)
1188             {
1189                 updateEdge
1190                 (
1191                     edgeI,
1192                     pointI,
1193                     neighbourWallInfo,
1194                     propagationTol_,
1195                     currentWallInfo
1196                 );
1197             }
1198         }
1199
1200         // Reset status of point
1201         changedPoint_[pointI] = false;
1202     }
1203
1204     // Handled all changed points by now
1205     nChangedPoints_ = 0;
1206
1207     if (debug)
1208     {
1209         Pout<< "Changed edges             : " << nChangedEdges_ << endl;
1210     }
1211
1212     // Sum nChangedPoints over all procs
1213     label totNChanged = nChangedEdges_;
1214
1215     reduce(totNChanged, sumOp<label>());
1216
1217     return totNChanged;
1218 }
1219
1220
1221 // Iterate
1222 template <class Type>
1223 Foam::label Foam::PointEdgeWave<Type>::iterate(const label maxIter)
1224 {
1225     if (nCyclicPatches_ > 0)
1226     {
1227         // Transfer changed points across cyclic halves
1228         handleCyclicPatches();
1229     }
1230     if (nGgiPatches_ > 0)
1231     {
1232         // Transfer changed points across ggi patches
1233         handleGgiPatches();
1234     }
1235     if (Pstream::parRun())
1236     {
1237         // Transfer changed points from neighbouring processors.
1238         handleProcPatches();
1239     }
1240
1241     nEvals_ = 0;
1242
1243     label iter = 0;
1244
1245     while(iter < maxIter)
1246     {
1247         if (debug)
1248         {
1249             Pout<< "Iteration " << iter << endl;
1250         }
1251
1252         label nEdges = pointToEdge();
1253
1254         if (debug)
1255         {
1256             Pout<< "Total changed edges       : " << nEdges << endl;
1257         }
1258
1259         if (nEdges == 0)
1260         {
1261             break;
1262         }
1263
1264         label nPoints = edgeToPoint();
1265
1266         if (debug)
1267         {
1268             Pout<< "Total changed points      : " << nPoints << endl;
1269
1270             Pout<< "Total evaluations         : " << nEvals_ << endl;
1271
1272             Pout<< "Remaining unvisited points: " << nUnvisitedPoints_ << endl;
1273
1274             Pout<< "Remaining unvisited edges : " << nUnvisitedEdges_ << endl;
1275
1276             Pout<< endl;
1277         }
1278
1279         if (nPoints == 0)
1280         {
1281             break;
1282         }
1283
1284         iter++;
1285     }
1286
1287     return iter;
1288 }
1289
1290
1291 // ************************************************************************* //