src/foam/meshes/polyMesh/globalMeshData/globalPoints.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 "globalPoints.H"
  27 #include "processorPolyPatch.H"
  28 #include "cyclicPolyPatch.H"
  29 #include "polyMesh.H"
  30
  31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  32
  33 defineTypeNameAndDebug(Foam::globalPoints, 0);
  34
  35
  36 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  37
  38 // Total number of points on processor patches. Is upper limit for number
  39 // of shared points
  40 Foam::label Foam::globalPoints::countPatchPoints
  41 (
  42     const polyBoundaryMesh& patches
  43 )
  44 {
  45     label nTotPoints = 0;
  46
  47     forAll(patches, patchI)
  48     {
  49         const polyPatch& pp = patches[patchI];
  50
  51         if
  52         (
  53             (Pstream::parRun() && isA<processorPolyPatch>(pp))
  54          || isA<cyclicPolyPatch>(pp)
  55         )
  56         {
  57             nTotPoints += pp.nPoints();
  58         }
  59     }
  60
  61     return nTotPoints;
  62 }
  63
  64
  65 // Collect all topological information about a point on a patch.
  66 // (this information is the patch faces using the point and the relative
  67 // position of the point in the face)
  68 void Foam::globalPoints::addToSend
  69 (
  70     const primitivePatch& pp,
  71     const label patchPointI,
  72     const procPointList& knownInfo,
  73
  74     DynamicList<label>& patchFaces,
  75     DynamicList<label>& indexInFace,
  76     DynamicList<procPointList>& allInfo
  77 )
  78 {
  79     label meshPointI = pp.meshPoints()[patchPointI];
  80
  81     // Add all faces using the point so we are sure we find it on the
  82     // other side.
  83     const labelList& pFaces = pp.pointFaces()[patchPointI];
  84
  85     forAll(pFaces, i)
  86     {
  87         label patchFaceI = pFaces[i];
  88
  89         const face& f = pp[patchFaceI];
  90
  91         patchFaces.append(patchFaceI);
  92         indexInFace.append(findIndex(f, meshPointI));
  93         allInfo.append(knownInfo);
  94     }
  95 }
  96
  97
  98 // Add nbrInfo to myInfo. Return true if anything changed.
  99 // nbrInfo is for a point a list of all the processors using it (and the
 100 // meshPoint label on that processor)
 101 bool Foam::globalPoints::mergeInfo
 102 (
 103     const procPointList& nbrInfo,
 104     procPointList& myInfo
 105 )
 106 {
 107     // Indices of entries in nbrInfo not yet in myInfo.
 108     DynamicList<label> newInfo(nbrInfo.size());
 109
 110     forAll(nbrInfo, i)
 111     {
 112         const procPoint& info = nbrInfo[i];
 113
 114         // Check if info already in myInfo.
 115         label index = -1;
 116
 117         forAll(myInfo, j)
 118         {
 119             if (myInfo[j] == info)
 120             {
 121                 // Already have information for processor/point combination
 122                 // in my list so skip.
 123                 index = j;
 124
 125                 break;
 126             }
 127         }
 128
 129         if (index == -1)
 130         {
 131             // Mark this information as being not yet in myInfo
 132             newInfo.append(i);
 133         }
 134     }
 135
 136     newInfo.shrink();
 137
 138     // Append all nbrInfos referenced in newInfo to myInfo.
 139
 140     label index = myInfo.size();
 141
 142     myInfo.setSize(index + newInfo.size());
 143
 144     forAll(newInfo, i)
 145     {
 146         myInfo[index++] = nbrInfo[newInfo[i]];
 147     }
 148
 149     // Did anything change?
 150     return newInfo.size() > 0;
 151 }
 152
 153
 154 // Updates database of current information on meshpoints with nbrInfo.
 155 // Uses mergeInfo above. Returns true if data kept for meshPointI changed.
 156 bool Foam::globalPoints::storeInfo
 157 (
 158     const procPointList& nbrInfo,
 159     const label meshPointI
 160 )
 161 {
 162     label infoChanged = false;
 163
 164     // Get the index into the procPoints list.
 165     Map<label>::iterator iter = meshToProcPoint_.find(meshPointI);
 166
 167     if (iter != meshToProcPoint_.end())
 168     {
 169         procPointList& knownInfo = procPoints_[iter()];
 170
 171         if (mergeInfo(nbrInfo, knownInfo))
 172         {
 173             infoChanged = true;
 174         }
 175     }
 176     else
 177     {
 178         procPointList knownInfo(1);
 179         knownInfo[0][0] = Pstream::myProcNo();
 180         knownInfo[0][1] = meshPointI;
 181
 182         if (mergeInfo(nbrInfo, knownInfo))
 183         {
 184             // Update addressing from into procPoints
 185             meshToProcPoint_.insert(meshPointI, procPoints_.size());
 186             // Insert into list of equivalences.
 187             procPoints_.append(knownInfo);
 188
 189             infoChanged = true;
 190         }
 191     }
 192     return infoChanged;
 193 }
 194
 195
 196 // Insert my own points into structure and mark as changed.
 197 void Foam::globalPoints::initOwnPoints
 198 (
 199     const bool allPoints,
 200     labelHashSet& changedPoints
 201 )
 202 {
 203     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 204
 205     forAll(patches, patchI)
 206     {
 207         const polyPatch& pp = patches[patchI];
 208
 209         if
 210         (
 211             (Pstream::parRun() && isA<processorPolyPatch>(pp))
 212          || isA<cyclicPolyPatch>(pp)
 213         )
 214         {
 215             const labelList& meshPoints = pp.meshPoints();
 216
 217             if (allPoints)
 218             {
 219                 // All points on patch
 220                 forAll(meshPoints, i)
 221                 {
 222                     label meshPointI = meshPoints[i];
 223
 224                     procPointList knownInfo(1);
 225                     knownInfo[0][0] = Pstream::myProcNo();
 226                     knownInfo[0][1] = meshPointI;
 227
 228                     // Update addressing from meshpoint to index in procPoints
 229                     meshToProcPoint_.insert(meshPointI, procPoints_.size());
 230                     // Insert into list of equivalences.
 231                     procPoints_.append(knownInfo);
 232
 233                     // Update changedpoints info.
 234                     changedPoints.insert(meshPointI);
 235                 }
 236             }
 237             else
 238             {
 239                 // Boundary points only
 240                 const labelList& boundaryPoints = pp.boundaryPoints();
 241
 242                 forAll(boundaryPoints, i)
 243                 {
 244                     label meshPointI = meshPoints[boundaryPoints[i]];
 245
 246                     procPointList knownInfo(1);
 247                     knownInfo[0][0] = Pstream::myProcNo();
 248                     knownInfo[0][1] = meshPointI;
 249
 250                     // Update addressing from meshpoint to index in procPoints
 251                     meshToProcPoint_.insert(meshPointI, procPoints_.size());
 252                     // Insert into list of equivalences.
 253                     procPoints_.append(knownInfo);
 254
 255                     // Update changedpoints info.
 256                     changedPoints.insert(meshPointI);
 257                 }
 258             }
 259         }
 260     }
 261 }
 262
 263
 264 // Send all my info on changedPoints_ to my neighbours.
 265 void Foam::globalPoints::sendPatchPoints(const labelHashSet& changedPoints)
 266  const
 267 {
 268     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 269
 270     forAll(patches, patchI)
 271     {
 272         const polyPatch& pp = patches[patchI];
 273
 274         if (Pstream::parRun() && isA<processorPolyPatch>(pp))
 275         {
 276             // Information to send:
 277             // patch face
 278             DynamicList<label> patchFaces(pp.nPoints());
 279             // index in patch face
 280             DynamicList<label> indexInFace(pp.nPoints());
 281             // all information I currently hold about this patchPoint
 282             DynamicList<procPointList> allInfo(pp.nPoints());
 283
 284
 285             // Now collect information on all mesh points mentioned in
 286             // changedPoints. Note that these points only should occur on
 287             // processorPatches (or rather this is a limitation!).
 288
 289             const labelList& meshPoints = pp.meshPoints();
 290
 291             forAll(meshPoints, patchPointI)
 292             {
 293                 label meshPointI = meshPoints[patchPointI];
 294
 295                 if (changedPoints.found(meshPointI))
 296                 {
 297                     label index = meshToProcPoint_[meshPointI];
 298
 299                     const procPointList& knownInfo = procPoints_[index];
 300
 301                     // Add my information about meshPointI to the send buffers
 302                     addToSend
 303                     (
 304                         pp,
 305                         patchPointI,
 306                         knownInfo,
 307
 308                         patchFaces,
 309                         indexInFace,
 310                         allInfo
 311                     );
 312                 }
 313             }
 314             patchFaces.shrink();
 315             indexInFace.shrink();
 316             allInfo.shrink();
 317
 318             // Send to neighbour
 319             {
 320                 const processorPolyPatch& procPatch =
 321                     refCast<const processorPolyPatch>(pp);
 322
 323                 if (debug)
 324                 {
 325                     Pout<< " Sending to "
 326                         << procPatch.neighbProcNo() << "   point information:"
 327                         << patchFaces.size() << endl;
 328                 }
 329
 330                 OPstream toNeighbour
 331                 (
 332                     Pstream::blocking,
 333                     procPatch.neighbProcNo()
 334                 );
 335
 336                 toNeighbour << patchFaces << indexInFace << allInfo;
 337             }
 338         }
 339     }
 340 }
 341
 342
 343 // Receive all my neighbours' information and merge with mine.
 344 // After finishing will have updated
 345 // - procPoints_ : all neighbour information merged in.
 346 // - meshToProcPoint_
 347 // - changedPoints: all meshPoints for which something changed.
 348 void Foam::globalPoints::receivePatchPoints(labelHashSet& changedPoints)
 349 {
 350     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 351
 352     // Reset changed points
 353     changedPoints.clear();
 354
 355     forAll(patches, patchI)
 356     {
 357         const polyPatch& pp = patches[patchI];
 358
 359         if (Pstream::parRun() && isA<processorPolyPatch>(pp))
 360         {
 361             const processorPolyPatch& procPatch =
 362                 refCast<const processorPolyPatch>(pp);
 363
 364             labelList patchFaces;
 365             labelList indexInFace;
 366             List<procPointList> nbrInfo;
 367
 368             {
 369                 IPstream fromNeighbour
 370                 (
 371                     Pstream::blocking,
 372                     procPatch.neighbProcNo()
 373                 );
 374                 fromNeighbour >> patchFaces >> indexInFace >> nbrInfo;
 375             }
 376
 377             if (debug)
 378             {
 379                 Pout<< " Received from "
 380                     << procPatch.neighbProcNo() << "   point information:"
 381                     << patchFaces.size() << endl;
 382             }
 383
 384             forAll(patchFaces, i)
 385             {
 386                 const face& f = pp[patchFaces[i]];
 387
 388                 // Get index in this face from index on face on other side.
 389                 label index = (f.size() - indexInFace[i]) % f.size();
 390
 391                 // Get the meshpoint on my side
 392                 label meshPointI = f[index];
 393
 394                 //const procPointList& info = nbrInfo[i];
 395                 //Pout << "Received for my coord "
 396                 //    << mesh_.points()[meshPointI];
 397                 //
 398                 //forAll(info, j)
 399                 //{
 400                 //    Pout<< ' ' <<info[j];
 401                 //}
 402                 //Pout<< endl;
 403
 404                 if (storeInfo(nbrInfo[i], meshPointI))
 405                 {
 406                     changedPoints.insert(meshPointI);
 407                 }
 408             }
 409         }
 410         else if (isA<cyclicPolyPatch>(pp))
 411         {
 412             // Handle cyclics: send lower half to upper half and vice versa.
 413             // Or since they both are in memory just do it point by point.
 414
 415             const cyclicPolyPatch& cycPatch =
 416                 refCast<const cyclicPolyPatch>(pp);
 417
 418             const labelList& meshPoints = pp.meshPoints();
 419
 420             //const edgeList& connections = cycPatch.coupledPoints();
 421             const edgeList connections(coupledPoints(cycPatch));
 422
 423             forAll(connections, i)
 424             {
 425                 const edge& e = connections[i];
 426
 427                 label meshPointA = meshPoints[e[0]];
 428                 label meshPointB = meshPoints[e[1]];
 429
 430                 // Do we have information on meshPointA?
 431                 Map<label>::iterator procPointA =
 432                     meshToProcPoint_.find(meshPointA);
 433
 434                 if (procPointA != meshToProcPoint_.end())
 435                 {
 436                     // Store A info onto pointB
 437                     if (storeInfo(procPoints_[procPointA()], meshPointB))
 438                     {
 439                         changedPoints.insert(meshPointB);
 440                     }
 441                 }
 442
 443                 // Same for info on pointB
 444                 Map<label>::iterator procPointB =
 445                     meshToProcPoint_.find(meshPointB);
 446
 447                 if (procPointB != meshToProcPoint_.end())
 448                 {
 449                     // Store B info onto pointA
 450                     if (storeInfo(procPoints_[procPointB()], meshPointA))
 451                     {
 452                         changedPoints.insert(meshPointA);
 453                     }
 454                 }
 455             }
 456         }
 457     }
 458 }
 459
 460
 461 // Remove entries which are handled by normal face-face communication. I.e.
 462 // those points where the equivalence list is only me and my (face)neighbour
 463 void Foam::globalPoints::remove(const Map<label>& directNeighbours)
 464 {
 465     // Save old ones.
 466     Map<label> oldMeshToProcPoint(meshToProcPoint_);
 467     meshToProcPoint_.clear();
 468
 469     List<procPointList> oldProcPoints;
 470     oldProcPoints.transfer(procPoints_);
 471
 472     // Go through all equivalences
 473     for
 474     (
 475         Map<label>::const_iterator iter = oldMeshToProcPoint.begin();
 476         iter != oldMeshToProcPoint.end();
 477         ++iter
 478     )
 479     {
 480         label meshPointI = iter.key();
 481         const procPointList& pointInfo = oldProcPoints[iter()];
 482
 483         if (pointInfo.size() == 2)
 484         {
 485             // I will be in this equivalence list.
 486             // Check whether my direct (=face) neighbour
 487             // is in it. This would be an ordinary connection and can be
 488             // handled by normal face-face connectivity.
 489
 490             const procPoint& a = pointInfo[0];
 491             const procPoint& b = pointInfo[1];
 492
 493             if
 494             (
 495                 (a[0] == Pstream::myProcNo() && directNeighbours.found(a[1]))
 496              || (b[0] == Pstream::myProcNo() && directNeighbours.found(b[1]))
 497             )
 498             {
 499                 // Normal faceNeighbours
 500                 if (a[0] == Pstream::myProcNo())
 501                 {
 502                     //Pout<< "Removing direct neighbour:"
 503                     //    << mesh_.points()[a[1]]
 504                     //    << endl;
 505                 }
 506                 else if (b[0] == Pstream::myProcNo())
 507                 {
 508                     //Pout<< "Removing direct neighbour:"
 509                     //    << mesh_.points()[b[1]]
 510                     //    << endl;
 511                 }
 512             }
 513             else
 514             {
 515                 // This condition will be very rare: points are used by
 516                 // two processors which are not face-face connected.
 517                 // e.g.
 518                 // +------+------+
 519                 // | wall |  B   |
 520                 // +------+------+
 521                 // |   A  | wall |
 522                 // +------+------+
 523                 // Processor A and B share a point. Note that this only will
 524                 // be found if the two domains are face connected at all
 525                 // (not shown in the picture)
 526
 527                 meshToProcPoint_.insert(meshPointI, procPoints_.size());
 528                 procPoints_.append(pointInfo);
 529             }
 530         }
 531         else if (pointInfo.size() == 1)
 532         {
 533             // This happens for 'wedge' like cyclics where the two halves
 534             // come together in the same point so share the same meshPoint.
 535             // So this meshPoint will have info of size one only.
 536             if
 537             (
 538                 pointInfo[0][0] != Pstream::myProcNo()
 539              || !directNeighbours.found(pointInfo[0][1])
 540             )
 541             {
 542                 meshToProcPoint_.insert(meshPointI, procPoints_.size());
 543                 procPoints_.append(pointInfo);
 544             }
 545         }
 546         else
 547         {
 548             meshToProcPoint_.insert(meshPointI, procPoints_.size());
 549             procPoints_.append(pointInfo);
 550         }
 551     }
 552
 553     procPoints_.shrink();
 554 }
 555
 556
 557 // Get (indices of) points for which I am master (= lowest numbered point on
 558 // lowest numbered processor).
 559 // (equivalence lists should be complete by now)
 560 Foam::labelList Foam::globalPoints::getMasterPoints() const
 561 {
 562     labelList masterPoints(nPatchPoints_);
 563     label nMaster = 0;
 564
 565     // Go through all equivalences and determine meshPoints where I am master.
 566     for
 567     (
 568         Map<label>::const_iterator iter = meshToProcPoint_.begin();
 569         iter != meshToProcPoint_.end();
 570         ++iter
 571     )
 572     {
 573         label meshPointI = iter.key();
 574         const procPointList& pointInfo = procPoints_[iter()];
 575
 576         if (pointInfo.size() < 2)
 577         {
 578             // Points should have an equivalence list >= 2 since otherwise
 579             // they would be direct neighbours and have been removed in the
 580             // call to 'remove'.
 581             Pout<< "MeshPoint:" << meshPointI
 582                 << " coord:" << mesh_.points()[meshPointI]
 583                 << " has no corresponding point on a neighbouring processor"
 584                 << endl;
 585             FatalErrorIn("globalPoints::getMasterPoints()")
 586                 << '[' << Pstream::myProcNo() << ']'
 587                 << " MeshPoint:" << meshPointI
 588                 << " coord:" << mesh_.points()[meshPointI]
 589                 << " has no corresponding point on a neighbouring processor"
 590                 << abort(FatalError);
 591         }
 592         else
 593         {
 594             // Find lowest processor and lowest mesh point on this processor.
 595             label lowestProcI = labelMax;
 596             label lowestPointI = labelMax;
 597
 598             forAll(pointInfo, i)
 599             {
 600                 label proc = pointInfo[i][0];
 601
 602                 if (proc < lowestProcI)
 603                 {
 604                     lowestProcI = proc;
 605                     lowestPointI = pointInfo[i][1];
 606                 }
 607                 else if (proc == lowestProcI)
 608                 {
 609                     lowestPointI = min(lowestPointI, pointInfo[i][1]);
 610                 }
 611             }
 612
 613             if
 614             (
 615                 lowestProcI == Pstream::myProcNo()
 616              && lowestPointI == meshPointI
 617             )
 618             {
 619                 // I am lowest numbered processor and point. Add to my list.
 620                 masterPoints[nMaster++] = meshPointI;
 621             }
 622         }
 623     }
 624
 625     masterPoints.setSize(nMaster);
 626
 627     return masterPoints;
 628 }
 629
 630
 631 // Send subset of lists
 632 void Foam::globalPoints::sendSharedPoints(const labelList& changedIndices) const
 633 {
 634     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 635
 636     forAll(patches, patchI)
 637     {
 638         const polyPatch& pp = patches[patchI];
 639
 640         if (Pstream::parRun() && isA<processorPolyPatch>(pp))
 641         {
 642             const processorPolyPatch& procPatch =
 643                 refCast<const processorPolyPatch>(pp);
 644
 645             OPstream toNeighbour(Pstream::blocking, procPatch.neighbProcNo());
 646
 647             if (debug)
 648             {
 649                 Pout<< "Sending to " << procPatch.neighbProcNo()
 650                     << "  changed sharedPoints info:"
 651                     << changedIndices.size() << endl;
 652             }
 653
 654             toNeighbour
 655                 << UIndirectList<label>(sharedPointAddr_, changedIndices)()
 656                 << UIndirectList<label>(sharedPointLabels_, changedIndices)();
 657         }
 658     }
 659 }
 660
 661
 662 // Receive shared point indices for all my shared points. Note that since
 663 // there are only a few here we can build a reverse map using the meshpoint
 664 // instead of doing all this relative point indexing (patch face + index in
 665 // face) as in send/receivePatchPoints
 666 void Foam::globalPoints::receiveSharedPoints(labelList& changedIndices)
 667 {
 668     changedIndices.setSize(sharedPointAddr_.size());
 669     label nChanged = 0;
 670
 671     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 672
 673     // Receive and set shared points
 674     forAll(patches, patchI)
 675     {
 676         const polyPatch& pp = patches[patchI];
 677
 678         if (Pstream::parRun() && isA<processorPolyPatch>(pp))
 679         {
 680             const processorPolyPatch& procPatch =
 681                 refCast<const processorPolyPatch>(pp);
 682
 683             // Map from neighbouring meshPoint to sharedPoint)
 684             Map<label> nbrSharedPoints(sharedPointAddr_.size());
 685
 686             {
 687                 // Receive meshPoints on neighbour and sharedPoints and build
 688                 // map from it. Note that we could have built the map on the
 689                 // neighbour and sent it over.
 690                 labelList nbrSharedPointAddr;
 691                 labelList nbrSharedPointLabels;
 692
 693                 {
 694                     IPstream fromNeighbour
 695                     (
 696                         Pstream::blocking,
 697                         procPatch.neighbProcNo()
 698                     );
 699                     fromNeighbour >> nbrSharedPointAddr >> nbrSharedPointLabels;
 700                 }
 701
 702                 // Insert into to map
 703                 forAll(nbrSharedPointLabels, i)
 704                 {
 705                     nbrSharedPoints.insert
 706                     (
 707                         nbrSharedPointLabels[i], // meshpoint on neighbour
 708                         nbrSharedPointAddr[i]    // sharedPoint label
 709                     );
 710                 }
 711             }
 712
 713
 714             // Merge into whatever information I hold.
 715             for
 716             (
 717                 Map<label>::const_iterator iter = meshToProcPoint_.begin();
 718                 iter != meshToProcPoint_.end();
 719                 ++iter
 720             )
 721             {
 722                 label meshPointI = iter.key();
 723                 label index = iter();
 724
 725                 if (sharedPointAddr_[index] == -1)
 726                 {
 727                     // No shared point known yet for this meshPoint.
 728                     // See if was received from neighbour.
 729                     const procPointList& knownInfo = procPoints_[index];
 730
 731                     // Check through the whole equivalence list for any
 732                     // point from the neighbour.
 733                     forAll(knownInfo, j)
 734                     {
 735                         const procPoint& info = knownInfo[j];
 736
 737                         if
 738                         (
 739                             (info[0] == procPatch.neighbProcNo())
 740                          && nbrSharedPoints.found(info[1])
 741                         )
 742                         {
 743                             // So this knownInfo contains the neighbour point
 744                             label sharedPointI = nbrSharedPoints[info[1]];
 745
 746                             sharedPointAddr_[index] = sharedPointI;
 747                             sharedPointLabels_[index] = meshPointI;
 748                             changedIndices[nChanged++] = index;
 749
 750                             break;
 751                         }
 752                     }
 753                 }
 754             }
 755         }
 756         else if (isA<cyclicPolyPatch>(pp))
 757         {
 758             const cyclicPolyPatch& cycPatch =
 759                 refCast<const cyclicPolyPatch>(pp);
 760
 761             // Build map from meshPoint to sharedPoint
 762             Map<label> meshToSharedPoint(sharedPointAddr_.size());
 763             forAll(sharedPointLabels_, i)
 764             {
 765                 label meshPointI = sharedPointLabels_[i];
 766
 767                 meshToSharedPoint.insert(meshPointI, sharedPointAddr_[i]);
 768             }
 769
 770             // Sync all info.
 771             //const edgeList& connections = cycPatch.coupledPoints();
 772             const edgeList connections(coupledPoints(cycPatch));
 773
 774             forAll(connections, i)
 775             {
 776                 const edge& e = connections[i];
 777
 778                 label meshPointA = pp.meshPoints()[e[0]];
 779                 label meshPointB = pp.meshPoints()[e[1]];
 780
 781                 // Do we already have shared point for meshPointA?
 782                 Map<label>::iterator fndA = meshToSharedPoint.find(meshPointA);
 783                 Map<label>::iterator fndB = meshToSharedPoint.find(meshPointB);
 784
 785                 if (fndA != meshToSharedPoint.end())
 786                 {
 787                     if (fndB != meshToSharedPoint.end())
 788                     {
 789                         if (fndA() != fndB())
 790                         {
 791                             FatalErrorIn
 792                             (
 793                                 "globalPoints::receiveSharedPoints"
 794                                 "(labelList&)"
 795                             )   << "On patch " << pp.name()
 796                                 << " connected points " << meshPointA
 797                                 << ' ' << mesh_.points()[meshPointA]
 798                                 << " and " << meshPointB
 799                                 << ' ' << mesh_.points()[meshPointB]
 800                                 << " are mapped to different shared points: "
 801                                 << fndA() << " and " << fndB()
 802                                 << abort(FatalError);
 803                         }
 804                     }
 805                     else
 806                     {
 807                         // No shared point yet for B.
 808                         label sharedPointI = fndA();
 809
 810                         // Store shared point for meshPointB
 811                         label index = meshToProcPoint_[meshPointB];
 812
 813                         sharedPointAddr_[index] = sharedPointI;
 814                         sharedPointLabels_[index] = meshPointB;
 815                         changedIndices[nChanged++] = index;
 816                     }
 817                 }
 818                 else
 819                 {
 820                     // No shared point yet for A.
 821                     if (fndB != meshToSharedPoint.end())
 822                     {
 823                         label sharedPointI = fndB();
 824
 825                         // Store shared point for meshPointA
 826                         label index = meshToProcPoint_[meshPointA];
 827
 828                         sharedPointAddr_[index] = sharedPointI;
 829                         sharedPointLabels_[index] = meshPointA;
 830                         changedIndices[nChanged++] = index;
 831                     }
 832                 }
 833             }
 834         }
 835     }
 836
 837     changedIndices.setSize(nChanged);
 838 }
 839
 840
 841 Foam::edgeList Foam::globalPoints::coupledPoints(const cyclicPolyPatch& pp)
 842 {
 843     // Look at cyclic patch as two halves, A and B.
 844     // Now all we know is that relative face index in halfA is same
 845     // as coupled face in halfB and also that the 0th vertex
 846     // corresponds.
 847
 848     // From halfA point to halfB or -1.
 849     labelList coupledPoint(pp.nPoints(), -1);
 850
 851     for (label patchFaceA = 0; patchFaceA < pp.size()/2; patchFaceA++)
 852     {
 853         const face& fA = pp.localFaces()[patchFaceA];
 854
 855         forAll(fA, indexA)
 856         {
 857             label patchPointA = fA[indexA];
 858
 859             if (coupledPoint[patchPointA] == -1)
 860             {
 861                 const face& fB = pp.localFaces()[patchFaceA + pp.size()/2];
 862
 863                 label indexB = (fB.size() - indexA) % fB.size();
 864
 865                 coupledPoint[patchPointA] = fB[indexB];
 866             }
 867         }
 868     }
 869
 870     edgeList connected(pp.nPoints());
 871
 872     // Extract coupled points.
 873     label connectedI = 0;
 874
 875     forAll(coupledPoint, i)
 876     {
 877         if (coupledPoint[i] != -1)
 878         {
 879             connected[connectedI++] = edge(i, coupledPoint[i]);
 880         }
 881     }
 882
 883     connected.setSize(connectedI);
 884
 885     return connected;
 886 }
 887
 888
 889 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 890
 891 // Construct from components
 892 Foam::globalPoints::globalPoints(const polyMesh& mesh)
 893 :
 894     mesh_(mesh),
 895     nPatchPoints_(countPatchPoints(mesh.boundaryMesh())),
 896     procPoints_(nPatchPoints_),
 897     meshToProcPoint_(nPatchPoints_),
 898     sharedPointAddr_(0),
 899     sharedPointLabels_(0),
 900     nGlobalPoints_(0)
 901 {
 902     if (debug)
 903     {
 904         Pout<< "globalPoints::globalPoints(const polyMesh&) : "
 905             << "doing processor to processor communication to get sharedPoints"
 906             << endl;
 907     }
 908
 909     labelHashSet changedPoints(nPatchPoints_);
 910
 911     // Initialize procPoints with my patch points. Keep track of points
 912     // inserted (in changedPoints)
 913     // There are two possible forms of this:
 914     // - initialize with all patch points (allPoints = true). This causes all
 915     //   patch points to be exchanged so a lot of information gets stored and
 916     //   transferred. This all gets filtered out later when removing the
 917     //   equivalence lists of size 2.
 918     // - initialize with boundary points of patches only (allPoints = false).
 919     //   This should work for all decompositions except extreme ones where a
 920     //   shared point is not on the boundary of any processor patches using it.
 921     //   This would happen if a domain was pinched such that two patches share
 922     //   a point or edge.
 923     initOwnPoints(true, changedPoints);
 924
 925     // Do one exchange iteration to get neighbour points.
 926     sendPatchPoints(changedPoints);
 927     receivePatchPoints(changedPoints);
 928
 929
 930     // Save neighbours reachable through face-face communication.
 931     Map<label> neighbourList(meshToProcPoint_);
 932
 933
 934     // Exchange until nothing changes on all processors.
 935     bool changed = false;
 936
 937     do
 938     {
 939         sendPatchPoints(changedPoints);
 940         receivePatchPoints(changedPoints);
 941
 942         changed = changedPoints.size() > 0;
 943         reduce(changed, orOp<bool>());
 944
 945     } while (changed);
 946
 947
 948     // Remove direct neighbours from point equivalences.
 949     remove(neighbourList);
 950
 951
 952     //Pout<< "After removing locally connected points:" << endl;
 953     //for
 954     //(
 955     //    Map<label>::const_iterator iter = meshToProcPoint_.begin();
 956     //    iter != meshToProcPoint_.end();
 957     //    ++iter
 958     //)
 959     //{
 960     //    label meshPointI = iter.key();
 961     //    const procPointList& pointInfo = procPoints_[iter()];
 962     //
 963     //    forAll(pointInfo, i)
 964     //    {
 965     //        Pout<< "    pointI:" << meshPointI << ' '
 966     //            << mesh.points()[meshPointI]
 967     //            << " connected to proc " << pointInfo[i][0]
 968     //            << " point:" << pointInfo[i][1]
 969     //        << endl;
 970     //    }
 971     //}
 972
 973
 974     // We now have - in procPoints_ - a list of points which are shared between
 975     // multiple processors. These are the ones for which are sharedPoint
 976     // needs to be determined. This is done by having the lowest numbered
 977     // processor in the equivalence list 'ask' for a sharedPoint number
 978     // and then distribute it across processor patches to the non-master
 979     // processors. Note: below piece of coding is not very efficient. Uses
 980     // a Map where possibly it shouldn't
 981
 982     // Initialize sharedPoint addressing. Is for every entry in procPoints_
 983     // the sharedPoint.
 984     sharedPointAddr_.setSize(meshToProcPoint_.size());
 985     sharedPointAddr_ = -1;
 986     sharedPointLabels_.setSize(meshToProcPoint_.size());
 987     sharedPointLabels_ = -1;
 988
 989
 990     // Get points for which I am master (lowest numbered proc)
 991     labelList masterPoints(getMasterPoints());
 992
 993
 994     // Determine number of master points on all processors.
 995     labelList sharedPointSizes(Pstream::nProcs());
 996     sharedPointSizes[Pstream::myProcNo()] = masterPoints.size();
 997
 998     Pstream::gatherList(sharedPointSizes);
 999     Pstream::scatterList(sharedPointSizes);
1000
1001     if (debug)
1002     {
1003         Pout<< "sharedPointSizes:" << sharedPointSizes << endl;
1004     }
1005
1006     // Get total number of shared points
1007     nGlobalPoints_ = 0;
1008     forAll(sharedPointSizes, procI)
1009     {
1010         nGlobalPoints_ += sharedPointSizes[procI];
1011     }
1012
1013     // Assign sharedPoint labels. Every processor gets assigned consecutive
1014     // numbers for its master points.
1015     // These are assigned in processor order so processor0 gets
1016     // 0..masterPoints.size()-1 etc.
1017
1018     // My point labels start after those of lower numbered processors
1019     label sharedPointI = 0;
1020     for (label procI = 0; procI < Pstream::myProcNo(); procI++)
1021     {
1022         sharedPointI += sharedPointSizes[procI];
1023     }
1024
1025     forAll(masterPoints, i)
1026     {
1027         label meshPointI = masterPoints[i];
1028
1029         label index = meshToProcPoint_[meshPointI];
1030
1031         sharedPointLabels_[index] = meshPointI;
1032         sharedPointAddr_[index] = sharedPointI++;
1033     }
1034
1035
1036     // Now we have a sharedPointLabel for some of the entries in procPoints.
1037     // Send this information to neighbours. Receive their information.
1038     // Loop until nothing changes.
1039
1040     // Initial subset to send is points for which I have sharedPoints
1041     labelList changedIndices(sharedPointAddr_.size());
1042     label nChanged = 0;
1043
1044     forAll(sharedPointAddr_, i)
1045     {
1046         if (sharedPointAddr_[i] != -1)
1047         {
1048             changedIndices[nChanged++] = i;
1049         }
1050     }
1051     changedIndices.setSize(nChanged);
1052
1053     changed = false;
1054
1055     do
1056     {
1057         if (debug)
1058         {
1059             Pout<< "Determined " << changedIndices.size() << " shared points."
1060                 << " Exchanging them" << endl;
1061         }
1062         sendSharedPoints(changedIndices);
1063         receiveSharedPoints(changedIndices);
1064
1065         changed = changedIndices.size() > 0;
1066         reduce(changed, orOp<bool>());
1067
1068     } while (changed);
1069
1070
1071     forAll(sharedPointLabels_, i)
1072     {
1073         if (sharedPointLabels_[i] == -1)
1074         {
1075             FatalErrorIn("globalPoints::globalPoints(const polyMesh& mesh)")
1076                 << "Problem: shared point on processor " << Pstream::myProcNo()
1077                 << " not set at index " << sharedPointLabels_[i] << endl
1078                 << "This might mean the individual processor domains are not"
1079                 << " connected and the overall domain consists of multiple"
1080                 << " regions. You can check this with checkMesh"
1081                 << abort(FatalError);
1082         }
1083     }
1084
1085     if (debug)
1086     {
1087         Pout<< "globalPoints::globalPoints(const polyMesh&) : "
1088             << "Finished global points" << endl;
1089     }
1090 }
1091
1092
1093 // ************************************************************************* //