src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.C

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