src/meshTools/PointEdgeWave/PointEdgeWave.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "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" << abort(FatalError);
 172
 173         forAll(pointInfo, i)
 174         {
 175             pointInfo[i].transform(rotTensor[i]);
 176         }
 177     }
 178 }
 179
 180
 181 // Update info for pointI, at position pt, with information from
 182 // neighbouring edge.
 183 // Updates:
 184 //      - changedPoint_, changedPoints_, nChangedPoints_,
 185 //      - statistics: nEvals_, nUnvisitedPoints_
 186 template <class Type>
 187 bool Foam::PointEdgeWave<Type>::updatePoint
 188 (
 189     const label pointI,
 190     const label neighbourEdgeI,
 191     const Type& neighbourInfo,
 192     const scalar tol,
 193     Type& pointInfo
 194 )
 195 {
 196     nEvals_++;
 197
 198     bool wasValid = pointInfo.valid();
 199
 200     bool propagate =
 201         pointInfo.updatePoint
 202         (
 203             mesh_,
 204             pointI,
 205             neighbourEdgeI,
 206             neighbourInfo,
 207             tol
 208         );
 209
 210     if (propagate)
 211     {
 212         if (!changedPoint_[pointI])
 213         {
 214             changedPoint_[pointI] = true;
 215             changedPoints_[nChangedPoints_++] = pointI;
 216         }
 217     }
 218
 219     if (!wasValid && pointInfo.valid())
 220     {
 221         --nUnvisitedPoints_;
 222     }
 223
 224     return propagate;
 225 }
 226
 227
 228 // Update info for pointI, at position pt, with information from
 229 // same point.
 230 // Updates:
 231 //      - changedPoint_, changedPoints_, nChangedPoints_,
 232 //      - statistics: nEvals_, nUnvisitedPoints_
 233 template <class Type>
 234 bool Foam::PointEdgeWave<Type>::updatePoint
 235 (
 236     const label pointI,
 237     const Type& neighbourInfo,
 238     const scalar tol,
 239     Type& pointInfo
 240 )
 241 {
 242     nEvals_++;
 243
 244     bool wasValid = pointInfo.valid();
 245
 246     bool propagate =
 247         pointInfo.updatePoint
 248         (
 249             mesh_,
 250             pointI,
 251             neighbourInfo,
 252             tol
 253         );
 254
 255     if (propagate)
 256     {
 257         if (!changedPoint_[pointI])
 258         {
 259             changedPoint_[pointI] = true;
 260             changedPoints_[nChangedPoints_++] = pointI;
 261         }
 262     }
 263
 264     if (!wasValid && pointInfo.valid())
 265     {
 266         --nUnvisitedPoints_;
 267     }
 268
 269     return propagate;
 270 }
 271
 272
 273 // Update info for edgeI, at position pt, with information from
 274 // neighbouring point.
 275 // Updates:
 276 //      - changedEdge_, changedEdges_, nChangedEdges_,
 277 //      - statistics: nEvals_, nUnvisitedEdge_
 278 template <class Type>
 279 bool Foam::PointEdgeWave<Type>::updateEdge
 280 (
 281     const label edgeI,
 282     const label neighbourPointI,
 283     const Type& neighbourInfo,
 284     const scalar tol,
 285     Type& edgeInfo
 286 )
 287 {
 288     nEvals_++;
 289
 290     bool wasValid = edgeInfo.valid();
 291
 292     bool propagate =
 293         edgeInfo.updateEdge
 294         (
 295             mesh_,
 296             edgeI,
 297             neighbourPointI,
 298             neighbourInfo,
 299             tol
 300         );
 301
 302     if (propagate)
 303     {
 304         if (!changedEdge_[edgeI])
 305         {
 306             changedEdge_[edgeI] = true;
 307             changedEdges_[nChangedEdges_++] = edgeI;
 308         }
 309     }
 310
 311     if (!wasValid && edgeInfo.valid())
 312     {
 313         --nUnvisitedEdges_;
 314     }
 315
 316     return propagate;
 317 }
 318
 319
 320 // Check if patches of given type name are present
 321 template <class Type>
 322 template <class PatchType>
 323 Foam::label Foam::PointEdgeWave<Type>::countPatchType() const
 324 {
 325     label nPatches = 0;
 326
 327     forAll(mesh_.boundaryMesh(), patchI)
 328     {
 329         if (isA<PatchType>(mesh_.boundaryMesh()[patchI]))
 330         {
 331             nPatches++;
 332         }
 333     }
 334     return nPatches;
 335 }
 336
 337
 338 // Collect changed patch points
 339 template <class Type>
 340 void Foam::PointEdgeWave<Type>::getChangedPatchPoints
 341 (
 342     const primitivePatch& patch,
 343
 344     DynamicList<Type>& patchInfo,
 345     DynamicList<label>& patchPoints,
 346     DynamicList<label>& owner,
 347     DynamicList<label>& ownerIndex
 348 ) const
 349 {
 350     const labelList& meshPoints = patch.meshPoints();
 351     const faceList& localFaces = patch.localFaces();
 352     const labelListList& pointFaces = patch.pointFaces();
 353
 354     forAll(meshPoints, patchPointI)
 355     {
 356         label meshPointI = meshPoints[patchPointI];
 357
 358         if (changedPoint_[meshPointI])
 359         {
 360             patchInfo.append(allPointInfo_[meshPointI]);
 361             patchPoints.append(patchPointI);
 362
 363             label patchFaceI = pointFaces[patchPointI][0];
 364
 365             const face& f = localFaces[patchFaceI];
 366
 367             label index = findIndex(f, patchPointI);
 368
 369             owner.append(patchFaceI);
 370             ownerIndex.append(index);
 371         }
 372     }
 373
 374     patchInfo.shrink();
 375     patchPoints.shrink();
 376     owner.shrink();
 377     ownerIndex.shrink();
 378 }
 379
 380
 381 // Update overall for changed patch points
 382 template <class Type>
 383 void Foam::PointEdgeWave<Type>::updateFromPatchInfo
 384 (
 385     const polyPatch& meshPatch,
 386     const primitivePatch& patch,
 387     const labelList& owner,
 388     const labelList& ownerIndex,
 389     List<Type>& patchInfo
 390 )
 391 {
 392     const faceList& localFaces = patch.localFaces();
 393     const labelList& meshPoints = patch.meshPoints();
 394
 395     // Get patch and mesh points.
 396     labelList changedPatchPoints(patchInfo.size());
 397     labelList changedMeshPoints(patchInfo.size());
 398
 399     forAll(owner, i)
 400     {
 401         label faceI = owner[i];
 402
 403         const face& f = localFaces[faceI];
 404
 405         label index = (f.size() - ownerIndex[i]) % f.size();
 406
 407         changedPatchPoints[i] = f[index];
 408         changedMeshPoints[i] = meshPoints[f[index]];
 409     }
 410
 411     // Adapt for entering domain
 412     enterDomain(meshPatch, patch, changedPatchPoints, patchInfo);
 413
 414     // Merge received info
 415     forAll(patchInfo, i)
 416     {
 417         updatePoint
 418         (
 419             changedMeshPoints[i],
 420             patchInfo[i],
 421             propagationTol_,
 422             allPointInfo_[changedMeshPoints[i]]
 423         );
 424     }
 425 }
 426
 427
 428
 429 // Transfer all the information to/from neighbouring processors
 430 template <class Type>
 431 void Foam::PointEdgeWave<Type>::handleProcPatches()
 432 {
 433     // 1. Send all point info on processor patches. Send as
 434     // face label + offset in face.
 435
 436     forAll(mesh_.boundaryMesh(), patchI)
 437     {
 438         const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 439
 440         if (isA<processorPolyPatch>(patch))
 441         {
 442             // Get all changed points in relative addressing
 443
 444             DynamicList<Type> patchInfo(patch.nPoints());
 445             DynamicList<label> patchPoints(patch.nPoints());
 446             DynamicList<label> owner(patch.nPoints());
 447             DynamicList<label> ownerIndex(patch.nPoints());
 448
 449             getChangedPatchPoints
 450             (
 451                 patch,
 452                 patchInfo,
 453                 patchPoints,
 454                 owner,
 455                 ownerIndex
 456             );
 457
 458             // Adapt for leaving domain
 459             leaveDomain(patch, patch, patchPoints, patchInfo);
 460
 461             const processorPolyPatch& procPatch =
 462                 refCast<const processorPolyPatch>(patch);
 463
 464             if (debug)
 465             {
 466                 Pout<< "Processor patch " << patchI << ' ' << patch.name()
 467                     << " communicating with " << procPatch.neighbProcNo()
 468                     << "  Sending:" << patchInfo.size() << endl;
 469             }
 470
 471             {
 472                 OPstream toNeighbour
 473                 (
 474                     Pstream::blocking,
 475                     procPatch.neighbProcNo()
 476                 );
 477
 478                 toNeighbour << owner << ownerIndex << patchInfo;
 479             }
 480         }
 481     }
 482
 483
 484     //
 485     // 2. Receive all point info on processor patches.
 486     //
 487
 488     forAll(mesh_.boundaryMesh(), patchI)
 489     {
 490         const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 491
 492         if (isA<processorPolyPatch>(patch))
 493         {
 494             const processorPolyPatch& procPatch =
 495                 refCast<const processorPolyPatch>(patch);
 496
 497             List<Type> patchInfo;
 498             labelList owner;
 499             labelList ownerIndex;
 500             {
 501                 IPstream fromNeighbour
 502                 (
 503                     Pstream::blocking,
 504                     procPatch.neighbProcNo()
 505                 );
 506
 507                 fromNeighbour >> owner >> ownerIndex >> patchInfo;
 508             }
 509
 510             if (debug)
 511             {
 512                 Pout<< "Processor patch " << patchI << ' ' << patch.name()
 513                     << " communicating with " << procPatch.neighbProcNo()
 514                     << "  Received:" << patchInfo.size() << endl;
 515             }
 516
 517             // Apply transform to received data for non-parallel planes
 518             if (!procPatch.parallel())
 519             {
 520                 transform(procPatch.reverseT(), patchInfo);
 521             }
 522
 523             updateFromPatchInfo
 524             (
 525                 patch,
 526                 patch,
 527                 owner,
 528                 ownerIndex,
 529                 patchInfo
 530             );
 531         }
 532     }
 533
 534
 535
 536     //
 537     // 3. Handle all shared points
 538     //    (Note:irrespective if changed or not for now)
 539     //
 540
 541     const globalMeshData& pd = mesh_.globalData();
 542
 543     List<Type> sharedData(pd.nGlobalPoints());
 544
 545     forAll(pd.sharedPointLabels(), i)
 546     {
 547         label meshPointI = pd.sharedPointLabels()[i];
 548
 549         // Fill my entries in the shared points
 550         sharedData[pd.sharedPointAddr()[i]] = allPointInfo_[meshPointI];
 551     }
 552
 553     // Combine on master. Reduce operator has to handle a list and call
 554     // Type.updatePoint for all elements
 555     combineReduce(sharedData, listUpdateOp<Type>());
 556
 557     forAll(pd.sharedPointLabels(), i)
 558     {
 559         label meshPointI = pd.sharedPointLabels()[i];
 560
 561         // Retrieve my entries from the shared points
 562         updatePoint
 563         (
 564             meshPointI,
 565             sharedData[pd.sharedPointAddr()[i]],
 566             propagationTol_,
 567             allPointInfo_[meshPointI]
 568         );
 569     }
 570 }
 571
 572
 573 template <class Type>
 574 void Foam::PointEdgeWave<Type>::handleCyclicPatches()
 575 {
 576     // 1. Send all point info on cyclic patches. Send as
 577     // face label + offset in face.
 578
 579     label cycHalf = 0;
 580
 581     forAll(mesh_.boundaryMesh(), patchI)
 582     {
 583         const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 584
 585         if (isA<cyclicPolyPatch>(patch))
 586         {
 587             const primitivePatch& halfA = cycHalves_[cycHalf++];
 588             const primitivePatch& halfB = cycHalves_[cycHalf++];
 589
 590             // HalfA : get all changed points in relative addressing
 591
 592             DynamicList<Type> halfAInfo(halfA.nPoints());
 593             DynamicList<label> halfAPoints(halfA.nPoints());
 594             DynamicList<label> halfAOwner(halfA.nPoints());
 595             DynamicList<label> halfAIndex(halfA.nPoints());
 596
 597             getChangedPatchPoints
 598             (
 599                 halfA,
 600                 halfAInfo,
 601                 halfAPoints,
 602                 halfAOwner,
 603                 halfAIndex
 604             );
 605
 606             // HalfB : get all changed points in relative addressing
 607
 608             DynamicList<Type> halfBInfo(halfB.nPoints());
 609             DynamicList<label> halfBPoints(halfB.nPoints());
 610             DynamicList<label> halfBOwner(halfB.nPoints());
 611             DynamicList<label> halfBIndex(halfB.nPoints());
 612
 613             getChangedPatchPoints
 614             (
 615                 halfB,
 616                 halfBInfo,
 617                 halfBPoints,
 618                 halfBOwner,
 619                 halfBIndex
 620             );
 621
 622
 623             // HalfA : adapt for leaving domain
 624             leaveDomain(patch, halfA, halfAPoints, halfAInfo);
 625
 626             // HalfB : adapt for leaving domain
 627             leaveDomain(patch, halfB, halfBPoints, halfBInfo);
 628
 629
 630             // Apply rotation for non-parallel planes
 631             const cyclicPolyPatch& cycPatch =
 632                 refCast<const cyclicPolyPatch>(patch);
 633
 634             if (!cycPatch.parallel())
 635             {
 636                 // received data from half1
 637                 transform(cycPatch.forwardT(), halfAInfo);
 638
 639                 // received data from half2
 640                 transform(cycPatch.reverseT(), halfBInfo);
 641             }
 642
 643             if (debug)
 644             {
 645                 Pout<< "Cyclic patch " << patchI << ' ' << patch.name()
 646                     << "  Changed on first half : " << halfAInfo.size()
 647                     << "  Changed on second half : " << halfBInfo.size()
 648                     << endl;
 649             }
 650
 651             // Half1: update with data from halfB
 652             updateFromPatchInfo
 653             (
 654                 patch,
 655                 halfA,
 656                 halfBOwner,
 657                 halfBIndex,
 658                 halfBInfo
 659             );
 660
 661             // Half2: update with data from halfA
 662             updateFromPatchInfo
 663             (
 664                 patch,
 665                 halfB,
 666                 halfAOwner,
 667                 halfAIndex,
 668                 halfAInfo
 669             );
 670
 671             if (debug)
 672             {
 673                 //checkCyclic(patch);
 674             }
 675         }
 676     }
 677 }
 678
 679
 680 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 681
 682 // Iterate, propagating changedPointsInfo across mesh, until no change (or
 683 // maxIter reached). Initial point values specified.
 684 template <class Type>
 685 Foam::PointEdgeWave<Type>::PointEdgeWave
 686 (
 687     const polyMesh& mesh,
 688     const labelList& changedPoints,
 689     const List<Type>& changedPointsInfo,
 690
 691     List<Type>& allPointInfo,
 692     List<Type>& allEdgeInfo,
 693
 694     const label maxIter
 695 )
 696 :
 697     mesh_(mesh),
 698     allPointInfo_(allPointInfo),
 699     allEdgeInfo_(allEdgeInfo),
 700     changedPoint_(mesh_.nPoints(), false),
 701     changedPoints_(mesh_.nPoints()),
 702     nChangedPoints_(0),
 703     changedEdge_(mesh_.nEdges(), false),
 704     changedEdges_(mesh_.nEdges()),
 705     nChangedEdges_(0),
 706     nCyclicPatches_(countPatchType<cyclicPolyPatch>()),
 707     cycHalves_(2*nCyclicPatches_),
 708     nEvals_(0),
 709     nUnvisitedPoints_(mesh_.nPoints()),
 710     nUnvisitedEdges_(mesh_.nEdges())
 711 {
 712     if (allPointInfo_.size() != mesh_.nPoints())
 713     {
 714         FatalErrorIn
 715         (
 716             "PointEdgeWave<Type>::PointEdgeWave"
 717             "(const polyMesh&, const labelList&, const List<Type>,"
 718             " List<Type>&, List<Type>&, const label maxIter)"
 719         )   << "size of pointInfo work array is not equal to the number"
 720             << " of points in the mesh" << endl
 721             << "    pointInfo   :" << allPointInfo_.size() << endl
 722             << "    mesh.nPoints:" << mesh_.nPoints()
 723             << exit(FatalError);
 724     }
 725     if (allEdgeInfo_.size() != mesh_.nEdges())
 726     {
 727         FatalErrorIn
 728         (
 729             "PointEdgeWave<Type>::PointEdgeWave"
 730             "(const polyMesh&, const labelList&, const List<Type>,"
 731             " List<Type>&, List<Type>&, const label maxIter)"
 732         )   << "size of edgeInfo work array is not equal to the number"
 733             << " of edges in the mesh" << endl
 734             << "    edgeInfo   :" << allEdgeInfo_.size() << endl
 735             << "    mesh.nEdges:" << mesh_.nEdges()
 736             << exit(FatalError);
 737     }
 738
 739
 740     // Calculate cyclic halves addressing.
 741     if (nCyclicPatches_ > 0)
 742     {
 743         calcCyclicAddressing();
 744     }
 745
 746
 747     // Set from initial changed points data
 748     setPointInfo(changedPoints, changedPointsInfo);
 749
 750     if (debug)
 751     {
 752         Pout<< "Seed points               : " << nChangedPoints_ << endl;
 753     }
 754
 755     // Iterate until nothing changes
 756     label iter = iterate(maxIter);
 757
 758     if ((maxIter > 0) && (iter >= maxIter))
 759     {
 760         FatalErrorIn
 761         (
 762             "PointEdgeWave<Type>::PointEdgeWave"
 763             "(const polyMesh&, const labelList&, const List<Type>,"
 764             " List<Type>&, List<Type>&, const label maxIter)"
 765         )   << "Maximum number of iterations reached. Increase maxIter." << endl
 766             << "    maxIter:" << maxIter << endl
 767             << "    nChangedPoints:" << nChangedPoints_ << endl
 768             << "    nChangedEdges:" << nChangedEdges_ << endl
 769             << exit(FatalError);
 770     }
 771 }
 772
 773
 774 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 775
 776 template <class Type>
 777 Foam::PointEdgeWave<Type>::~PointEdgeWave()
 778 {}
 779
 780
 781 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 782
 783
 784 template <class Type>
 785 Foam::label Foam::PointEdgeWave<Type>::getUnsetPoints() const
 786 {
 787     return nUnvisitedPoints_;
 788 }
 789
 790
 791 template <class Type>
 792 Foam::label Foam::PointEdgeWave<Type>::getUnsetEdges() const
 793 {
 794     return nUnvisitedEdges_;
 795 }
 796
 797
 798 // Copy point information into member data
 799 template <class Type>
 800 void Foam::PointEdgeWave<Type>::setPointInfo
 801 (
 802     const labelList& changedPoints,
 803     const List<Type>& changedPointsInfo
 804 )
 805 {
 806     forAll(changedPoints, changedPointI)
 807     {
 808         label pointI = changedPoints[changedPointI];
 809
 810         bool wasValid = allPointInfo_[pointI].valid();
 811
 812         // Copy info for pointI
 813         allPointInfo_[pointI] = changedPointsInfo[changedPointI];
 814
 815         // Maintain count of unset points
 816         if (!wasValid && allPointInfo_[pointI].valid())
 817         {
 818             --nUnvisitedPoints_;
 819         }
 820
 821         // Mark pointI as changed, both on list and on point itself.
 822
 823         if (!changedPoint_[pointI])
 824         {
 825             changedPoint_[pointI] = true;
 826             changedPoints_[nChangedPoints_++] = pointI;
 827         }
 828     }
 829 }
 830
 831
 832 // Propagate information from edge to point. Return number of points changed.
 833 template <class Type>
 834 Foam::label Foam::PointEdgeWave<Type>::edgeToPoint()
 835 {
 836     for
 837     (
 838         label changedEdgeI = 0;
 839         changedEdgeI < nChangedEdges_;
 840         changedEdgeI++
 841     )
 842     {
 843         label edgeI = changedEdges_[changedEdgeI];
 844
 845         if (!changedEdge_[edgeI])
 846         {
 847             FatalErrorIn("PointEdgeWave<Type>::edgeToPoint()")
 848                 << "edge " << edgeI
 849                 << " not marked as having been changed" << nl
 850                 << "This might be caused by multiple occurences of the same"
 851                 << " seed point." << abort(FatalError);
 852         }
 853
 854
 855         const Type& neighbourWallInfo = allEdgeInfo_[edgeI];
 856
 857         // Evaluate all connected points (= edge endpoints)
 858         const edge& e = mesh_.edges()[edgeI];
 859
 860         forAll(e, eI)
 861         {
 862             Type& currentWallInfo = allPointInfo_[e[eI]];
 863
 864             if (currentWallInfo != neighbourWallInfo)
 865             {
 866                 updatePoint
 867                 (
 868                     e[eI],
 869                     edgeI,
 870                     neighbourWallInfo,
 871                     propagationTol_,
 872                     currentWallInfo
 873                 );
 874             }
 875         }
 876
 877         // Reset status of edge
 878         changedEdge_[edgeI] = false;
 879     }
 880
 881     // Handled all changed edges by now
 882     nChangedEdges_ = 0;
 883
 884     if (nCyclicPatches_ > 0)
 885     {
 886         // Transfer changed points across cyclic halves
 887         handleCyclicPatches();
 888     }
 889     if (Pstream::parRun())
 890     {
 891         // Transfer changed points from neighbouring processors.
 892         handleProcPatches();
 893     }
 894
 895     if (debug)
 896     {
 897         Pout<< "Changed points            : " << nChangedPoints_ << endl;
 898     }
 899
 900     // Sum nChangedPoints over all procs
 901     label totNChanged = nChangedPoints_;
 902
 903     reduce(totNChanged, sumOp<label>());
 904
 905     return totNChanged;
 906 }
 907
 908
 909 // Propagate information from point to edge. Return number of edges changed.
 910 template <class Type>
 911 Foam::label Foam::PointEdgeWave<Type>::pointToEdge()
 912 {
 913     const labelListList& pointEdges = mesh_.pointEdges();
 914
 915     for
 916     (
 917         label changedPointI = 0;
 918         changedPointI < nChangedPoints_;
 919         changedPointI++
 920     )
 921     {
 922         label pointI = changedPoints_[changedPointI];
 923
 924         if (!changedPoint_[pointI])
 925         {
 926             FatalErrorIn("PointEdgeWave<Type>::pointToEdge()")
 927                 << "Point " << pointI
 928                 << " not marked as having been changed" << nl
 929                 << "This might be caused by multiple occurences of the same"
 930                 << " seed point." << abort(FatalError);
 931         }
 932
 933         const Type& neighbourWallInfo = allPointInfo_[pointI];
 934
 935         // Evaluate all connected edges
 936
 937         const labelList& edgeLabels = pointEdges[pointI];
 938         forAll(edgeLabels, edgeLabelI)
 939         {
 940             label edgeI = edgeLabels[edgeLabelI];
 941
 942             Type& currentWallInfo = allEdgeInfo_[edgeI];
 943
 944             if (currentWallInfo != neighbourWallInfo)
 945             {
 946                 updateEdge
 947                 (
 948                     edgeI,
 949                     pointI,
 950                     neighbourWallInfo,
 951                     propagationTol_,
 952                     currentWallInfo
 953                 );
 954             }
 955         }
 956
 957         // Reset status of point
 958         changedPoint_[pointI] = false;
 959     }
 960
 961     // Handled all changed points by now
 962     nChangedPoints_ = 0;
 963
 964     if (debug)
 965     {
 966         Pout<< "Changed edges             : " << nChangedEdges_ << endl;
 967     }
 968
 969     // Sum nChangedPoints over all procs
 970     label totNChanged = nChangedEdges_;
 971
 972     reduce(totNChanged, sumOp<label>());
 973
 974     return totNChanged;
 975 }
 976
 977
 978 // Iterate
 979 template <class Type>
 980 Foam::label Foam::PointEdgeWave<Type>::iterate(const label maxIter)
 981 {
 982     if (nCyclicPatches_ > 0)
 983     {
 984         // Transfer changed points across cyclic halves
 985         handleCyclicPatches();
 986     }
 987     if (Pstream::parRun())
 988     {
 989         // Transfer changed points from neighbouring processors.
 990         handleProcPatches();
 991     }
 992
 993     nEvals_ = 0;
 994
 995     label iter = 0;
 996
 997     while(iter < maxIter)
 998     {
 999         if (debug)
1000         {
1001             Pout<< "Iteration " << iter << endl;
1002         }
1003
1004         label nEdges = pointToEdge();
1005
1006         if (debug)
1007         {
1008             Pout<< "Total changed edges       : " << nEdges << endl;
1009         }
1010
1011         if (nEdges == 0)
1012         {
1013             break;
1014         }
1015
1016         label nPoints = edgeToPoint();
1017
1018         if (debug)
1019         {
1020             Pout<< "Total changed points      : " << nPoints << endl;
1021
1022             Pout<< "Total evaluations         : " << nEvals_ << endl;
1023
1024             Pout<< "Remaining unvisited points: " << nUnvisitedPoints_ << endl;
1025
1026             Pout<< "Remaining unvisited edges : " << nUnvisitedEdges_ << endl;
1027
1028             Pout<< endl;
1029         }
1030
1031         if (nPoints == 0)
1032         {
1033             break;
1034         }
1035
1036         iter++;
1037     }
1038
1039     return iter;
1040 }
1041
1042 // ************************************************************************* //