src/meshTools/PointEdgeWave/PointEdgeWave.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 "PointEdgeWave.H"
  28 #include "polyMesh.H"
  29 #include "processorPolyPatch.H"
  30 #include "cyclicPolyPatch.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     DynamicList<label>& patchPoints,
 347     DynamicList<label>& owner,
 348     DynamicList<label>& 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             patchPoints.append(patchPointI);
 363
 364             label patchFaceI = pointFaces[patchPointI][0];
 365
 366             const face& f = localFaces[patchFaceI];
 367
 368             label index = findIndex(f, patchPointI);
 369
 370             owner.append(patchFaceI);
 371             ownerIndex.append(index);
 372         }
 373     }
 374
 375     patchInfo.shrink();
 376     patchPoints.shrink();
 377     owner.shrink();
 378     ownerIndex.shrink();
 379 }
 380
 381
 382 // Update overall for changed patch points
 383 template <class Type>
 384 void Foam::PointEdgeWave<Type>::updateFromPatchInfo
 385 (
 386     const polyPatch& meshPatch,
 387     const primitivePatch& patch,
 388     const labelList& owner,
 389     const labelList& ownerIndex,
 390     List<Type>& patchInfo
 391 )
 392 {
 393     const faceList& localFaces = patch.localFaces();
 394     const labelList& meshPoints = patch.meshPoints();
 395
 396     // Get patch and mesh points.
 397     labelList changedPatchPoints(patchInfo.size());
 398     labelList changedMeshPoints(patchInfo.size());
 399
 400     forAll(owner, i)
 401     {
 402         label faceI = owner[i];
 403
 404         const face& f = localFaces[faceI];
 405
 406         label index = (f.size() - ownerIndex[i]) % f.size();
 407
 408         changedPatchPoints[i] = f[index];
 409         changedMeshPoints[i] = meshPoints[f[index]];
 410     }
 411
 412     // Adapt for entering domain
 413     enterDomain(meshPatch, patch, changedPatchPoints, patchInfo);
 414
 415     // Merge received info
 416     forAll(patchInfo, i)
 417     {
 418         updatePoint
 419         (
 420             changedMeshPoints[i],
 421             patchInfo[i],
 422             propagationTol_,
 423             allPointInfo_[changedMeshPoints[i]]
 424         );
 425     }
 426 }
 427
 428
 429
 430 // Transfer all the information to/from neighbouring processors
 431 template <class Type>
 432 void Foam::PointEdgeWave<Type>::handleProcPatches()
 433 {
 434     // 1. Send all point info on processor patches. Send as
 435     // face label + offset in face.
 436
 437     forAll(mesh_.boundaryMesh(), patchI)
 438     {
 439         const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 440
 441         if (isA<processorPolyPatch>(patch))
 442         {
 443             // Get all changed points in relative addressing
 444
 445             DynamicList<Type> patchInfo(patch.nPoints());
 446             DynamicList<label> patchPoints(patch.nPoints());
 447             DynamicList<label> owner(patch.nPoints());
 448             DynamicList<label> ownerIndex(patch.nPoints());
 449
 450             getChangedPatchPoints
 451             (
 452                 patch,
 453                 patchInfo,
 454                 patchPoints,
 455                 owner,
 456                 ownerIndex
 457             );
 458
 459             // Adapt for leaving domain
 460             leaveDomain(patch, patch, patchPoints, patchInfo);
 461
 462             const processorPolyPatch& procPatch =
 463                 refCast<const processorPolyPatch>(patch);
 464
 465             if (debug)
 466             {
 467                 Pout<< "Processor patch " << patchI << ' ' << patch.name()
 468                     << " communicating with " << procPatch.neighbProcNo()
 469                     << "  Sending:" << patchInfo.size() << endl;
 470             }
 471
 472             {
 473                 OPstream toNeighbour
 474                 (
 475                     Pstream::blocking,
 476                     procPatch.neighbProcNo()
 477                 );
 478
 479                 toNeighbour << owner << ownerIndex << patchInfo;
 480             }
 481         }
 482     }
 483
 484
 485     //
 486     // 2. Receive all point info on processor patches.
 487     //
 488
 489     forAll(mesh_.boundaryMesh(), patchI)
 490     {
 491         const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 492
 493         if (isA<processorPolyPatch>(patch))
 494         {
 495             const processorPolyPatch& procPatch =
 496                 refCast<const processorPolyPatch>(patch);
 497
 498             List<Type> patchInfo;
 499             labelList owner;
 500             labelList ownerIndex;
 501             {
 502                 IPstream fromNeighbour
 503                 (
 504                     Pstream::blocking,
 505                     procPatch.neighbProcNo()
 506                 );
 507
 508                 fromNeighbour >> owner >> ownerIndex >> patchInfo;
 509             }
 510
 511             if (debug)
 512             {
 513                 Pout<< "Processor patch " << patchI << ' ' << patch.name()
 514                     << " communicating with " << procPatch.neighbProcNo()
 515                     << "  Received:" << patchInfo.size() << endl;
 516             }
 517
 518             // Apply transform to received data for non-parallel planes
 519             if (!procPatch.parallel())
 520             {
 521                 transform(procPatch.forwardT(), patchInfo);
 522             }
 523
 524             updateFromPatchInfo
 525             (
 526                 patch,
 527                 patch,
 528                 owner,
 529                 ownerIndex,
 530                 patchInfo
 531             );
 532         }
 533     }
 534
 535
 536
 537     //
 538     // 3. Handle all shared points
 539     //    (Note:irrespective if changed or not for now)
 540     //
 541
 542     const globalMeshData& pd = mesh_.globalData();
 543
 544     List<Type> sharedData(pd.nGlobalPoints());
 545
 546     forAll(pd.sharedPointLabels(), i)
 547     {
 548         label meshPointI = pd.sharedPointLabels()[i];
 549
 550         // Fill my entries in the shared points
 551         sharedData[pd.sharedPointAddr()[i]] = allPointInfo_[meshPointI];
 552     }
 553
 554     // Combine on master. Reduce operator has to handle a list and call
 555     // Type.updatePoint for all elements
 556     combineReduce(sharedData, listUpdateOp<Type>());
 557
 558     forAll(pd.sharedPointLabels(), i)
 559     {
 560         label meshPointI = pd.sharedPointLabels()[i];
 561
 562         // Retrieve my entries from the shared points
 563         updatePoint
 564         (
 565             meshPointI,
 566             sharedData[pd.sharedPointAddr()[i]],
 567             propagationTol_,
 568             allPointInfo_[meshPointI]
 569         );
 570     }
 571 }
 572
 573
 574 template <class Type>
 575 void Foam::PointEdgeWave<Type>::handleCyclicPatches()
 576 {
 577     // 1. Send all point info on cyclic patches. Send as
 578     // face label + offset in face.
 579
 580     label cycHalf = 0;
 581
 582     forAll(mesh_.boundaryMesh(), patchI)
 583     {
 584         const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 585
 586         if (isA<cyclicPolyPatch>(patch))
 587         {
 588             const primitivePatch& halfA = cycHalves_[cycHalf++];
 589             const primitivePatch& halfB = cycHalves_[cycHalf++];
 590
 591             // HalfA : get all changed points in relative addressing
 592
 593             DynamicList<Type> halfAInfo(halfA.nPoints());
 594             DynamicList<label> halfAPoints(halfA.nPoints());
 595             DynamicList<label> halfAOwner(halfA.nPoints());
 596             DynamicList<label> halfAIndex(halfA.nPoints());
 597
 598             getChangedPatchPoints
 599             (
 600                 halfA,
 601                 halfAInfo,
 602                 halfAPoints,
 603                 halfAOwner,
 604                 halfAIndex
 605             );
 606
 607             // HalfB : get all changed points in relative addressing
 608
 609             DynamicList<Type> halfBInfo(halfB.nPoints());
 610             DynamicList<label> halfBPoints(halfB.nPoints());
 611             DynamicList<label> halfBOwner(halfB.nPoints());
 612             DynamicList<label> halfBIndex(halfB.nPoints());
 613
 614             getChangedPatchPoints
 615             (
 616                 halfB,
 617                 halfBInfo,
 618                 halfBPoints,
 619                 halfBOwner,
 620                 halfBIndex
 621             );
 622
 623
 624             // HalfA : adapt for leaving domain
 625             leaveDomain(patch, halfA, halfAPoints, halfAInfo);
 626
 627             // HalfB : adapt for leaving domain
 628             leaveDomain(patch, halfB, halfBPoints, halfBInfo);
 629
 630
 631             // Apply rotation for non-parallel planes
 632             const cyclicPolyPatch& cycPatch =
 633                 refCast<const cyclicPolyPatch>(patch);
 634
 635             if (!cycPatch.parallel())
 636             {
 637                 // received data from half1
 638                 transform(cycPatch.forwardT(), halfAInfo);
 639
 640                 // received data from half2
 641                 transform(cycPatch.forwardT(), halfBInfo);
 642             }
 643
 644             if (debug)
 645             {
 646                 Pout<< "Cyclic patch " << patchI << ' ' << patch.name()
 647                     << "  Changed on first half : " << halfAInfo.size()
 648                     << "  Changed on second half : " << halfBInfo.size()
 649                     << endl;
 650             }
 651
 652             // Half1: update with data from halfB
 653             updateFromPatchInfo
 654             (
 655                 patch,
 656                 halfA,
 657                 halfBOwner,
 658                 halfBIndex,
 659                 halfBInfo
 660             );
 661
 662             // Half2: update with data from halfA
 663             updateFromPatchInfo
 664             (
 665                 patch,
 666                 halfB,
 667                 halfAOwner,
 668                 halfAIndex,
 669                 halfAInfo
 670             );
 671
 672             if (debug)
 673             {
 674                 //checkCyclic(patch);
 675             }
 676         }
 677     }
 678 }
 679
 680
 681 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 682
 683 // Iterate, propagating changedPointsInfo across mesh, until no change (or
 684 // maxIter reached). Initial point values specified.
 685 template <class Type>
 686 Foam::PointEdgeWave<Type>::PointEdgeWave
 687 (
 688     const polyMesh& mesh,
 689     const labelList& changedPoints,
 690     const List<Type>& changedPointsInfo,
 691
 692     List<Type>& allPointInfo,
 693     List<Type>& allEdgeInfo,
 694
 695     const label maxIter
 696 )
 697 :
 698     mesh_(mesh),
 699     allPointInfo_(allPointInfo),
 700     allEdgeInfo_(allEdgeInfo),
 701     changedPoint_(mesh_.nPoints(), false),
 702     changedPoints_(mesh_.nPoints()),
 703     nChangedPoints_(0),
 704     changedEdge_(mesh_.nEdges(), false),
 705     changedEdges_(mesh_.nEdges()),
 706     nChangedEdges_(0),
 707     nCyclicPatches_(countPatchType<cyclicPolyPatch>()),
 708     cycHalves_(2*nCyclicPatches_),
 709     nEvals_(0),
 710     nUnvisitedPoints_(mesh_.nPoints()),
 711     nUnvisitedEdges_(mesh_.nEdges())
 712 {
 713     if (allPointInfo_.size() != mesh_.nPoints())
 714     {
 715         FatalErrorIn
 716         (
 717             "PointEdgeWave<Type>::PointEdgeWave"
 718             "(const polyMesh&, const labelList&, const List<Type>,"
 719             " List<Type>&, List<Type>&, const label maxIter)"
 720         )   << "size of pointInfo work array is not equal to the number"
 721             << " of points in the mesh" << endl
 722             << "    pointInfo   :" << allPointInfo_.size() << endl
 723             << "    mesh.nPoints:" << mesh_.nPoints()
 724             << exit(FatalError);
 725     }
 726     if (allEdgeInfo_.size() != mesh_.nEdges())
 727     {
 728         FatalErrorIn
 729         (
 730             "PointEdgeWave<Type>::PointEdgeWave"
 731             "(const polyMesh&, const labelList&, const List<Type>,"
 732             " List<Type>&, List<Type>&, const label maxIter)"
 733         )   << "size of edgeInfo work array is not equal to the number"
 734             << " of edges in the mesh" << endl
 735             << "    edgeInfo   :" << allEdgeInfo_.size() << endl
 736             << "    mesh.nEdges:" << mesh_.nEdges()
 737             << exit(FatalError);
 738     }
 739
 740
 741     // Calculate cyclic halves addressing.
 742     if (nCyclicPatches_ > 0)
 743     {
 744         calcCyclicAddressing();
 745     }
 746
 747
 748     // Set from initial changed points data
 749     setPointInfo(changedPoints, changedPointsInfo);
 750
 751     if (debug)
 752     {
 753         Pout<< "Seed points               : " << nChangedPoints_ << endl;
 754     }
 755
 756     // Iterate until nothing changes
 757     label iter = iterate(maxIter);
 758
 759     if ((maxIter > 0) && (iter >= maxIter))
 760     {
 761         FatalErrorIn
 762         (
 763             "PointEdgeWave<Type>::PointEdgeWave"
 764             "(const polyMesh&, const labelList&, const List<Type>,"
 765             " List<Type>&, List<Type>&, const label maxIter)"
 766         )   << "Maximum number of iterations reached. Increase maxIter." << nl
 767             << "    maxIter:" << maxIter << endl
 768             << "    nChangedPoints:" << nChangedPoints_ << endl
 769             << "    nChangedEdges:" << nChangedEdges_ << endl
 770             << exit(FatalError);
 771     }
 772 }
 773
 774
 775 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 776
 777 template <class Type>
 778 Foam::PointEdgeWave<Type>::~PointEdgeWave()
 779 {}
 780
 781
 782 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 783
 784
 785 template <class Type>
 786 Foam::label Foam::PointEdgeWave<Type>::getUnsetPoints() const
 787 {
 788     return nUnvisitedPoints_;
 789 }
 790
 791
 792 template <class Type>
 793 Foam::label Foam::PointEdgeWave<Type>::getUnsetEdges() const
 794 {
 795     return nUnvisitedEdges_;
 796 }
 797
 798
 799 // Copy point information into member data
 800 template <class Type>
 801 void Foam::PointEdgeWave<Type>::setPointInfo
 802 (
 803     const labelList& changedPoints,
 804     const List<Type>& changedPointsInfo
 805 )
 806 {
 807     forAll(changedPoints, changedPointI)
 808     {
 809         label pointI = changedPoints[changedPointI];
 810
 811         bool wasValid = allPointInfo_[pointI].valid();
 812
 813         // Copy info for pointI
 814         allPointInfo_[pointI] = changedPointsInfo[changedPointI];
 815
 816         // Maintain count of unset points
 817         if (!wasValid && allPointInfo_[pointI].valid())
 818         {
 819             --nUnvisitedPoints_;
 820         }
 821
 822         // Mark pointI as changed, both on list and on point itself.
 823
 824         if (!changedPoint_[pointI])
 825         {
 826             changedPoint_[pointI] = true;
 827             changedPoints_[nChangedPoints_++] = pointI;
 828         }
 829     }
 830 }
 831
 832
 833 // Propagate information from edge to point. Return number of points changed.
 834 template <class Type>
 835 Foam::label Foam::PointEdgeWave<Type>::edgeToPoint()
 836 {
 837     for
 838     (
 839         label changedEdgeI = 0;
 840         changedEdgeI < nChangedEdges_;
 841         changedEdgeI++
 842     )
 843     {
 844         label edgeI = changedEdges_[changedEdgeI];
 845
 846         if (!changedEdge_[edgeI])
 847         {
 848             FatalErrorIn("PointEdgeWave<Type>::edgeToPoint()")
 849                 << "edge " << edgeI
 850                 << " not marked as having been changed" << nl
 851                 << "This might be caused by multiple occurences of the same"
 852                 << " seed point." << abort(FatalError);
 853         }
 854
 855
 856         const Type& neighbourWallInfo = allEdgeInfo_[edgeI];
 857
 858         // Evaluate all connected points (= edge endpoints)
 859         const edge& e = mesh_.edges()[edgeI];
 860
 861         forAll(e, eI)
 862         {
 863             Type& currentWallInfo = allPointInfo_[e[eI]];
 864
 865             if (currentWallInfo != neighbourWallInfo)
 866             {
 867                 updatePoint
 868                 (
 869                     e[eI],
 870                     edgeI,
 871                     neighbourWallInfo,
 872                     propagationTol_,
 873                     currentWallInfo
 874                 );
 875             }
 876         }
 877
 878         // Reset status of edge
 879         changedEdge_[edgeI] = false;
 880     }
 881
 882     // Handled all changed edges by now
 883     nChangedEdges_ = 0;
 884
 885     if (nCyclicPatches_ > 0)
 886     {
 887         // Transfer changed points across cyclic halves
 888         handleCyclicPatches();
 889     }
 890     if (Pstream::parRun())
 891     {
 892         // Transfer changed points from neighbouring processors.
 893         handleProcPatches();
 894     }
 895
 896     if (debug)
 897     {
 898         Pout<< "Changed points            : " << nChangedPoints_ << endl;
 899     }
 900
 901     // Sum nChangedPoints over all procs
 902     label totNChanged = nChangedPoints_;
 903
 904     reduce(totNChanged, sumOp<label>());
 905
 906     return totNChanged;
 907 }
 908
 909
 910 // Propagate information from point to edge. Return number of edges changed.
 911 template <class Type>
 912 Foam::label Foam::PointEdgeWave<Type>::pointToEdge()
 913 {
 914     const labelListList& pointEdges = mesh_.pointEdges();
 915
 916     for
 917     (
 918         label changedPointI = 0;
 919         changedPointI < nChangedPoints_;
 920         changedPointI++
 921     )
 922     {
 923         label pointI = changedPoints_[changedPointI];
 924
 925         if (!changedPoint_[pointI])
 926         {
 927             FatalErrorIn("PointEdgeWave<Type>::pointToEdge()")
 928                 << "Point " << pointI
 929                 << " not marked as having been changed" << nl
 930                 << "This might be caused by multiple occurences of the same"
 931                 << " seed point." << abort(FatalError);
 932         }
 933
 934         const Type& neighbourWallInfo = allPointInfo_[pointI];
 935
 936         // Evaluate all connected edges
 937
 938         const labelList& edgeLabels = pointEdges[pointI];
 939         forAll(edgeLabels, edgeLabelI)
 940         {
 941             label edgeI = edgeLabels[edgeLabelI];
 942
 943             Type& currentWallInfo = allEdgeInfo_[edgeI];
 944
 945             if (currentWallInfo != neighbourWallInfo)
 946             {
 947                 updateEdge
 948                 (
 949                     edgeI,
 950                     pointI,
 951                     neighbourWallInfo,
 952                     propagationTol_,
 953                     currentWallInfo
 954                 );
 955             }
 956         }
 957
 958         // Reset status of point
 959         changedPoint_[pointI] = false;
 960     }
 961
 962     // Handled all changed points by now
 963     nChangedPoints_ = 0;
 964
 965     if (debug)
 966     {
 967         Pout<< "Changed edges             : " << nChangedEdges_ << endl;
 968     }
 969
 970     // Sum nChangedPoints over all procs
 971     label totNChanged = nChangedEdges_;
 972
 973     reduce(totNChanged, sumOp<label>());
 974
 975     return totNChanged;
 976 }
 977
 978
 979 // Iterate
 980 template <class Type>
 981 Foam::label Foam::PointEdgeWave<Type>::iterate(const label maxIter)
 982 {
 983     if (nCyclicPatches_ > 0)
 984     {
 985         // Transfer changed points across cyclic halves
 986         handleCyclicPatches();
 987     }
 988     if (Pstream::parRun())
 989     {
 990         // Transfer changed points from neighbouring processors.
 991         handleProcPatches();
 992     }
 993
 994     nEvals_ = 0;
 995
 996     label iter = 0;
 997
 998     while(iter < maxIter)
 999     {
1000         if (debug)
1001         {
1002             Pout<< "Iteration " << iter << endl;
1003         }
1004
1005         label nEdges = pointToEdge();
1006
1007         if (debug)
1008         {
1009             Pout<< "Total changed edges       : " << nEdges << endl;
1010         }
1011
1012         if (nEdges == 0)
1013         {
1014             break;
1015         }
1016
1017         label nPoints = edgeToPoint();
1018
1019         if (debug)
1020         {
1021             Pout<< "Total changed points      : " << nPoints << endl;
1022
1023             Pout<< "Total evaluations         : " << nEvals_ << endl;
1024
1025             Pout<< "Remaining unvisited points: " << nUnvisitedPoints_ << endl;
1026
1027             Pout<< "Remaining unvisited edges : " << nUnvisitedEdges_ << endl;
1028
1029             Pout<< endl;
1030         }
1031
1032         if (nPoints == 0)
1033         {
1034             break;
1035         }
1036
1037         iter++;
1038     }
1039
1040     return iter;
1041 }
1042
1043
1044 // ************************************************************************* //