src/foam/meshes/polyMesh/syncTools/syncToolsTemplates.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 "syncTools.H"
  27 #include "polyMesh.H"
  28 #include "processorPolyPatch.H"
  29 #include "cyclicPolyPatch.H"
  30 #include "globalMeshData.H"
  31 #include "contiguous.H"
  32 #include "transform.H"
  33
  34 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  35
  36 template <class T>
  37 void Foam::syncTools::separateList
  38 (
  39     const vectorField& separation,
  40     UList<T>& field
  41 )
  42 {}
  43
  44
  45 template <class T>
  46 void Foam::syncTools::separateList
  47 (
  48     const vectorField& separation,
  49     Map<T>& field
  50 )
  51 {}
  52
  53
  54 template <class T>
  55 void Foam::syncTools::separateList
  56 (
  57     const vectorField& separation,
  58     EdgeMap<T>& field
  59 )
  60 {}
  61
  62
  63 // Combine val with existing value at index
  64 template <class T, class CombineOp>
  65 void Foam::syncTools::combine
  66 (
  67     Map<T>& pointValues,
  68     const CombineOp& cop,
  69     const label index,
  70     const T& val
  71 )
  72 {
  73     typename Map<T>::iterator iter = pointValues.find(index);
  74
  75     if (iter != pointValues.end())
  76     {
  77         cop(iter(), val);
  78     }
  79     else
  80     {
  81         pointValues.insert(index, val);
  82     }
  83 }
  84
  85
  86 // Combine val with existing value at (implicit index) e.
  87 template <class T, class CombineOp>
  88 void Foam::syncTools::combine
  89 (
  90     EdgeMap<T>& edgeValues,
  91     const CombineOp& cop,
  92     const edge& index,
  93     const T& val
  94 )
  95 {
  96     typename EdgeMap<T>::iterator iter = edgeValues.find(index);
  97
  98     if (iter != edgeValues.end())
  99     {
 100         cop(iter(), val);
 101     }
 102     else
 103     {
 104         edgeValues.insert(index, val);
 105     }
 106 }
 107
 108
 109 template <class T, class CombineOp>
 110 void Foam::syncTools::syncPointMap
 111 (
 112     const polyMesh& mesh,
 113     Map<T>& pointValues,        // from mesh point label to value
 114     const CombineOp& cop,
 115     const bool applySeparation
 116 )
 117 {
 118     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 119
 120     if (!hasCouples(patches))
 121     {
 122         return;
 123     }
 124
 125     // Is there any coupled patch with transformation?
 126     bool hasTransformation = false;
 127
 128     if (Pstream::parRun())
 129     {
 130         // Send
 131
 132         forAll(patches, patchI)
 133         {
 134             if
 135             (
 136                 isA<processorPolyPatch>(patches[patchI])
 137              && patches[patchI].nPoints() > 0
 138             )
 139             {
 140                 const processorPolyPatch& procPatch =
 141                     refCast<const processorPolyPatch>(patches[patchI]);
 142
 143                 // Get data per patchPoint in neighbouring point numbers.
 144
 145                 const labelList& meshPts = procPatch.meshPoints();
 146                 const labelList& nbrPts = procPatch.neighbPoints();
 147
 148                 // Extract local values. Create map from nbrPoint to value.
 149                 // Note: how small initial size?
 150                 Map<T> patchInfo(meshPts.size() / 20);
 151
 152                 forAll(meshPts, i)
 153                 {
 154                     typename Map<T>::const_iterator iter =
 155                         pointValues.find(meshPts[i]);
 156
 157                     if (iter != pointValues.end())
 158                     {
 159                         if (nbrPts[i] >= 0)
 160                         {
 161                             patchInfo.insert(nbrPts[i], iter());
 162                         }
 163                     }
 164                 }
 165
 166                 OPstream toNeighb(Pstream::blocking, procPatch.neighbProcNo());
 167                 toNeighb << patchInfo;
 168             }
 169         }
 170
 171
 172         // Receive and combine.
 173
 174         forAll(patches, patchI)
 175         {
 176             if
 177             (
 178                 isA<processorPolyPatch>(patches[patchI])
 179              && patches[patchI].nPoints() > 0
 180             )
 181             {
 182                 const processorPolyPatch& procPatch =
 183                     refCast<const processorPolyPatch>(patches[patchI]);
 184                 checkTransform(procPatch, applySeparation);
 185
 186                 IPstream fromNb(Pstream::blocking, procPatch.neighbProcNo());
 187                 Map<T> nbrPatchInfo(fromNb);
 188
 189                 if (!procPatch.parallel())
 190                 {
 191                     hasTransformation = true;
 192                     transformList(procPatch.forwardT(), nbrPatchInfo);
 193                 }
 194                 else if (applySeparation && procPatch.separated())
 195                 {
 196                     hasTransformation = true;
 197                     separateList(-procPatch.separation(), nbrPatchInfo);
 198                 }
 199
 200                 const labelList& meshPts = procPatch.meshPoints();
 201
 202                 // Only update those values which come from neighbour
 203
 204                 forAllConstIter
 205                 (
 206                     typename Map<T>,
 207                     nbrPatchInfo,
 208                     nbrIter
 209                 )
 210                 {
 211                     combine
 212                     (
 213                         pointValues,
 214                         cop,
 215                         meshPts[nbrIter.key()],
 216                         nbrIter()
 217                     );
 218                 }
 219             }
 220         }
 221     }
 222
 223     // Do the cyclics.
 224     forAll(patches, patchI)
 225     {
 226         if (isA<cyclicPolyPatch>(patches[patchI]))
 227         {
 228             const cyclicPolyPatch& cycPatch =
 229                 refCast<const cyclicPolyPatch>(patches[patchI]);
 230             checkTransform(cycPatch, applySeparation);
 231
 232             const edgeList& coupledPoints = cycPatch.coupledPoints();
 233             const labelList& meshPts = cycPatch.meshPoints();
 234
 235             // Extract local values. Create map from nbrPoint to value.
 236             Map<T> half0Values(meshPts.size() / 20);
 237             Map<T> half1Values(meshPts.size() / 20);
 238
 239             forAll(coupledPoints, i)
 240             {
 241                 const edge& e = coupledPoints[i];
 242
 243                 typename Map<T>::const_iterator point0Fnd =
 244                     pointValues.find(meshPts[e[0]]);
 245
 246                 if (point0Fnd != pointValues.end())
 247                 {
 248                     half0Values.insert(i, point0Fnd());
 249                 }
 250
 251                 typename Map<T>::const_iterator point1Fnd =
 252                     pointValues.find(meshPts[e[1]]);
 253
 254                 if (point1Fnd != pointValues.end())
 255                 {
 256                     half1Values.insert(i, point1Fnd());
 257                 }
 258             }
 259
 260             if (!cycPatch.parallel())
 261             {
 262                 hasTransformation = true;
 263                 transformList(cycPatch.reverseT(), half0Values);
 264                 transformList(cycPatch.forwardT(), half1Values);
 265             }
 266             else if (applySeparation && cycPatch.separated())
 267             {
 268                 hasTransformation = true;
 269
 270                 const vectorField& v = cycPatch.coupledPolyPatch::separation();
 271                 separateList(v, half0Values);
 272                 separateList(-v, half1Values);
 273             }
 274
 275             forAll(coupledPoints, i)
 276             {
 277                 const edge& e = coupledPoints[i];
 278
 279                 typename Map<T>::const_iterator half1Fnd = half1Values.find(i);
 280
 281                 if (half1Fnd != half1Values.end())
 282                 {
 283                     combine
 284                     (
 285                         pointValues,
 286                         cop,
 287                         meshPts[e[0]],
 288                         half1Fnd()
 289                     );
 290                 }
 291
 292                 typename Map<T>::const_iterator half0Fnd = half0Values.find(i);
 293
 294                 if (half0Fnd != half0Values.end())
 295                 {
 296                     combine
 297                     (
 298                         pointValues,
 299                         cop,
 300                         meshPts[e[1]],
 301                         half0Fnd()
 302                     );
 303                 }
 304             }
 305         }
 306     }
 307
 308     // Note: hasTransformation is only used for warning messages so
 309     // reduction not strictly nessecary.
 310     // reduce(hasTransformation, orOp<bool>());
 311
 312     // Synchronize multiple shared points.
 313     const globalMeshData& pd = mesh.globalData();
 314
 315     if (pd.nGlobalPoints() > 0)
 316     {
 317         if (hasTransformation)
 318         {
 319             WarningIn
 320             (
 321                 "syncTools<class T, class CombineOp>::syncPointMap"
 322                 "(const polyMesh&, Map<T>&, const CombineOp&"
 323                 ", const bool)"
 324             )   << "There are decomposed cyclics in this mesh with"
 325                 << " transformations." << endl
 326                 << "This is not supported. The result will be incorrect"
 327                 << endl;
 328         }
 329         // meshPoint per local index
 330         const labelList& sharedPtLabels = pd.sharedPointLabels();
 331         // global shared index per local index
 332         const labelList& sharedPtAddr = pd.sharedPointAddr();
 333
 334         // Values on shared points. Keyed on global shared index.
 335         Map<T> sharedPointValues(sharedPtAddr.size());
 336
 337
 338         // Fill my entries in the shared points
 339         forAll(sharedPtLabels, i)
 340         {
 341             label meshPointI = sharedPtLabels[i];
 342
 343             typename Map<T>::const_iterator fnd =
 344                 pointValues.find(meshPointI);
 345
 346             if (fnd != pointValues.end())
 347             {
 348                 combine
 349                 (
 350                     sharedPointValues,
 351                     cop,
 352                     sharedPtAddr[i],    // index
 353                     fnd()               // value
 354                 );
 355             }
 356         }
 357
 358         // Reduce on master.
 359
 360         if (Pstream::parRun())
 361         {
 362             if (Pstream::master())
 363             {
 364                 // Receive the edges using shared points from the slave.
 365                 for
 366                 (
 367                     int slave=Pstream::firstSlave();
 368                     slave<=Pstream::lastSlave();
 369                     slave++
 370                 )
 371                 {
 372                     IPstream fromSlave(Pstream::blocking, slave);
 373                     Map<T> nbrValues(fromSlave);
 374
 375                     // Merge neighbouring values with my values
 376                     forAllConstIter(typename Map<T>, nbrValues, iter)
 377                     {
 378                         combine
 379                         (
 380                             sharedPointValues,
 381                             cop,
 382                             iter.key(), // edge
 383                             iter()      // value
 384                         );
 385                     }
 386                 }
 387
 388                 // Send back
 389                 for
 390                 (
 391                     int slave=Pstream::firstSlave();
 392                     slave<=Pstream::lastSlave();
 393                     slave++
 394                 )
 395                 {
 396                     OPstream toSlave(Pstream::blocking, slave);
 397                     toSlave << sharedPointValues;
 398                 }
 399             }
 400             else
 401             {
 402                 // Send to master
 403                 {
 404                     OPstream toMaster
 405                     (
 406                         Pstream::blocking,
 407                         Pstream::masterNo()
 408                     );
 409                     toMaster << sharedPointValues;
 410                 }
 411                 // Receive merged values
 412                 {
 413                     IPstream fromMaster
 414                     (
 415                         Pstream::blocking,
 416                         Pstream::masterNo()
 417                     );
 418                     fromMaster >> sharedPointValues;
 419                 }
 420             }
 421         }
 422
 423
 424         // Merge sharedPointValues (keyed on sharedPointAddr) into
 425         // pointValues (keyed on mesh points).
 426
 427         // Map from global shared index to meshpoint
 428         Map<label> sharedToMeshPoint(2*sharedPtAddr.size());
 429         forAll(sharedPtAddr, i)
 430         {
 431             sharedToMeshPoint.insert(sharedPtAddr[i], sharedPtLabels[i]);
 432         }
 433
 434         forAllConstIter(Map<label>, sharedToMeshPoint, iter)
 435         {
 436             // Do I have a value for my shared point
 437             typename Map<T>::const_iterator sharedFnd =
 438                 sharedPointValues.find(iter.key());
 439
 440             if (sharedFnd != sharedPointValues.end())
 441             {
 442                 combine
 443                 (
 444                     pointValues,
 445                     cop,
 446                     iter(),     // index
 447                     sharedFnd() // value
 448                 );
 449             }
 450         }
 451     }
 452 }
 453
 454
 455 template <class T, class CombineOp>
 456 void Foam::syncTools::syncEdgeMap
 457 (
 458     const polyMesh& mesh,
 459     EdgeMap<T>& edgeValues,
 460     const CombineOp& cop,
 461     const bool applySeparation
 462 )
 463 {
 464     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 465
 466     if (!hasCouples(patches))
 467     {
 468         return;
 469     }
 470
 471
 472     // Do synchronisation without constructing globalEdge addressing
 473     // (since this constructs mesh edge addressing)
 474
 475
 476     // Swap proc patch info
 477     // ~~~~~~~~~~~~~~~~~~~~
 478
 479     if (Pstream::parRun())
 480     {
 481         // Send
 482
 483         forAll(patches, patchI)
 484         {
 485             if
 486             (
 487                 isA<processorPolyPatch>(patches[patchI])
 488              && patches[patchI].nEdges() > 0
 489             )
 490             {
 491                 const processorPolyPatch& procPatch =
 492                     refCast<const processorPolyPatch>(patches[patchI]);
 493
 494                 // Get data per patch edge in neighbouring edge.
 495
 496                 const edgeList& edges = procPatch.edges();
 497                 const labelList& meshPts = procPatch.meshPoints();
 498                 const labelList& nbrPts = procPatch.neighbPoints();
 499
 500                 EdgeMap<T> patchInfo(edges.size() / 20);
 501
 502                 forAll(edges, i)
 503                 {
 504                     const edge& e = edges[i];
 505                     const edge meshEdge(meshPts[e[0]], meshPts[e[1]]);
 506
 507                     typename EdgeMap<T>::const_iterator iter =
 508                         edgeValues.find(meshEdge);
 509
 510                     if (iter != edgeValues.end())
 511                     {
 512                         const edge nbrEdge(nbrPts[e[0]], nbrPts[e[1]]);
 513
 514                         if (nbrEdge[0] >= 0 && nbrEdge[1] >= 0)
 515                         {
 516                             patchInfo.insert(nbrEdge, iter());
 517                         }
 518                     }
 519                 }
 520
 521                 OPstream toNeighb(Pstream::blocking, procPatch.neighbProcNo());
 522                 toNeighb << patchInfo;
 523             }
 524         }
 525
 526
 527         // Receive and combine.
 528
 529         forAll(patches, patchI)
 530         {
 531             if
 532             (
 533                 isA<processorPolyPatch>(patches[patchI])
 534              && patches[patchI].nEdges() > 0
 535             )
 536             {
 537                 const processorPolyPatch& procPatch =
 538                     refCast<const processorPolyPatch>(patches[patchI]);
 539                 checkTransform(procPatch, applySeparation);
 540
 541                 const labelList& meshPts = procPatch.meshPoints();
 542
 543                 IPstream fromNbr(Pstream::blocking, procPatch.neighbProcNo());
 544                 EdgeMap<T> nbrPatchInfo(fromNbr);
 545
 546                 if (!procPatch.parallel())
 547                 {
 548                     transformList(procPatch.forwardT(), nbrPatchInfo);
 549                 }
 550                 else if (applySeparation && procPatch.separated())
 551                 {
 552                     separateList(-procPatch.separation(), nbrPatchInfo);
 553                 }
 554
 555                 // Only update those values which come from neighbour
 556
 557                 forAllConstIter
 558                 (
 559                     typename EdgeMap<T>,
 560                     nbrPatchInfo,
 561                     nbrIter
 562                 )
 563                 {
 564                     const edge& e = nbrIter.key();
 565                     const edge meshEdge(meshPts[e[0]], meshPts[e[1]]);
 566
 567                     combine
 568                     (
 569                         edgeValues,
 570                         cop,
 571                         meshEdge,   // edge
 572                         nbrIter()   // value
 573                     );
 574                 }
 575             }
 576         }
 577     }
 578
 579
 580     // Swap cyclic info
 581     // ~~~~~~~~~~~~~~~~
 582
 583     forAll(patches, patchI)
 584     {
 585         if (isA<cyclicPolyPatch>(patches[patchI]))
 586         {
 587             const cyclicPolyPatch& cycPatch =
 588                 refCast<const cyclicPolyPatch>(patches[patchI]);
 589             checkTransform(cycPatch, applySeparation);
 590
 591             const edgeList& coupledEdges = cycPatch.coupledEdges();
 592             const labelList& meshPts = cycPatch.meshPoints();
 593             const edgeList& edges = cycPatch.edges();
 594
 595             // Extract local values. Create map from nbrPoint to value.
 596             Map<T> half0Values(meshPts.size() / 20);
 597             Map<T> half1Values(meshPts.size() / 20);
 598
 599             forAll(coupledEdges, i)
 600             {
 601                 const edge& twoEdges = coupledEdges[i];
 602
 603                 {
 604                     const edge& e0 = edges[twoEdges[0]];
 605                     const edge meshEdge0(meshPts[e0[0]], meshPts[e0[1]]);
 606
 607                     typename EdgeMap<T>::const_iterator iter =
 608                         edgeValues.find(meshEdge0);
 609
 610                     if (iter != edgeValues.end())
 611                     {
 612                         half0Values.insert(i, iter());
 613                     }
 614                 }
 615                 {
 616                     const edge& e1 = edges[twoEdges[1]];
 617                     const edge meshEdge1(meshPts[e1[0]], meshPts[e1[1]]);
 618
 619                     typename EdgeMap<T>::const_iterator iter =
 620                         edgeValues.find(meshEdge1);
 621
 622                     if (iter != edgeValues.end())
 623                     {
 624                         half1Values.insert(i, iter());
 625                     }
 626                 }
 627             }
 628
 629
 630             // Transform
 631
 632             if (!cycPatch.parallel())
 633             {
 634                 transformList(cycPatch.reverseT(), half0Values);
 635                 transformList(cycPatch.forwardT(), half1Values);
 636             }
 637             else if (applySeparation && cycPatch.separated())
 638             {
 639                 const vectorField& v = cycPatch.coupledPolyPatch::separation();
 640                 separateList(v, half0Values);
 641                 separateList(-v, half1Values);
 642             }
 643
 644
 645             // Extract and combine information
 646
 647             forAll(coupledEdges, i)
 648             {
 649                 const edge& twoEdges = coupledEdges[i];
 650
 651                 typename Map<T>::const_iterator half1Fnd =
 652                     half1Values.find(i);
 653
 654                 if (half1Fnd != half1Values.end())
 655                 {
 656                     const edge& e0 = edges[twoEdges[0]];
 657                     const edge meshEdge0(meshPts[e0[0]], meshPts[e0[1]]);
 658
 659                     combine
 660                     (
 661                         edgeValues,
 662                         cop,
 663                         meshEdge0,  // edge
 664                         half1Fnd()  // value
 665                     );
 666                 }
 667
 668                 typename Map<T>::const_iterator half0Fnd =
 669                     half0Values.find(i);
 670                 if (half0Fnd != half0Values.end())
 671                 {
 672                     const edge& e1 = edges[twoEdges[1]];
 673                     const edge meshEdge1(meshPts[e1[0]], meshPts[e1[1]]);
 674
 675                     combine
 676                     (
 677                         edgeValues,
 678                         cop,
 679                         meshEdge1,  // edge
 680                         half0Fnd()  // value
 681                     );
 682                 }
 683             }
 684         }
 685     }
 686
 687     // Synchronize multiple shared points.
 688     // Problem is that we don't want to construct shared edges so basically
 689     // we do here like globalMeshData but then using sparse edge representation
 690     // (EdgeMap instead of mesh.edges())
 691
 692     const globalMeshData& pd = mesh.globalData();
 693     const labelList& sharedPtAddr = pd.sharedPointAddr();
 694     const labelList& sharedPtLabels = pd.sharedPointLabels();
 695
 696     // 1. Create map from meshPoint to globalShared index.
 697     Map<label> meshToShared(2*sharedPtLabels.size());
 698     forAll(sharedPtLabels, i)
 699     {
 700         meshToShared.insert(sharedPtLabels[i], sharedPtAddr[i]);
 701     }
 702
 703     // Values on shared points. From two sharedPtAddr indices to a value.
 704     EdgeMap<T> sharedEdgeValues(meshToShared.size());
 705
 706     // From shared edge to mesh edge. Used for merging later on.
 707     EdgeMap<edge> potentialSharedEdge(meshToShared.size());
 708
 709     // 2. Find any edges using two global shared points. These will always be
 710     // on the outside of the mesh. (though might not be on coupled patch
 711     // if is single edge and not on coupled face)
 712     // Store value (if any) on sharedEdgeValues
 713     for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
 714     {
 715         const face& f = mesh.faces()[faceI];
 716
 717         forAll(f, fp)
 718         {
 719             label v0 = f[fp];
 720             label v1 = f[f.fcIndex(fp)];
 721
 722             Map<label>::const_iterator v0Fnd = meshToShared.find(v0);
 723
 724             if (v0Fnd != meshToShared.end())
 725             {
 726                 Map<label>::const_iterator v1Fnd = meshToShared.find(v1);
 727
 728                 if (v1Fnd != meshToShared.end())
 729                 {
 730                     const edge meshEdge(v0, v1);
 731
 732                     // edge in shared point labels
 733                     const edge sharedEdge(v0Fnd(), v1Fnd());
 734
 735                     // Store mesh edge as a potential shared edge.
 736                     potentialSharedEdge.insert(sharedEdge, meshEdge);
 737
 738                     typename EdgeMap<T>::const_iterator edgeFnd =
 739                         edgeValues.find(meshEdge);
 740
 741                     if (edgeFnd != edgeValues.end())
 742                     {
 743                         // edge exists in edgeValues. See if already in map
 744                         // (since on same processor, e.g. cyclic)
 745                         combine
 746                         (
 747                             sharedEdgeValues,
 748                             cop,
 749                             sharedEdge, // edge
 750                             edgeFnd()   // value
 751                         );
 752                     }
 753                 }
 754             }
 755         }
 756     }
 757
 758
 759     // Now sharedEdgeValues will contain per potential sharedEdge the value.
 760     // (potential since an edge having two shared points is not nessecary a
 761     //  shared edge).
 762     // Reduce this on the master.
 763
 764     if (Pstream::parRun())
 765     {
 766         if (Pstream::master())
 767         {
 768             // Receive the edges using shared points from the slave.
 769             for
 770             (
 771                 int slave=Pstream::firstSlave();
 772                 slave<=Pstream::lastSlave();
 773                 slave++
 774             )
 775             {
 776                 IPstream fromSlave(Pstream::blocking, slave);
 777                 EdgeMap<T> nbrValues(fromSlave);
 778
 779                 // Merge neighbouring values with my values
 780                 forAllConstIter(typename EdgeMap<T>, nbrValues, iter)
 781                 {
 782                     combine
 783                     (
 784                         sharedEdgeValues,
 785                         cop,
 786                         iter.key(), // edge
 787                         iter()      // value
 788                     );
 789                 }
 790             }
 791
 792             // Send back
 793             for
 794             (
 795                 int slave=Pstream::firstSlave();
 796                 slave<=Pstream::lastSlave();
 797                 slave++
 798             )
 799             {
 800
 801                 OPstream toSlave(Pstream::blocking, slave);
 802                 toSlave << sharedEdgeValues;
 803             }
 804         }
 805         else
 806         {
 807             // Send to master
 808             {
 809                 OPstream toMaster(Pstream::blocking, Pstream::masterNo());
 810                 toMaster << sharedEdgeValues;
 811             }
 812             // Receive merged values
 813             {
 814                 IPstream fromMaster(Pstream::blocking, Pstream::masterNo());
 815                 fromMaster >> sharedEdgeValues;
 816             }
 817         }
 818     }
 819
 820
 821     // Merge sharedEdgeValues (keyed on sharedPointAddr) into edgeValues
 822     // (keyed on mesh points).
 823
 824     // Loop over all my shared edges.
 825     forAllConstIter(typename EdgeMap<edge>, potentialSharedEdge, iter)
 826     {
 827         const edge& sharedEdge = iter.key();
 828         const edge& meshEdge = iter();
 829
 830         // Do I have a value for the shared edge?
 831         typename EdgeMap<T>::const_iterator sharedFnd =
 832             sharedEdgeValues.find(sharedEdge);
 833
 834         if (sharedFnd != sharedEdgeValues.end())
 835         {
 836             combine
 837             (
 838                 edgeValues,
 839                 cop,
 840                 meshEdge,       // edge
 841                 sharedFnd()     // value
 842             );
 843         }
 844     }
 845 }
 846
 847
 848 template <class T, class CombineOp>
 849 void Foam::syncTools::syncPointList
 850 (
 851     const polyMesh& mesh,
 852     UList<T>& pointValues,
 853     const CombineOp& cop,
 854     const T& nullValue,
 855     const bool applySeparation
 856 )
 857 {
 858     if (pointValues.size() != mesh.nPoints())
 859     {
 860         FatalErrorIn
 861         (
 862             "syncTools<class T, class CombineOp>::syncPointList"
 863             "(const polyMesh&, UList<T>&, const CombineOp&, const T&"
 864             ", const bool)"
 865         )   << "Number of values " << pointValues.size()
 866             << " is not equal to the number of points in the mesh "
 867             << mesh.nPoints() << abort(FatalError);
 868     }
 869
 870     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 871
 872     if (!hasCouples(patches))
 873     {
 874         return;
 875     }
 876
 877     // Is there any coupled patch with transformation?
 878     bool hasTransformation = false;
 879
 880     if (Pstream::parRun())
 881     {
 882         // Send
 883
 884         forAll(patches, patchI)
 885         {
 886             if
 887             (
 888                 isA<processorPolyPatch>(patches[patchI])
 889              && patches[patchI].nPoints() > 0
 890             )
 891             {
 892                 const processorPolyPatch& procPatch =
 893                     refCast<const processorPolyPatch>(patches[patchI]);
 894
 895                 // Get data per patchPoint in neighbouring point numbers.
 896                 List<T> patchInfo(procPatch.nPoints(), nullValue);
 897
 898                 const labelList& meshPts = procPatch.meshPoints();
 899                 const labelList& nbrPts = procPatch.neighbPoints();
 900
 901                 forAll(nbrPts, pointI)
 902                 {
 903                     label nbrPointI = nbrPts[pointI];
 904                     if (nbrPointI >= 0 && nbrPointI < patchInfo.size())
 905                     {
 906                         patchInfo[nbrPointI] = pointValues[meshPts[pointI]];
 907                     }
 908                 }
 909
 910                 OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
 911                 toNbr << patchInfo;
 912             }
 913         }
 914
 915
 916         // Receive and combine.
 917
 918         forAll(patches, patchI)
 919         {
 920             if
 921             (
 922                 isA<processorPolyPatch>(patches[patchI])
 923              && patches[patchI].nPoints() > 0
 924             )
 925             {
 926                 const processorPolyPatch& procPatch =
 927                     refCast<const processorPolyPatch>(patches[patchI]);
 928                 checkTransform(procPatch, applySeparation);
 929
 930                 List<T> nbrPatchInfo(procPatch.nPoints());
 931                 {
 932                     // We do not know the number of points on the other side
 933                     // so cannot use Pstream::read.
 934                     IPstream fromNbr
 935                     (
 936                         Pstream::blocking,
 937                         procPatch.neighbProcNo()
 938                     );
 939                     fromNbr >> nbrPatchInfo;
 940                 }
 941                 // Null any value which is not on neighbouring processor
 942                 nbrPatchInfo.setSize(procPatch.nPoints(), nullValue);
 943
 944                 if (!procPatch.parallel())
 945                 {
 946                     hasTransformation = true;
 947                     transformList(procPatch.forwardT(), nbrPatchInfo);
 948                 }
 949                 else if (applySeparation && procPatch.separated())
 950                 {
 951                     hasTransformation = true;
 952                     separateList(-procPatch.separation(), nbrPatchInfo);
 953                 }
 954
 955                 const labelList& meshPts = procPatch.meshPoints();
 956
 957                 forAll(meshPts, pointI)
 958                 {
 959                     label meshPointI = meshPts[pointI];
 960                     cop(pointValues[meshPointI], nbrPatchInfo[pointI]);
 961                 }
 962             }
 963         }
 964     }
 965
 966     // Do the cyclics.
 967     forAll(patches, patchI)
 968     {
 969         if (isA<cyclicPolyPatch>(patches[patchI]))
 970         {
 971             const cyclicPolyPatch& cycPatch =
 972                 refCast<const cyclicPolyPatch>(patches[patchI]);
 973
 974             checkTransform(cycPatch, applySeparation);
 975
 976             const edgeList& coupledPoints = cycPatch.coupledPoints();
 977             const labelList& meshPts = cycPatch.meshPoints();
 978
 979             List<T> half0Values(coupledPoints.size());
 980             List<T> half1Values(coupledPoints.size());
 981
 982             forAll(coupledPoints, i)
 983             {
 984                 const edge& e = coupledPoints[i];
 985
 986                 label point0 = meshPts[e[0]];
 987                 label point1 = meshPts[e[1]];
 988
 989                 half0Values[i] = pointValues[point0];
 990                 half1Values[i] = pointValues[point1];
 991             }
 992
 993             if (!cycPatch.parallel())
 994             {
 995                 hasTransformation = true;
 996                 transformList(cycPatch.reverseT(), half0Values);
 997                 transformList(cycPatch.forwardT(), half1Values);
 998             }
 999             else if (applySeparation && cycPatch.separated())
1000             {
1001                 hasTransformation = true;
1002                 const vectorField& v = cycPatch.coupledPolyPatch::separation();
1003                 separateList(v, half0Values);
1004                 separateList(-v, half1Values);
1005             }
1006
1007             forAll(coupledPoints, i)
1008             {
1009                 const edge& e = coupledPoints[i];
1010
1011                 label point0 = meshPts[e[0]];
1012                 label point1 = meshPts[e[1]];
1013
1014                 cop(pointValues[point0], half1Values[i]);
1015                 cop(pointValues[point1], half0Values[i]);
1016             }
1017         }
1018     }
1019
1020     //- Note: hasTransformation is only used for warning messages so
1021     //  reduction not strictly nessecary.
1022     //reduce(hasTransformation, orOp<bool>());
1023
1024     // Synchronize multiple shared points.
1025     const globalMeshData& pd = mesh.globalData();
1026
1027     if (pd.nGlobalPoints() > 0)
1028     {
1029         if (hasTransformation)
1030         {
1031             WarningIn
1032             (
1033                 "syncTools<class T, class CombineOp>::syncPointList"
1034                 "(const polyMesh&, UList<T>&, const CombineOp&, const T&"
1035                 ", const bool)"
1036             )   << "There are decomposed cyclics in this mesh with"
1037                 << " transformations." << endl
1038                 << "This is not supported. The result will be incorrect"
1039                 << endl;
1040         }
1041
1042
1043         // Values on shared points.
1044         List<T> sharedPts(pd.nGlobalPoints(), nullValue);
1045
1046         forAll(pd.sharedPointLabels(), i)
1047         {
1048             label meshPointI = pd.sharedPointLabels()[i];
1049             // Fill my entries in the shared points
1050             sharedPts[pd.sharedPointAddr()[i]] = pointValues[meshPointI];
1051         }
1052
1053         // Combine on master.
1054         Pstream::listCombineGather(sharedPts, cop);
1055         Pstream::listCombineScatter(sharedPts);
1056
1057         // Now we will all have the same information. Merge it back with
1058         // my local information.
1059         forAll(pd.sharedPointLabels(), i)
1060         {
1061             label meshPointI = pd.sharedPointLabels()[i];
1062             pointValues[meshPointI] = sharedPts[pd.sharedPointAddr()[i]];
1063         }
1064     }
1065 }
1066
1067
1068 template <class T, class CombineOp>
1069 void Foam::syncTools::syncPointList
1070 (
1071     const polyMesh& mesh,
1072     const labelList& meshPoints,
1073     UList<T>& pointValues,
1074     const CombineOp& cop,
1075     const T& nullValue,
1076     const bool applySeparation
1077 )
1078 {
1079     if (pointValues.size() != meshPoints.size())
1080     {
1081         FatalErrorIn
1082         (
1083             "syncTools<class T, class CombineOp>::syncPointList"
1084             "(const polyMesh&, const labelList&, UList<T>&, const CombineOp&"
1085             ", const T&, const bool)"
1086         )   << "Number of values " << pointValues.size()
1087             << " is not equal to the number of points "
1088             << meshPoints.size() << abort(FatalError);
1089     }
1090
1091     if (!hasCouples(mesh.boundaryMesh()))
1092     {
1093         return;
1094     }
1095
1096     List<T> meshValues(mesh.nPoints(), nullValue);
1097
1098     forAll(meshPoints, i)
1099     {
1100         meshValues[meshPoints[i]] = pointValues[i];
1101     }
1102
1103     syncTools::syncPointList
1104     (
1105         mesh,
1106         meshValues,
1107         cop,            // combine op
1108         nullValue,      // null value
1109         applySeparation // separation
1110     );
1111
1112     forAll(meshPoints, i)
1113     {
1114         pointValues[i] = meshValues[meshPoints[i]];
1115     }
1116 }
1117
1118
1119 template <class T, class CombineOp>
1120 void Foam::syncTools::syncEdgeList
1121 (
1122     const polyMesh& mesh,
1123     UList<T>& edgeValues,
1124     const CombineOp& cop,
1125     const T& nullValue,
1126     const bool applySeparation
1127 )
1128 {
1129     if (edgeValues.size() != mesh.nEdges())
1130     {
1131         FatalErrorIn
1132         (
1133             "syncTools<class T, class CombineOp>::syncEdgeList"
1134             "(const polyMesh&, UList<T>&, const CombineOp&, const T&"
1135             ", const bool)"
1136         )   << "Number of values " << edgeValues.size()
1137             << " is not equal to the number of edges in the mesh "
1138             << mesh.nEdges() << abort(FatalError);
1139     }
1140
1141     const polyBoundaryMesh& patches = mesh.boundaryMesh();
1142
1143     if (!hasCouples(patches))
1144     {
1145         return;
1146     }
1147
1148     // Is there any coupled patch with transformation?
1149     bool hasTransformation = false;
1150
1151     if (Pstream::parRun())
1152     {
1153         // Send
1154
1155         forAll(patches, patchI)
1156         {
1157             if
1158             (
1159                 isA<processorPolyPatch>(patches[patchI])
1160              && patches[patchI].nEdges() > 0
1161             )
1162             {
1163                 const processorPolyPatch& procPatch =
1164                     refCast<const processorPolyPatch>(patches[patchI]);
1165
1166                 const labelList& meshEdges = procPatch.meshEdges();
1167                 const labelList& neighbEdges = procPatch.neighbEdges();
1168
1169                 // Get region per patch edge in neighbouring edge numbers.
1170                 List<T> patchInfo(procPatch.nEdges(), nullValue);
1171
1172                 forAll(neighbEdges, edgeI)
1173                 {
1174                     label nbrEdgeI = neighbEdges[edgeI];
1175
1176                     if (nbrEdgeI >= 0 && nbrEdgeI < patchInfo.size())
1177                     {
1178                         patchInfo[nbrEdgeI] = edgeValues[meshEdges[edgeI]];
1179                     }
1180                 }
1181
1182                 OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
1183                 toNbr << patchInfo;
1184             }
1185         }
1186
1187         // Receive and combine.
1188
1189         forAll(patches, patchI)
1190         {
1191             if
1192             (
1193                 isA<processorPolyPatch>(patches[patchI])
1194              && patches[patchI].nEdges() > 0
1195             )
1196             {
1197                 const processorPolyPatch& procPatch =
1198                     refCast<const processorPolyPatch>(patches[patchI]);
1199
1200                 checkTransform(procPatch, applySeparation);
1201
1202                 const labelList& meshEdges = procPatch.meshEdges();
1203
1204                 // Receive from neighbour. Is per patch edge the region of the
1205                 // neighbouring patch edge.
1206                 List<T> nbrPatchInfo(procPatch.nEdges());
1207
1208                 {
1209                     IPstream fromNeighb
1210                     (
1211                         Pstream::blocking,
1212                         procPatch.neighbProcNo()
1213                     );
1214                     fromNeighb >> nbrPatchInfo;
1215                 }
1216                 // Null any value which is not on neighbouring processor
1217                 nbrPatchInfo.setSize(procPatch.nEdges(), nullValue);
1218
1219                 if (!procPatch.parallel())
1220                 {
1221                     hasTransformation = true;
1222                     transformList(procPatch.forwardT(), nbrPatchInfo);
1223                 }
1224                 else if (applySeparation && procPatch.separated())
1225                 {
1226                     hasTransformation = true;
1227                     separateList(-procPatch.separation(), nbrPatchInfo);
1228                 }
1229
1230                 forAll(meshEdges, edgeI)
1231                 {
1232                     label meshEdgeI = meshEdges[edgeI];
1233
1234                     cop(edgeValues[meshEdgeI], nbrPatchInfo[edgeI]);
1235                 }
1236             }
1237         }
1238     }
1239
1240     // Do the cyclics.
1241     forAll(patches, patchI)
1242     {
1243         if (isA<cyclicPolyPatch>(patches[patchI]))
1244         {
1245             const cyclicPolyPatch& cycPatch =
1246                 refCast<const cyclicPolyPatch>(patches[patchI]);
1247
1248             checkTransform(cycPatch, applySeparation);
1249
1250             const edgeList& coupledEdges = cycPatch.coupledEdges();
1251             const labelList& meshEdges = cycPatch.meshEdges();
1252
1253             List<T> half0Values(coupledEdges.size());
1254             List<T> half1Values(coupledEdges.size());
1255
1256             forAll(coupledEdges, i)
1257             {
1258                 const edge& e = coupledEdges[i];
1259
1260                 label meshEdge0 = meshEdges[e[0]];
1261                 label meshEdge1 = meshEdges[e[1]];
1262
1263                 half0Values[i] = edgeValues[meshEdge0];
1264                 half1Values[i] = edgeValues[meshEdge1];
1265             }
1266
1267             if (!cycPatch.parallel())
1268             {
1269                 hasTransformation = true;
1270                 transformList(cycPatch.reverseT(), half0Values);
1271                 transformList(cycPatch.forwardT(), half1Values);
1272             }
1273             else if (applySeparation && cycPatch.separated())
1274             {
1275                 hasTransformation = true;
1276
1277                 const vectorField& v = cycPatch.coupledPolyPatch::separation();
1278                 separateList(v, half0Values);
1279                 separateList(-v, half1Values);
1280             }
1281
1282             forAll(coupledEdges, i)
1283             {
1284                 const edge& e = coupledEdges[i];
1285
1286                 label meshEdge0 = meshEdges[e[0]];
1287                 label meshEdge1 = meshEdges[e[1]];
1288
1289                 cop(edgeValues[meshEdge0], half1Values[i]);
1290                 cop(edgeValues[meshEdge1], half0Values[i]);
1291             }
1292         }
1293     }
1294
1295     //- Note: hasTransformation is only used for warning messages so
1296     //  reduction not strictly nessecary.
1297     //reduce(hasTransformation, orOp<bool>());
1298
1299     // Do the multiple shared edges
1300     const globalMeshData& pd = mesh.globalData();
1301
1302     if (pd.nGlobalEdges() > 0)
1303     {
1304         if (hasTransformation)
1305         {
1306             WarningIn
1307             (
1308                 "syncTools<class T, class CombineOp>::syncEdgeList"
1309                 "(const polyMesh&, UList<T>&, const CombineOp&, const T&"
1310                 ", const bool)"
1311             )   << "There are decomposed cyclics in this mesh with"
1312                 << " transformations." << endl
1313                 << "This is not supported. The result will be incorrect"
1314                 << endl;
1315         }
1316
1317         // Values on shared edges.
1318         List<T> sharedPts(pd.nGlobalEdges(), nullValue);
1319
1320         forAll(pd.sharedEdgeLabels(), i)
1321         {
1322             label meshEdgeI = pd.sharedEdgeLabels()[i];
1323
1324             // Fill my entries in the shared edges
1325             sharedPts[pd.sharedEdgeAddr()[i]] = edgeValues[meshEdgeI];
1326         }
1327
1328         // Combine on master.
1329         Pstream::listCombineGather(sharedPts, cop);
1330         Pstream::listCombineScatter(sharedPts);
1331
1332         // Now we will all have the same information. Merge it back with
1333         // my local information.
1334         forAll(pd.sharedEdgeLabels(), i)
1335         {
1336             label meshEdgeI = pd.sharedEdgeLabels()[i];
1337             edgeValues[meshEdgeI] = sharedPts[pd.sharedEdgeAddr()[i]];
1338         }
1339     }
1340 }
1341
1342
1343 template <class T, class CombineOp>
1344 void Foam::syncTools::syncBoundaryFaceList
1345 (
1346     const polyMesh& mesh,
1347     UList<T>& faceValues,
1348     const CombineOp& cop,
1349     const bool applySeparation
1350 )
1351 {
1352     const label nBFaces = mesh.nFaces() - mesh.nInternalFaces();
1353
1354     if (faceValues.size() != nBFaces)
1355     {
1356         FatalErrorIn
1357         (
1358             "syncTools<class T, class CombineOp>::syncBoundaryFaceList"
1359             "(const polyMesh&, UList<T>&, const CombineOp&"
1360             ", const bool)"
1361         )   << "Number of values " << faceValues.size()
1362             << " is not equal to the number of boundary faces in the mesh "
1363             << nBFaces << abort(FatalError);
1364     }
1365
1366     const polyBoundaryMesh& patches = mesh.boundaryMesh();
1367
1368     if (!hasCouples(patches))
1369     {
1370         return;
1371     }
1372
1373
1374     if (Pstream::parRun())
1375     {
1376         // Send
1377
1378         forAll(patches, patchI)
1379         {
1380             if
1381             (
1382                 isA<processorPolyPatch>(patches[patchI])
1383              && patches[patchI].size() > 0
1384             )
1385             {
1386                 const processorPolyPatch& procPatch =
1387                     refCast<const processorPolyPatch>(patches[patchI]);
1388
1389                 label patchStart = procPatch.start()-mesh.nInternalFaces();
1390
1391                 if (contiguous<T>())
1392                 {
1393                     OPstream::write
1394                     (
1395                         Pstream::blocking,
1396                         procPatch.neighbProcNo(),
1397                         reinterpret_cast<const char*>(&faceValues[patchStart]),
1398                         procPatch.size()*sizeof(T)
1399                     );
1400                 }
1401                 else
1402                 {
1403                     OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
1404                     toNbr <<
1405                         SubList<T>(faceValues, procPatch.size(), patchStart);
1406                 }
1407             }
1408         }
1409
1410
1411         // Receive and combine.
1412
1413         forAll(patches, patchI)
1414         {
1415             if
1416             (
1417                 isA<processorPolyPatch>(patches[patchI])
1418              && patches[patchI].size() > 0
1419             )
1420             {
1421                 const processorPolyPatch& procPatch =
1422                     refCast<const processorPolyPatch>(patches[patchI]);
1423
1424                 List<T> nbrPatchInfo(procPatch.size());
1425
1426                 if (contiguous<T>())
1427                 {
1428                     IPstream::read
1429                     (
1430                         Pstream::blocking,
1431                         procPatch.neighbProcNo(),
1432                         reinterpret_cast<char*>(nbrPatchInfo.begin()),
1433                         nbrPatchInfo.byteSize()
1434                     );
1435                 }
1436                 else
1437                 {
1438                     IPstream fromNeighb
1439                     (
1440                         Pstream::blocking,
1441                         procPatch.neighbProcNo()
1442                     );
1443                     fromNeighb >> nbrPatchInfo;
1444                 }
1445
1446                 if (!procPatch.parallel())
1447                 {
1448                     transformList(procPatch.forwardT(), nbrPatchInfo);
1449                 }
1450                 else if (applySeparation && procPatch.separated())
1451                 {
1452                     separateList(-procPatch.separation(), nbrPatchInfo);
1453                 }
1454
1455
1456                 label bFaceI = procPatch.start()-mesh.nInternalFaces();
1457
1458                 forAll(nbrPatchInfo, i)
1459                 {
1460                     cop(faceValues[bFaceI++], nbrPatchInfo[i]);
1461                 }
1462             }
1463         }
1464     }
1465
1466     // Do the cyclics.
1467     forAll(patches, patchI)
1468     {
1469         if (isA<cyclicPolyPatch>(patches[patchI]))
1470         {
1471             const cyclicPolyPatch& cycPatch =
1472                 refCast<const cyclicPolyPatch>(patches[patchI]);
1473
1474             label patchStart = cycPatch.start()-mesh.nInternalFaces();
1475
1476             label half = cycPatch.size()/2;
1477             label half1Start = patchStart+half;
1478
1479             List<T> half0Values(SubList<T>(faceValues, half, patchStart));
1480             List<T> half1Values(SubList<T>(faceValues, half, half1Start));
1481
1482             if (!cycPatch.parallel())
1483             {
1484                 transformList(cycPatch.reverseT(), half0Values);
1485                 transformList(cycPatch.forwardT(), half1Values);
1486             }
1487             else if (applySeparation && cycPatch.separated())
1488             {
1489                 const vectorField& v = cycPatch.coupledPolyPatch::separation();
1490                 separateList(v, half0Values);
1491                 separateList(-v, half1Values);
1492             }
1493
1494             label i0 = patchStart;
1495             forAll(half1Values, i)
1496             {
1497                 cop(faceValues[i0++], half1Values[i]);
1498             }
1499
1500             label i1 = half1Start;
1501             forAll(half0Values, i)
1502             {
1503                 cop(faceValues[i1++], half0Values[i]);
1504             }
1505         }
1506     }
1507 }
1508
1509
1510 template <class T, class CombineOp>
1511 void Foam::syncTools::syncFaceList
1512 (
1513     const polyMesh& mesh,
1514     UList<T>& faceValues,
1515     const CombineOp& cop,
1516     const bool applySeparation
1517 )
1518 {
1519     if (faceValues.size() != mesh.nFaces())
1520     {
1521         FatalErrorIn
1522         (
1523             "syncTools<class T, class CombineOp>::syncFaceList"
1524             "(const polyMesh&, UList<T>&, const CombineOp&"
1525             ", const bool)"
1526         )   << "Number of values " << faceValues.size()
1527             << " is not equal to the number of faces in the mesh "
1528             << mesh.nFaces() << abort(FatalError);
1529     }
1530
1531     SubList<T> bndValues
1532     (
1533         faceValues,
1534         mesh.nFaces()-mesh.nInternalFaces(),
1535         mesh.nInternalFaces()
1536     );
1537
1538     syncBoundaryFaceList
1539     (
1540         mesh,
1541         bndValues,
1542         cop,
1543         applySeparation
1544     );
1545 }
1546
1547
1548 template <class T>
1549 void Foam::syncTools::swapBoundaryFaceList
1550 (
1551     const polyMesh& mesh,
1552     UList<T>& faceValues,
1553     const bool applySeparation
1554 )
1555 {
1556     syncBoundaryFaceList(mesh, faceValues, eqOp<T>(), applySeparation);
1557 }
1558
1559
1560 template <class T>
1561 void Foam::syncTools::swapFaceList
1562 (
1563     const polyMesh& mesh,
1564     UList<T>& faceValues,
1565     const bool applySeparation
1566 )
1567 {
1568     syncFaceList(mesh, faceValues, eqOp<T>(), applySeparation);
1569 }
1570
1571
1572 template <unsigned nBits, class CombineOp>
1573 void Foam::syncTools::syncFaceList
1574 (
1575     const polyMesh& mesh,
1576     PackedList<nBits>& faceValues,
1577     const CombineOp& cop
1578 )
1579 {
1580     if (faceValues.size() != mesh.nFaces())
1581     {
1582         FatalErrorIn
1583         (
1584             "syncTools<unsigned nBits, class CombineOp>::syncFaceList"
1585             "(const polyMesh&, PackedList<nBits>&, const CombineOp&)"
1586         )   << "Number of values " << faceValues.size()
1587             << " is not equal to the number of faces in the mesh "
1588             << mesh.nFaces() << abort(FatalError);
1589     }
1590
1591     const polyBoundaryMesh& patches = mesh.boundaryMesh();
1592
1593     if (!hasCouples(patches))
1594     {
1595         return;
1596     }
1597
1598     // Patch data (proc patches only).
1599     List<List<unsigned int> > patchValues(patches.size());
1600
1601     if (Pstream::parRun())
1602     {
1603         // Send
1604
1605         forAll(patches, patchI)
1606         {
1607             if
1608             (
1609                 isA<processorPolyPatch>(patches[patchI])
1610              && patches[patchI].size() > 0
1611             )
1612             {
1613                 const processorPolyPatch& procPatch =
1614                     refCast<const processorPolyPatch>(patches[patchI]);
1615
1616                 patchValues[patchI].setSize(procPatch.size());
1617                 forAll(procPatch, i)
1618                 {
1619                     patchValues[patchI][i] =
1620                         faceValues.get(procPatch.start()+i);
1621                 }
1622
1623                 OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
1624                 toNbr << patchValues[patchI];
1625             }
1626         }
1627
1628
1629         // Receive and combine.
1630
1631         forAll(patches, patchI)
1632         {
1633             if
1634             (
1635                 isA<processorPolyPatch>(patches[patchI])
1636              && patches[patchI].size() > 0
1637             )
1638             {
1639                 const processorPolyPatch& procPatch =
1640                     refCast<const processorPolyPatch>(patches[patchI]);
1641
1642                 {
1643                     IPstream fromNbr
1644                     (
1645                         Pstream::blocking,
1646                         procPatch.neighbProcNo()
1647                     );
1648                     fromNbr >> patchValues[patchI];
1649                 }
1650
1651                 // Combine (bitwise)
1652                 forAll(procPatch, i)
1653                 {
1654                     unsigned int patchVal = patchValues[patchI][i];
1655                     label meshFaceI = procPatch.start()+i;
1656                     unsigned int faceVal = faceValues.get(meshFaceI);
1657                     cop(faceVal, patchVal);
1658                     faceValues.set(meshFaceI, faceVal);
1659                 }
1660             }
1661         }
1662     }
1663
1664     // Do the cyclics.
1665     forAll(patches, patchI)
1666     {
1667         if (isA<cyclicPolyPatch>(patches[patchI]))
1668         {
1669             const cyclicPolyPatch& cycPatch =
1670                 refCast<const cyclicPolyPatch>(patches[patchI]);
1671
1672             label half = cycPatch.size()/2;
1673
1674             for (label i = 0; i < half; i++)
1675             {
1676                 label meshFace0 = cycPatch.start()+i;
1677                 unsigned int val0 = faceValues.get(meshFace0);
1678                 label meshFace1 = meshFace0 + half;
1679                 unsigned int val1 = faceValues.get(meshFace1);
1680
1681                 unsigned int t = val0;
1682                 cop(t, val1);
1683                 faceValues.set(meshFace0, t);
1684
1685                 cop(val1, val0);
1686                 faceValues.set(meshFace1, val1);
1687             }
1688         }
1689     }
1690 }
1691
1692
1693 template <unsigned nBits>
1694 void Foam::syncTools::swapFaceList
1695 (
1696     const polyMesh& mesh,
1697     PackedList<nBits>& faceValues
1698 )
1699 {
1700     syncFaceList(mesh, faceValues, eqOp<unsigned int>());
1701 }
1702
1703
1704 template <unsigned nBits, class CombineOp>
1705 void Foam::syncTools::syncPointList
1706 (
1707     const polyMesh& mesh,
1708     PackedList<nBits>& pointValues,
1709     const CombineOp& cop,
1710     const unsigned int nullValue
1711 )
1712 {
1713     if (pointValues.size() != mesh.nPoints())
1714     {
1715         FatalErrorIn
1716         (
1717             "syncTools<unsigned nBits, class CombineOp>::syncPointList"
1718             "(const polyMesh&, PackedList<nBits>&, const CombineOp&"
1719             ", const unsigned int&)"
1720         )   << "Number of values " << pointValues.size()
1721             << " is not equal to the number of points in the mesh "
1722             << mesh.nPoints() << abort(FatalError);
1723     }
1724
1725     const polyBoundaryMesh& patches = mesh.boundaryMesh();
1726
1727     if (!hasCouples(patches))
1728     {
1729         return;
1730     }
1731
1732     // Patch data (proc patches only).
1733     List<List<unsigned int> > patchValues(patches.size());
1734
1735     if (Pstream::parRun())
1736     {
1737         // Send
1738
1739         forAll(patches, patchI)
1740         {
1741             if
1742             (
1743                 isA<processorPolyPatch>(patches[patchI])
1744              && patches[patchI].nPoints() > 0
1745             )
1746             {
1747                 const processorPolyPatch& procPatch =
1748                     refCast<const processorPolyPatch>(patches[patchI]);
1749
1750                 patchValues[patchI].setSize(procPatch.nPoints());
1751                 patchValues[patchI] = nullValue;
1752
1753                 const labelList& meshPts = procPatch.meshPoints();
1754                 const labelList& nbrPts = procPatch.neighbPoints();
1755
1756                 forAll(nbrPts, pointI)
1757                 {
1758                     label nbrPointI = nbrPts[pointI];
1759                     if (nbrPointI >= 0 && nbrPointI < procPatch.nPoints())
1760                     {
1761                         patchValues[patchI][nbrPointI] =
1762                             pointValues.get(meshPts[pointI]);
1763                     }
1764                 }
1765
1766                 OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
1767                 toNbr << patchValues[patchI];
1768             }
1769         }
1770
1771
1772         // Receive and combine.
1773
1774         forAll(patches, patchI)
1775         {
1776             if
1777             (
1778                 isA<processorPolyPatch>(patches[patchI])
1779              && patches[patchI].nPoints() > 0
1780             )
1781             {
1782                 const processorPolyPatch& procPatch =
1783                     refCast<const processorPolyPatch>(patches[patchI]);
1784
1785                 {
1786                     // We do not know the number of points on the other side
1787                     // so cannot use Pstream::read.
1788                     IPstream fromNbr
1789                     (
1790                         Pstream::blocking,
1791                         procPatch.neighbProcNo()
1792                     );
1793                     fromNbr >> patchValues[patchI];
1794                 }
1795
1796                 // Null any value which is not on neighbouring processor
1797                 patchValues[patchI].setSize(procPatch.nPoints(), nullValue);
1798
1799                 const labelList& meshPts = procPatch.meshPoints();
1800
1801                 forAll(meshPts, pointI)
1802                 {
1803                     label meshPointI = meshPts[pointI];
1804                     unsigned int pointVal = pointValues.get(meshPointI);
1805                     cop(pointVal, patchValues[patchI][pointI]);
1806                     pointValues.set(meshPointI, pointVal);
1807                 }
1808             }
1809         }
1810     }
1811
1812     // Do the cyclics.
1813     forAll(patches, patchI)
1814     {
1815         if (isA<cyclicPolyPatch>(patches[patchI]))
1816         {
1817             const cyclicPolyPatch& cycPatch =
1818                 refCast<const cyclicPolyPatch>(patches[patchI]);
1819
1820             const edgeList& coupledPoints = cycPatch.coupledPoints();
1821             const labelList& meshPts = cycPatch.meshPoints();
1822
1823             forAll(coupledPoints, i)
1824             {
1825                 const edge& e = coupledPoints[i];
1826
1827                 label point0 = meshPts[e[0]];
1828                 label point1 = meshPts[e[1]];
1829
1830                 unsigned int val0 = pointValues.get(point0);
1831                 unsigned int t = val0;
1832                 unsigned int val1 = pointValues.get(point1);
1833
1834                 cop(t, val1);
1835                 pointValues.set(point0, t);
1836                 cop(val1, val0);
1837                 pointValues.set(point1, val1);
1838             }
1839         }
1840     }
1841
1842     // Synchronize multiple shared points.
1843     const globalMeshData& pd = mesh.globalData();
1844
1845     if (pd.nGlobalPoints() > 0)
1846     {
1847         // Values on shared points. Use unpacked storage for ease!
1848         List<unsigned int> sharedPts(pd.nGlobalPoints(), nullValue);
1849
1850         forAll(pd.sharedPointLabels(), i)
1851         {
1852             label meshPointI = pd.sharedPointLabels()[i];
1853             // Fill my entries in the shared points
1854             sharedPts[pd.sharedPointAddr()[i]] = pointValues.get(meshPointI);
1855         }
1856
1857         // Combine on master.
1858         Pstream::listCombineGather(sharedPts, cop);
1859         Pstream::listCombineScatter(sharedPts);
1860
1861         // Now we will all have the same information. Merge it back with
1862         // my local information.
1863         forAll(pd.sharedPointLabels(), i)
1864         {
1865             label meshPointI = pd.sharedPointLabels()[i];
1866             pointValues.set(meshPointI, sharedPts[pd.sharedPointAddr()[i]]);
1867         }
1868     }
1869 }
1870
1871
1872 template <unsigned nBits, class CombineOp>
1873 void Foam::syncTools::syncEdgeList
1874 (
1875     const polyMesh& mesh,
1876     PackedList<nBits>& edgeValues,
1877     const CombineOp& cop,
1878     const unsigned int nullValue
1879 )
1880 {
1881     if (edgeValues.size() != mesh.nEdges())
1882     {
1883         FatalErrorIn
1884         (
1885             "syncTools<unsigned nBits, class CombineOp>::syncEdgeList"
1886             "(const polyMesh&, PackedList<nBits>&, const CombineOp&"
1887             ", const unsigned int&)"
1888         )   << "Number of values " << edgeValues.size()
1889             << " is not equal to the number of edges in the mesh "
1890             << mesh.nEdges() << abort(FatalError);
1891     }
1892
1893     const polyBoundaryMesh& patches = mesh.boundaryMesh();
1894
1895     if (!hasCouples(patches))
1896     {
1897         return;
1898     }
1899
1900     // Patch data (proc patches only).
1901     List<List<unsigned int> > patchValues(patches.size());
1902
1903     if (Pstream::parRun())
1904     {
1905         // Send
1906
1907         forAll(patches, patchI)
1908         {
1909             if
1910             (
1911                 isA<processorPolyPatch>(patches[patchI])
1912              && patches[patchI].nEdges() > 0
1913             )
1914             {
1915                 const processorPolyPatch& procPatch =
1916                     refCast<const processorPolyPatch>(patches[patchI]);
1917
1918                 patchValues[patchI].setSize(procPatch.nEdges(), nullValue);
1919
1920                 const labelList& meshEdges = procPatch.meshEdges();
1921                 const labelList& neighbEdges = procPatch.neighbEdges();
1922
1923                 forAll(neighbEdges, edgeI)
1924                 {
1925                     label nbrEdgeI = neighbEdges[edgeI];
1926                     if (nbrEdgeI >= 0 && nbrEdgeI < procPatch.nEdges())
1927                     {
1928                         patchValues[patchI][nbrEdgeI] =
1929                             edgeValues.get(meshEdges[edgeI]);
1930                     }
1931                 }
1932
1933                 OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
1934                 toNbr << patchValues[patchI];
1935             }
1936         }
1937
1938
1939         // Receive and combine.
1940
1941         forAll(patches, patchI)
1942         {
1943             if
1944             (
1945                 isA<processorPolyPatch>(patches[patchI])
1946              && patches[patchI].nEdges() > 0
1947             )
1948             {
1949                 const processorPolyPatch& procPatch =
1950                     refCast<const processorPolyPatch>(patches[patchI]);
1951
1952                 {
1953                     IPstream fromNeighb
1954                     (
1955                         Pstream::blocking,
1956                         procPatch.neighbProcNo()
1957                     );
1958                     fromNeighb >> patchValues[patchI];
1959                 }
1960
1961                 patchValues[patchI].setSize(procPatch.nEdges(), nullValue);
1962
1963                 const labelList& meshEdges = procPatch.meshEdges();
1964
1965                 forAll(meshEdges, edgeI)
1966                 {
1967                     unsigned int patchVal = patchValues[patchI][edgeI];
1968                     label meshEdgeI = meshEdges[edgeI];
1969                     unsigned int edgeVal = edgeValues.get(meshEdgeI);
1970                     cop(edgeVal, patchVal);
1971                     edgeValues.set(meshEdgeI, edgeVal);
1972                 }
1973             }
1974         }
1975     }
1976
1977     // Do the cyclics.
1978     forAll(patches, patchI)
1979     {
1980         if (isA<cyclicPolyPatch>(patches[patchI]))
1981         {
1982             const cyclicPolyPatch& cycPatch =
1983                 refCast<const cyclicPolyPatch>(patches[patchI]);
1984
1985             const edgeList& coupledEdges = cycPatch.coupledEdges();
1986             const labelList& meshEdges = cycPatch.meshEdges();
1987
1988             forAll(coupledEdges, i)
1989             {
1990                 const edge& e = coupledEdges[i];
1991
1992                 label edge0 = meshEdges[e[0]];
1993                 label edge1 = meshEdges[e[1]];
1994
1995                 unsigned int val0 = edgeValues.get(edge0);
1996                 unsigned int t = val0;
1997                 unsigned int val1 = edgeValues.get(edge1);
1998
1999                 cop(t, val1);
2000                 edgeValues.set(edge0, t);
2001                 cop(val1, val0);
2002                 edgeValues.set(edge1, val1);
2003             }
2004         }
2005     }
2006
2007     // Synchronize multiple shared edges.
2008     const globalMeshData& pd = mesh.globalData();
2009
2010     if (pd.nGlobalEdges() > 0)
2011     {
2012         // Values on shared edges. Use unpacked storage for ease!
2013         List<unsigned int> sharedPts(pd.nGlobalEdges(), nullValue);
2014
2015         forAll(pd.sharedEdgeLabels(), i)
2016         {
2017             label meshEdgeI = pd.sharedEdgeLabels()[i];
2018             // Fill my entries in the shared edges
2019             sharedPts[pd.sharedEdgeAddr()[i]] = edgeValues.get(meshEdgeI);
2020         }
2021
2022         // Combine on master.
2023         Pstream::listCombineGather(sharedPts, cop);
2024         Pstream::listCombineScatter(sharedPts);
2025
2026         // Now we will all have the same information. Merge it back with
2027         // my local information.
2028         forAll(pd.sharedEdgeLabels(), i)
2029         {
2030             label meshEdgeI = pd.sharedEdgeLabels()[i];
2031             edgeValues.set(meshEdgeI, sharedPts[pd.sharedEdgeAddr()[i]]);
2032         }
2033     }
2034 }
2035
2036
2037 // ************************************************************************* //