src/meshTools/PointEdgeWave/PointEdgeWave.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2010 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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, 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 "OPstream.H"
  31 #include "IPstream.H"
  32 #include "PstreamCombineReduceOps.H"
  33 #include "debug.H"
  34 #include "typeInfo.H"
  35 #include "globalMeshData.H"
  36 #include "pointFields.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 template <class Type, class TrackingData>
  41 Foam::scalar Foam::PointEdgeWave<Type, TrackingData>::propagationTol_ = 0.01;
  42
  43 template <class Type, class TrackingData>
  44 Foam::label Foam::PointEdgeWave<Type, TrackingData>::dummyTrackData_ = 12345;
  45
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 // Handle leaving domain. Implementation referred to Type
  50 template <class Type, class TrackingData>
  51 void Foam::PointEdgeWave<Type, TrackingData>::leaveDomain
  52 (
  53     const polyPatch& patch,
  54     const labelList& patchPointLabels,
  55     List<Type>& pointInfo
  56 ) const
  57 {
  58     const labelList& meshPoints = patch.meshPoints();
  59
  60     forAll(patchPointLabels, i)
  61     {
  62         label patchPointI = patchPointLabels[i];
  63
  64         const point& pt = patch.points()[meshPoints[patchPointI]];
  65
  66         pointInfo[i].leaveDomain(patch, patchPointI, pt, td_);
  67     }
  68 }
  69
  70
  71 // Handle entering domain. Implementation referred to Type
  72 template <class Type, class TrackingData>
  73 void Foam::PointEdgeWave<Type, TrackingData>::enterDomain
  74 (
  75     const polyPatch& patch,
  76     const labelList& patchPointLabels,
  77     List<Type>& pointInfo
  78 ) const
  79 {
  80     const labelList& meshPoints = patch.meshPoints();
  81
  82     forAll(patchPointLabels, i)
  83     {
  84         label patchPointI = patchPointLabels[i];
  85
  86         const point& pt = patch.points()[meshPoints[patchPointI]];
  87
  88         pointInfo[i].enterDomain(patch, patchPointI, pt, td_);
  89     }
  90 }
  91
  92
  93 // Transform. Implementation referred to Type
  94 template <class Type, class TrackingData>
  95 void Foam::PointEdgeWave<Type, TrackingData>::transform
  96 (
  97     const polyPatch& patch,
  98     const tensorField& rotTensor,
  99     List<Type>& pointInfo
 100 ) const
 101 {
 102     if (rotTensor.size() == 1)
 103     {
 104         const tensor& T = rotTensor[0];
 105
 106         forAll(pointInfo, i)
 107         {
 108             pointInfo[i].transform(T, td_);
 109         }
 110     }
 111     else
 112     {
 113         FatalErrorIn
 114         (
 115             "PointEdgeWave<Type, TrackingData>::transform"
 116             "(const tensorField&, List<Type>&)"
 117         )   << "Non-uniform transformation on patch " << patch.name()
 118             << " not supported for point fields"
 119             << abort(FatalError);
 120
 121         forAll(pointInfo, i)
 122         {
 123             pointInfo[i].transform(rotTensor[i], td_);
 124         }
 125     }
 126 }
 127
 128
 129 // Update info for pointI, at position pt, with information from
 130 // neighbouring edge.
 131 // Updates:
 132 //      - changedPoint_, changedPoints_, nChangedPoints_,
 133 //      - statistics: nEvals_, nUnvisitedPoints_
 134 template <class Type, class TrackingData>
 135 bool Foam::PointEdgeWave<Type, TrackingData>::updatePoint
 136 (
 137     const label pointI,
 138     const label neighbourEdgeI,
 139     const Type& neighbourInfo,
 140     Type& pointInfo
 141 )
 142 {
 143     nEvals_++;
 144
 145     bool wasValid = pointInfo.valid(td_);
 146
 147     bool propagate =
 148         pointInfo.updatePoint
 149         (
 150             mesh_,
 151             pointI,
 152             neighbourEdgeI,
 153             neighbourInfo,
 154             propagationTol_,
 155             td_
 156         );
 157
 158     if (propagate)
 159     {
 160         if (!changedPoint_[pointI])
 161         {
 162             changedPoint_[pointI] = true;
 163             changedPoints_[nChangedPoints_++] = pointI;
 164         }
 165     }
 166
 167     if (!wasValid && pointInfo.valid(td_))
 168     {
 169         --nUnvisitedPoints_;
 170     }
 171
 172     return propagate;
 173 }
 174
 175
 176 // Update info for pointI, at position pt, with information from
 177 // same point.
 178 // Updates:
 179 //      - changedPoint_, changedPoints_, nChangedPoints_,
 180 //      - statistics: nEvals_, nUnvisitedPoints_
 181 template <class Type, class TrackingData>
 182 bool Foam::PointEdgeWave<Type, TrackingData>::updatePoint
 183 (
 184     const label pointI,
 185     const Type& neighbourInfo,
 186     Type& pointInfo
 187 )
 188 {
 189     nEvals_++;
 190
 191     bool wasValid = pointInfo.valid(td_);
 192
 193     bool propagate =
 194         pointInfo.updatePoint
 195         (
 196             mesh_,
 197             pointI,
 198             neighbourInfo,
 199             propagationTol_,
 200             td_
 201         );
 202
 203     if (propagate)
 204     {
 205         if (!changedPoint_[pointI])
 206         {
 207             changedPoint_[pointI] = true;
 208             changedPoints_[nChangedPoints_++] = pointI;
 209         }
 210     }
 211
 212     if (!wasValid && pointInfo.valid(td_))
 213     {
 214         --nUnvisitedPoints_;
 215     }
 216
 217     return propagate;
 218 }
 219
 220
 221 // Update info for edgeI, at position pt, with information from
 222 // neighbouring point.
 223 // Updates:
 224 //      - changedEdge_, changedEdges_, nChangedEdges_,
 225 //      - statistics: nEvals_, nUnvisitedEdge_
 226 template <class Type, class TrackingData>
 227 bool Foam::PointEdgeWave<Type, TrackingData>::updateEdge
 228 (
 229     const label edgeI,
 230     const label neighbourPointI,
 231     const Type& neighbourInfo,
 232     Type& edgeInfo
 233 )
 234 {
 235     nEvals_++;
 236
 237     bool wasValid = edgeInfo.valid(td_);
 238
 239     bool propagate =
 240         edgeInfo.updateEdge
 241         (
 242             mesh_,
 243             edgeI,
 244             neighbourPointI,
 245             neighbourInfo,
 246             propagationTol_,
 247             td_
 248         );
 249
 250     if (propagate)
 251     {
 252         if (!changedEdge_[edgeI])
 253         {
 254             changedEdge_[edgeI] = true;
 255             changedEdges_[nChangedEdges_++] = edgeI;
 256         }
 257     }
 258
 259     if (!wasValid && edgeInfo.valid(td_))
 260     {
 261         --nUnvisitedEdges_;
 262     }
 263
 264     return propagate;
 265 }
 266
 267
 268 // Check if patches of given type name are present
 269 template <class Type, class TrackingData>
 270 template <class PatchType>
 271 Foam::label Foam::PointEdgeWave<Type, TrackingData>::countPatchType() const
 272 {
 273     label nPatches = 0;
 274
 275     forAll(mesh_.boundaryMesh(), patchI)
 276     {
 277         if (isA<PatchType>(mesh_.boundaryMesh()[patchI]))
 278         {
 279             nPatches++;
 280         }
 281     }
 282     return nPatches;
 283 }
 284
 285
 286 // Transfer all the information to/from neighbouring processors
 287 template <class Type, class TrackingData>
 288 void Foam::PointEdgeWave<Type, TrackingData>::handleProcPatches()
 289 {
 290     // 1. Send all point info on processor patches.
 291
 292     PstreamBuffers pBufs(Pstream::nonBlocking);
 293
 294     DynamicList<Type> patchInfo;
 295     DynamicList<label> thisPoints;
 296     DynamicList<label> nbrPoints;
 297
 298     forAll(mesh_.globalData().processorPatches(), i)
 299     {
 300         label patchI = mesh_.globalData().processorPatches()[i];
 301         const processorPolyPatch& procPatch =
 302             refCast<const processorPolyPatch>(mesh_.boundaryMesh()[patchI]);
 303
 304         patchInfo.clear();
 305         patchInfo.reserve(procPatch.nPoints());
 306         thisPoints.clear();
 307         thisPoints.reserve(procPatch.nPoints());
 308         nbrPoints.clear();
 309         nbrPoints.reserve(procPatch.nPoints());
 310
 311         // Get all changed points in reverse order
 312         const labelList& neighbPoints = procPatch.neighbPoints();
 313         forAll(neighbPoints, thisPointI)
 314         {
 315             label meshPointI = procPatch.meshPoints()[thisPointI];
 316             if (changedPoint_[meshPointI])
 317             {
 318                 patchInfo.append(allPointInfo_[meshPointI]);
 319                 thisPoints.append(thisPointI);
 320                 nbrPoints.append(neighbPoints[thisPointI]);
 321             }
 322         }
 323
 324         // Adapt for leaving domain
 325         leaveDomain(procPatch, thisPoints, patchInfo);
 326
 327         if (debug)
 328         {
 329             Pout<< "Processor patch " << patchI << ' ' << procPatch.name()
 330                 << " communicating with " << procPatch.neighbProcNo()
 331                 << "  Sending:" << patchInfo.size() << endl;
 332         }
 333
 334         UOPstream toNeighbour(procPatch.neighbProcNo(), pBufs);
 335         toNeighbour << nbrPoints << patchInfo;
 336     }
 337
 338
 339     pBufs.finishedSends();
 340
 341     //
 342     // 2. Receive all point info on processor patches.
 343     //
 344
 345     forAll(mesh_.globalData().processorPatches(), i)
 346     {
 347         label patchI = mesh_.globalData().processorPatches()[i];
 348         const processorPolyPatch& procPatch =
 349             refCast<const processorPolyPatch>(mesh_.boundaryMesh()[patchI]);
 350
 351         List<Type> patchInfo;
 352         labelList patchPoints;
 353
 354         {
 355             UIPstream fromNeighbour(procPatch.neighbProcNo(), pBufs);
 356             fromNeighbour >> patchPoints >> patchInfo;
 357         }
 358
 359         if (debug)
 360         {
 361             Pout<< "Processor patch " << patchI << ' ' << procPatch.name()
 362                 << " communicating with " << procPatch.neighbProcNo()
 363                 << "  Received:" << patchInfo.size() << endl;
 364         }
 365
 366         // Apply transform to received data for non-parallel planes
 367         if (!procPatch.parallel())
 368         {
 369             transform(procPatch, procPatch.forwardT(), patchInfo);
 370         }
 371
 372         // Adapt for entering domain
 373         enterDomain(procPatch, patchPoints, patchInfo);
 374
 375         // Merge received info
 376         const labelList& meshPoints = procPatch.meshPoints();
 377         forAll(patchInfo, i)
 378         {
 379             label meshPointI = meshPoints[patchPoints[i]];
 380
 381             if (!allPointInfo_[meshPointI].equal(patchInfo[i], td_))
 382             {
 383                 updatePoint
 384                 (
 385                     meshPointI,
 386                     patchInfo[i],
 387                     allPointInfo_[meshPointI]
 388                 );
 389             }
 390         }
 391     }
 392
 393
 394
 395     //
 396     // 3. Handle all shared points
 397     //    (Note:irrespective if changed or not for now)
 398     //
 399
 400     const globalMeshData& pd = mesh_.globalData();
 401
 402     List<Type> sharedData(pd.nGlobalPoints());
 403
 404     forAll(pd.sharedPointLabels(), i)
 405     {
 406         label meshPointI = pd.sharedPointLabels()[i];
 407
 408         // Fill my entries in the shared points
 409         sharedData[pd.sharedPointAddr()[i]] = allPointInfo_[meshPointI];
 410     }
 411
 412     // Combine on master. Reduce operator has to handle a list and call
 413     // Type.updatePoint for all elements
 414     combineReduce(sharedData, listUpdateOp<Type>(propagationTol_, td_));
 415
 416     forAll(pd.sharedPointLabels(), i)
 417     {
 418         label meshPointI = pd.sharedPointLabels()[i];
 419
 420         // Retrieve my entries from the shared points.
 421         const Type& nbrInfo = sharedData[pd.sharedPointAddr()[i]];
 422
 423         if (!allPointInfo_[meshPointI].equal(nbrInfo, td_))
 424         {
 425             updatePoint
 426             (
 427                 meshPointI,
 428                 nbrInfo,
 429                 allPointInfo_[meshPointI]
 430             );
 431         }
 432     }
 433 }
 434
 435
 436 template <class Type, class TrackingData>
 437 void Foam::PointEdgeWave<Type, TrackingData>::handleCyclicPatches()
 438 {
 439     // 1. Send all point info on cyclic patches.
 440
 441     DynamicList<Type> nbrInfo;
 442     DynamicList<label> nbrPoints;
 443     DynamicList<label> thisPoints;
 444
 445     forAll(mesh_.boundaryMesh(), patchI)
 446     {
 447         const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 448
 449         if (isA<cyclicPolyPatch>(patch))
 450         {
 451             const cyclicPolyPatch& cycPatch =
 452                 refCast<const cyclicPolyPatch>(patch);
 453
 454             nbrInfo.clear();
 455             nbrInfo.reserve(cycPatch.nPoints());
 456             nbrPoints.clear();
 457             nbrPoints.reserve(cycPatch.nPoints());
 458             thisPoints.clear();
 459             thisPoints.reserve(cycPatch.nPoints());
 460
 461             // Collect nbrPatch points that have changed
 462             {
 463                 const cyclicPolyPatch& nbrPatch = cycPatch.neighbPatch();
 464                 const edgeList& pairs = cycPatch.coupledPoints();
 465                 const labelList& meshPoints = nbrPatch.meshPoints();
 466
 467                 forAll(pairs, pairI)
 468                 {
 469                     label thisPointI = pairs[pairI][0];
 470                     label nbrPointI = pairs[pairI][1];
 471                     label meshPointI = meshPoints[nbrPointI];
 472
 473                     if (changedPoint_[meshPointI])
 474                     {
 475                         nbrInfo.append(allPointInfo_[meshPointI]);
 476                         nbrPoints.append(nbrPointI);
 477                         thisPoints.append(thisPointI);
 478                     }
 479                 }
 480
 481                 // nbr : adapt for leaving domain
 482                 leaveDomain(nbrPatch, nbrPoints, nbrInfo);
 483             }
 484
 485
 486             // Apply rotation for non-parallel planes
 487
 488             if (!cycPatch.parallel())
 489             {
 490                 // received data from half1
 491                 transform(cycPatch, cycPatch.forwardT(), nbrInfo);
 492             }
 493
 494             if (debug)
 495             {
 496                 Pout<< "Cyclic patch " << patchI << ' ' << patch.name()
 497                     << "  Changed : " << nbrInfo.size()
 498                     << endl;
 499             }
 500
 501             // Adapt for entering domain
 502             enterDomain(cycPatch, thisPoints, nbrInfo);
 503
 504             // Merge received info
 505             const labelList& meshPoints = cycPatch.meshPoints();
 506             forAll(nbrInfo, i)
 507             {
 508                 label meshPointI = meshPoints[thisPoints[i]];
 509
 510                 if (!allPointInfo_[meshPointI].equal(nbrInfo[i], td_))
 511                 {
 512                     updatePoint
 513                     (
 514                         meshPointI,
 515                         nbrInfo[i],
 516                         allPointInfo_[meshPointI]
 517                     );
 518                 }
 519             }
 520         }
 521     }
 522 }
 523
 524
 525 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 526
 527 // Iterate, propagating changedPointsInfo across mesh, until no change (or
 528 // maxIter reached). Initial point values specified.
 529 template <class Type, class TrackingData>
 530 Foam::PointEdgeWave<Type, TrackingData>::PointEdgeWave
 531 (
 532     const polyMesh& mesh,
 533     const labelList& changedPoints,
 534     const List<Type>& changedPointsInfo,
 535
 536     UList<Type>& allPointInfo,
 537     UList<Type>& allEdgeInfo,
 538
 539     const label maxIter,
 540     TrackingData& td
 541 )
 542 :
 543     mesh_(mesh),
 544     allPointInfo_(allPointInfo),
 545     allEdgeInfo_(allEdgeInfo),
 546     td_(td),
 547     changedPoint_(mesh_.nPoints(), false),
 548     changedPoints_(mesh_.nPoints()),
 549     nChangedPoints_(0),
 550     changedEdge_(mesh_.nEdges(), false),
 551     changedEdges_(mesh_.nEdges()),
 552     nChangedEdges_(0),
 553     nCyclicPatches_(countPatchType<cyclicPolyPatch>()),
 554     nEvals_(0),
 555     nUnvisitedPoints_(mesh_.nPoints()),
 556     nUnvisitedEdges_(mesh_.nEdges())
 557 {
 558     if (allPointInfo_.size() != mesh_.nPoints())
 559     {
 560         FatalErrorIn
 561         (
 562             "PointEdgeWave<Type, TrackingData>::PointEdgeWave"
 563             "(const polyMesh&, const labelList&, const List<Type>,"
 564             " List<Type>&, List<Type>&, const label maxIter)"
 565         )   << "size of pointInfo work array is not equal to the number"
 566             << " of points in the mesh" << endl
 567             << "    pointInfo   :" << allPointInfo_.size() << endl
 568             << "    mesh.nPoints:" << mesh_.nPoints()
 569             << exit(FatalError);
 570     }
 571     if (allEdgeInfo_.size() != mesh_.nEdges())
 572     {
 573         FatalErrorIn
 574         (
 575             "PointEdgeWave<Type, TrackingData>::PointEdgeWave"
 576             "(const polyMesh&, const labelList&, const List<Type>,"
 577             " List<Type>&, List<Type>&, const label maxIter)"
 578         )   << "size of edgeInfo work array is not equal to the number"
 579             << " of edges in the mesh" << endl
 580             << "    edgeInfo   :" << allEdgeInfo_.size() << endl
 581             << "    mesh.nEdges:" << mesh_.nEdges()
 582             << exit(FatalError);
 583     }
 584
 585
 586     // Set from initial changed points data
 587     setPointInfo(changedPoints, changedPointsInfo);
 588
 589     if (debug)
 590     {
 591         Pout<< "Seed points               : " << nChangedPoints_ << endl;
 592     }
 593
 594     // Iterate until nothing changes
 595     label iter = iterate(maxIter);
 596
 597     if ((maxIter > 0) && (iter >= maxIter))
 598     {
 599         FatalErrorIn
 600         (
 601             "PointEdgeWave<Type, TrackingData>::PointEdgeWave"
 602             "(const polyMesh&, const labelList&, const List<Type>,"
 603             " List<Type>&, List<Type>&, const label maxIter)"
 604         )   << "Maximum number of iterations reached. Increase maxIter." << endl
 605             << "    maxIter:" << maxIter << endl
 606             << "    nChangedPoints:" << nChangedPoints_ << endl
 607             << "    nChangedEdges:" << nChangedEdges_ << endl
 608             << exit(FatalError);
 609     }
 610 }
 611
 612
 613 template <class Type, class TrackingData>
 614 Foam::PointEdgeWave<Type, TrackingData>::PointEdgeWave
 615 (
 616     const polyMesh& mesh,
 617     UList<Type>& allPointInfo,
 618     UList<Type>& allEdgeInfo,
 619     TrackingData& td
 620 )
 621 :
 622     mesh_(mesh),
 623     allPointInfo_(allPointInfo),
 624     allEdgeInfo_(allEdgeInfo),
 625     td_(td),
 626     changedPoint_(mesh_.nPoints(), false),
 627     changedPoints_(mesh_.nPoints()),
 628     nChangedPoints_(0),
 629     changedEdge_(mesh_.nEdges(), false),
 630     changedEdges_(mesh_.nEdges()),
 631     nChangedEdges_(0),
 632     nCyclicPatches_(countPatchType<cyclicPolyPatch>()),
 633     nEvals_(0),
 634     nUnvisitedPoints_(mesh_.nPoints()),
 635     nUnvisitedEdges_(mesh_.nEdges())
 636 {}
 637
 638
 639 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 640
 641 template <class Type, class TrackingData>
 642 Foam::PointEdgeWave<Type, TrackingData>::~PointEdgeWave()
 643 {}
 644
 645
 646 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 647
 648
 649 template <class Type, class TrackingData>
 650 Foam::label Foam::PointEdgeWave<Type, TrackingData>::getUnsetPoints() const
 651 {
 652     return nUnvisitedPoints_;
 653 }
 654
 655
 656 template <class Type, class TrackingData>
 657 Foam::label Foam::PointEdgeWave<Type, TrackingData>::getUnsetEdges() const
 658 {
 659     return nUnvisitedEdges_;
 660 }
 661
 662
 663 // Copy point information into member data
 664 template <class Type, class TrackingData>
 665 void Foam::PointEdgeWave<Type, TrackingData>::setPointInfo
 666 (
 667     const labelList& changedPoints,
 668     const List<Type>& changedPointsInfo
 669 )
 670 {
 671     forAll(changedPoints, changedPointI)
 672     {
 673         label pointI = changedPoints[changedPointI];
 674
 675         bool wasValid = allPointInfo_[pointI].valid(td_);
 676
 677         // Copy info for pointI
 678         allPointInfo_[pointI] = changedPointsInfo[changedPointI];
 679
 680         // Maintain count of unset points
 681         if (!wasValid && allPointInfo_[pointI].valid(td_))
 682         {
 683             --nUnvisitedPoints_;
 684         }
 685
 686         // Mark pointI as changed, both on list and on point itself.
 687
 688         if (!changedPoint_[pointI])
 689         {
 690             changedPoint_[pointI] = true;
 691             changedPoints_[nChangedPoints_++] = pointI;
 692         }
 693     }
 694 }
 695
 696
 697 // Propagate information from edge to point. Return number of points changed.
 698 template <class Type, class TrackingData>
 699 Foam::label Foam::PointEdgeWave<Type, TrackingData>::edgeToPoint()
 700 {
 701     for
 702     (
 703         label changedEdgeI = 0;
 704         changedEdgeI < nChangedEdges_;
 705         changedEdgeI++
 706     )
 707     {
 708         label edgeI = changedEdges_[changedEdgeI];
 709
 710         if (!changedEdge_[edgeI])
 711         {
 712             FatalErrorIn("PointEdgeWave<Type, TrackingData>::edgeToPoint()")
 713                 << "edge " << edgeI
 714                 << " not marked as having been changed" << nl
 715                 << "This might be caused by multiple occurences of the same"
 716                 << " seed point." << abort(FatalError);
 717         }
 718
 719
 720         const Type& neighbourWallInfo = allEdgeInfo_[edgeI];
 721
 722         // Evaluate all connected points (= edge endpoints)
 723         const edge& e = mesh_.edges()[edgeI];
 724
 725         forAll(e, eI)
 726         {
 727             Type& currentWallInfo = allPointInfo_[e[eI]];
 728
 729             if (!currentWallInfo.equal(neighbourWallInfo, td_))
 730             {
 731                 updatePoint
 732                 (
 733                     e[eI],
 734                     edgeI,
 735                     neighbourWallInfo,
 736                     currentWallInfo
 737                 );
 738             }
 739         }
 740
 741         // Reset status of edge
 742         changedEdge_[edgeI] = false;
 743     }
 744
 745     // Handled all changed edges by now
 746     nChangedEdges_ = 0;
 747
 748     if (nCyclicPatches_ > 0)
 749     {
 750         // Transfer changed points across cyclic halves
 751         handleCyclicPatches();
 752     }
 753     if (Pstream::parRun())
 754     {
 755         // Transfer changed points from neighbouring processors.
 756         handleProcPatches();
 757     }
 758
 759     if (debug)
 760     {
 761         Pout<< "Changed points            : " << nChangedPoints_ << endl;
 762     }
 763
 764     // Sum nChangedPoints over all procs
 765     label totNChanged = nChangedPoints_;
 766
 767     reduce(totNChanged, sumOp<label>());
 768
 769     return totNChanged;
 770 }
 771
 772
 773 // Propagate information from point to edge. Return number of edges changed.
 774 template <class Type, class TrackingData>
 775 Foam::label Foam::PointEdgeWave<Type, TrackingData>::pointToEdge()
 776 {
 777     const labelListList& pointEdges = mesh_.pointEdges();
 778
 779     for
 780     (
 781         label changedPointI = 0;
 782         changedPointI < nChangedPoints_;
 783         changedPointI++
 784     )
 785     {
 786         label pointI = changedPoints_[changedPointI];
 787
 788         if (!changedPoint_[pointI])
 789         {
 790             FatalErrorIn("PointEdgeWave<Type, TrackingData>::pointToEdge()")
 791                 << "Point " << pointI
 792                 << " not marked as having been changed" << nl
 793                 << "This might be caused by multiple occurences of the same"
 794                 << " seed point." << abort(FatalError);
 795         }
 796
 797         const Type& neighbourWallInfo = allPointInfo_[pointI];
 798
 799         // Evaluate all connected edges
 800
 801         const labelList& edgeLabels = pointEdges[pointI];
 802         forAll(edgeLabels, edgeLabelI)
 803         {
 804             label edgeI = edgeLabels[edgeLabelI];
 805
 806             Type& currentWallInfo = allEdgeInfo_[edgeI];
 807
 808             if (!currentWallInfo.equal(neighbourWallInfo, td_))
 809             {
 810                 updateEdge
 811                 (
 812                     edgeI,
 813                     pointI,
 814                     neighbourWallInfo,
 815                     currentWallInfo
 816                 );
 817             }
 818         }
 819
 820         // Reset status of point
 821         changedPoint_[pointI] = false;
 822     }
 823
 824     // Handled all changed points by now
 825     nChangedPoints_ = 0;
 826
 827     if (debug)
 828     {
 829         Pout<< "Changed edges             : " << nChangedEdges_ << endl;
 830     }
 831
 832     // Sum nChangedPoints over all procs
 833     label totNChanged = nChangedEdges_;
 834
 835     reduce(totNChanged, sumOp<label>());
 836
 837     return totNChanged;
 838 }
 839
 840
 841 // Iterate
 842 template <class Type, class TrackingData>
 843 Foam::label Foam::PointEdgeWave<Type, TrackingData>::iterate
 844 (
 845     const label maxIter
 846 )
 847 {
 848     if (nCyclicPatches_ > 0)
 849     {
 850         // Transfer changed points across cyclic halves
 851         handleCyclicPatches();
 852     }
 853     if (Pstream::parRun())
 854     {
 855         // Transfer changed points from neighbouring processors.
 856         handleProcPatches();
 857     }
 858
 859     nEvals_ = 0;
 860
 861     label iter = 0;
 862
 863     while (iter < maxIter)
 864     {
 865         if (debug)
 866         {
 867             Pout<< "Iteration " << iter << endl;
 868         }
 869
 870         label nEdges = pointToEdge();
 871
 872         if (debug)
 873         {
 874             Pout<< "Total changed edges       : " << nEdges << endl;
 875         }
 876
 877         if (nEdges == 0)
 878         {
 879             break;
 880         }
 881
 882         label nPoints = edgeToPoint();
 883
 884         if (debug)
 885         {
 886             Pout<< "Total changed points      : " << nPoints << endl;
 887
 888             Pout<< "Total evaluations         : " << nEvals_ << endl;
 889
 890             Pout<< "Remaining unvisited points: " << nUnvisitedPoints_ << endl;
 891
 892             Pout<< "Remaining unvisited edges : " << nUnvisitedEdges_ << endl;
 893
 894             Pout<< endl;
 895         }
 896
 897         if (nPoints == 0)
 898         {
 899             break;
 900         }
 901
 902         iter++;
 903     }
 904
 905     return iter;
 906 }
 907
 908
 909 // ************************************************************************* //