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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "globalMeshData.H"
  27 #include "Time.H"
  28 #include "Pstream.H"
  29 #include "PstreamCombineReduceOps.H"
  30 #include "processorPolyPatch.H"
  31 #include "demandDrivenData.H"
  32 #include "globalPoints.H"
  33 #include "polyMesh.H"
  34 #include "mapDistribute.H"
  35 #include "labelIOList.H"
  36 #include "PackedList.H"
  37 #include "mergePoints.H"
  38 #include "matchPoints.H"
  39 #include "OFstream.H"
  40 #include "globalIndexAndTransform.H"
  41
  42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  43
  44 defineTypeNameAndDebug(Foam::globalMeshData, 0);
  45
  46 // Geometric matching tolerance. Factor of mesh bounding box.
  47 const Foam::scalar Foam::globalMeshData::matchTol_ = 1E-8;
  48
  49 namespace Foam
  50 {
  51 template<>
  52 class minEqOp<labelPair>
  53 {
  54 public:
  55     void operator()(labelPair& x, const labelPair& y) const
  56     {
  57         x[0] = min(x[0], y[0]);
  58         x[1] = min(x[1], y[1]);
  59     }
  60 };
  61 }
  62
  63
  64 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  65
  66 // Collect processor patch addressing.
  67 void Foam::globalMeshData::initProcAddr()
  68 {
  69     processorPatchIndices_.setSize(mesh_.boundaryMesh().size());
  70     processorPatchIndices_ = -1;
  71
  72     processorPatchNeighbours_.setSize(mesh_.boundaryMesh().size());
  73     processorPatchNeighbours_ = -1;
  74
  75     // Construct processor patch indexing. processorPatchNeighbours_ only
  76     // set if running in parallel!
  77     processorPatches_.setSize(mesh_.boundaryMesh().size());
  78
  79     label nNeighbours = 0;
  80
  81     forAll(mesh_.boundaryMesh(), patchi)
  82     {
  83         if (isA<processorPolyPatch>(mesh_.boundaryMesh()[patchi]))
  84         {
  85             processorPatches_[nNeighbours] = patchi;
  86             processorPatchIndices_[patchi] = nNeighbours++;
  87         }
  88     }
  89     processorPatches_.setSize(nNeighbours);
  90
  91
  92     if (Pstream::parRun())
  93     {
  94         PstreamBuffers pBufs(Pstream::nonBlocking);
  95
  96         // Send indices of my processor patches to my neighbours
  97         forAll(processorPatches_, i)
  98         {
  99             label patchi = processorPatches_[i];
 100
 101             UOPstream toNeighbour
 102             (
 103                 refCast<const processorPolyPatch>
 104                 (
 105                     mesh_.boundaryMesh()[patchi]
 106                 ).neighbProcNo(),
 107                 pBufs
 108             );
 109
 110             toNeighbour << processorPatchIndices_[patchi];
 111         }
 112
 113         pBufs.finishedSends();
 114
 115         forAll(processorPatches_, i)
 116         {
 117             label patchi = processorPatches_[i];
 118
 119             UIPstream fromNeighbour
 120             (
 121                 refCast<const processorPolyPatch>
 122                 (
 123                     mesh_.boundaryMesh()[patchi]
 124                 ).neighbProcNo(),
 125                 pBufs
 126             );
 127
 128             fromNeighbour >> processorPatchNeighbours_[patchi];
 129         }
 130     }
 131 }
 132
 133
 134 void Foam::globalMeshData::calcSharedPoints() const
 135 {
 136     if
 137     (
 138         nGlobalPoints_ != -1
 139      || sharedPointLabelsPtr_.valid()
 140      || sharedPointAddrPtr_.valid()
 141     )
 142     {
 143         FatalErrorIn("globalMeshData::calcSharedPoints()")
 144             << "Shared point addressing already done" << abort(FatalError);
 145     }
 146
 147     // Calculate all shared points (exclude points that are only
 148     // on two coupled patches). This does all the hard work.
 149     globalPoints parallelPoints(mesh_, false, true);
 150
 151     // Count the number of master points
 152     label nMaster = 0;
 153     forAll(parallelPoints.pointPoints(), i)
 154     {
 155         const labelList& pPoints = parallelPoints.pointPoints()[i];
 156         const labelList& transPPoints =
 157             parallelPoints.transformedPointPoints()[i];
 158
 159         if (pPoints.size()+transPPoints.size() > 0)
 160         {
 161             nMaster++;
 162         }
 163     }
 164
 165     // Allocate global numbers
 166     globalIndex masterNumbering(nMaster);
 167
 168     nGlobalPoints_ = masterNumbering.size();
 169
 170
 171     // Push master number to slaves
 172     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 173     // 1. Fill master and slave slots
 174     nMaster = 0;
 175     labelList master(parallelPoints.map().constructSize(), -1);
 176     forAll(parallelPoints.pointPoints(), i)
 177     {
 178         const labelList& pPoints = parallelPoints.pointPoints()[i];
 179         const labelList& transPPoints =
 180             parallelPoints.transformedPointPoints()[i];
 181
 182         if (pPoints.size()+transPPoints.size() > 0)
 183         {
 184             master[i] = masterNumbering.toGlobal(nMaster);
 185             forAll(pPoints, j)
 186             {
 187                 master[pPoints[j]] = master[i];
 188             }
 189             forAll(transPPoints, j)
 190             {
 191                 master[transPPoints[j]] = master[i];
 192             }
 193             nMaster++;
 194         }
 195     }
 196
 197
 198     // 2. Push slave slots back to local storage on originating processor
 199     // For all the four types of points:
 200     // - local master : already set
 201     // - local transformed slave point : the reverse transform at
 202     //   reverseDistribute will have copied it back to its originating local
 203     //   point
 204     // - remote untransformed slave point : sent back to originating processor
 205     // - remote transformed slave point : the reverse transform will
 206     //   copy it back into the remote slot which then gets sent back to
 207     //   originating processor
 208
 209     parallelPoints.map().reverseDistribute
 210     (
 211         parallelPoints.map().constructSize(),
 212         master
 213     );
 214
 215
 216     // Collect all points that are a master or refer to a master.
 217     nMaster = 0;
 218     forAll(parallelPoints.pointPoints(), i)
 219     {
 220         if (master[i] != -1)
 221         {
 222             nMaster++;
 223         }
 224     }
 225
 226     sharedPointLabelsPtr_.reset(new labelList(nMaster));
 227     labelList& sharedPointLabels = sharedPointLabelsPtr_();
 228     sharedPointAddrPtr_.reset(new labelList(nMaster));
 229     labelList& sharedPointAddr = sharedPointAddrPtr_();
 230     nMaster = 0;
 231
 232     forAll(parallelPoints.pointPoints(), i)
 233     {
 234         if (master[i] != -1)
 235         {
 236             // I am master or slave
 237             sharedPointLabels[nMaster] = i;
 238             sharedPointAddr[nMaster] = master[i];
 239             nMaster++;
 240         }
 241     }
 242
 243     if (debug)
 244     {
 245         Pout<< "globalMeshData : nGlobalPoints_:" << nGlobalPoints_ << nl
 246             << "globalMeshData : sharedPointLabels_:"
 247             << sharedPointLabelsPtr_().size() << nl
 248             << "globalMeshData : sharedPointAddr_:"
 249             << sharedPointAddrPtr_().size() << endl;
 250     }
 251 }
 252
 253
 254 // Given information about locally used edges allocate global shared edges.
 255 void Foam::globalMeshData::countSharedEdges
 256 (
 257     const EdgeMap<labelList>& procSharedEdges,
 258     EdgeMap<label>& globalShared,
 259     label& sharedEdgeI
 260 )
 261 {
 262     // Count occurrences of procSharedEdges in global shared edges table.
 263     forAllConstIter(EdgeMap<labelList>, procSharedEdges, iter)
 264     {
 265         const edge& e = iter.key();
 266
 267         EdgeMap<label>::iterator globalFnd = globalShared.find(e);
 268
 269         if (globalFnd == globalShared.end())
 270         {
 271             // First time occurrence of this edge. Check how many we are adding.
 272             if (iter().size() == 1)
 273             {
 274                 // Only one edge. Mark with special value.
 275                 globalShared.insert(e, -1);
 276             }
 277             else
 278             {
 279                 // Edge used more than once (even by local shared edges alone)
 280                 // so allocate proper shared edge label.
 281                 globalShared.insert(e, sharedEdgeI++);
 282             }
 283         }
 284         else
 285         {
 286             if (globalFnd() == -1)
 287             {
 288                 // Second time occurence of this edge. Assign proper
 289                 // edge label.
 290                 globalFnd() = sharedEdgeI++;
 291             }
 292         }
 293     }
 294 }
 295
 296
 297 // Shared edges are shared between multiple processors. By their nature both
 298 // of their endpoints are shared points. (but not all edges using two shared
 299 // points are shared edges! There might e.g. be an edge between two unrelated
 300 // clusters of shared points)
 301 void Foam::globalMeshData::calcSharedEdges() const
 302 {
 303     if
 304     (
 305         nGlobalEdges_ != -1
 306      || sharedEdgeLabelsPtr_.valid()
 307      || sharedEdgeAddrPtr_.valid()
 308     )
 309     {
 310         FatalErrorIn("globalMeshData::calcSharedEdges()")
 311             << "Shared edge addressing already done" << abort(FatalError);
 312     }
 313
 314
 315     const labelList& sharedPtAddr = sharedPointAddr();
 316     const labelList& sharedPtLabels = sharedPointLabels();
 317
 318     // Since don't want to construct pointEdges for whole mesh create
 319     // Map for all shared points.
 320     Map<label> meshToShared(2*sharedPtLabels.size());
 321     forAll(sharedPtLabels, i)
 322     {
 323         meshToShared.insert(sharedPtLabels[i], i);
 324     }
 325
 326     // Find edges using shared points. Store correspondence to local edge
 327     // numbering. Note that multiple local edges can have the same shared
 328     // points! (for cyclics or separated processor patches)
 329     EdgeMap<labelList> localShared(2*sharedPtAddr.size());
 330
 331     const edgeList& edges = mesh_.edges();
 332
 333     forAll(edges, edgeI)
 334     {
 335         const edge& e = edges[edgeI];
 336
 337         Map<label>::const_iterator e0Fnd = meshToShared.find(e[0]);
 338
 339         if (e0Fnd != meshToShared.end())
 340         {
 341             Map<label>::const_iterator e1Fnd = meshToShared.find(e[1]);
 342
 343             if (e1Fnd != meshToShared.end())
 344             {
 345                 // Found edge which uses shared points. Probably shared.
 346
 347                 // Construct the edge in shared points (or rather global indices
 348                 // of the shared points)
 349                 edge sharedEdge
 350                 (
 351                     sharedPtAddr[e0Fnd()],
 352                     sharedPtAddr[e1Fnd()]
 353                 );
 354
 355                 EdgeMap<labelList>::iterator iter =
 356                     localShared.find(sharedEdge);
 357
 358                 if (iter == localShared.end())
 359                 {
 360                     // First occurrence of this point combination. Store.
 361                     localShared.insert(sharedEdge, labelList(1, edgeI));
 362                 }
 363                 else
 364                 {
 365                     // Add this edge to list of edge labels.
 366                     labelList& edgeLabels = iter();
 367
 368                     label sz = edgeLabels.size();
 369                     edgeLabels.setSize(sz+1);
 370                     edgeLabels[sz] = edgeI;
 371                 }
 372             }
 373         }
 374     }
 375
 376
 377     // Now we have a table on every processors which gives its edges which use
 378     // shared points. Send this all to the master and have it allocate
 379     // global edge numbers for it. But only allocate a global edge number for
 380     // edge if it is used more than once!
 381     // Note that we are now sending the whole localShared to the master whereas
 382     // we only need the local count (i.e. the number of times a global edge is
 383     // used). But then this only gets done once so not too bothered about the
 384     // extra global communication.
 385
 386     EdgeMap<label> globalShared(nGlobalPoints());
 387
 388     if (Pstream::master())
 389     {
 390         label sharedEdgeI = 0;
 391
 392         // Merge my shared edges into the global list
 393         if (debug)
 394         {
 395             Pout<< "globalMeshData::calcSharedEdges : Merging in from proc0 : "
 396                 << localShared.size() << endl;
 397         }
 398         countSharedEdges(localShared, globalShared, sharedEdgeI);
 399
 400         // Receive data from slaves and insert
 401         if (Pstream::parRun())
 402         {
 403             for
 404             (
 405                 int slave=Pstream::firstSlave();
 406                 slave<=Pstream::lastSlave();
 407                 slave++
 408             )
 409             {
 410                 // Receive the edges using shared points from the slave.
 411                 IPstream fromSlave(Pstream::blocking, slave);
 412                 EdgeMap<labelList> procSharedEdges(fromSlave);
 413
 414                 if (debug)
 415                 {
 416                     Pout<< "globalMeshData::calcSharedEdges : "
 417                         << "Merging in from proc"
 418                         << Foam::name(slave) << " : " << procSharedEdges.size()
 419                         << endl;
 420                 }
 421                 countSharedEdges(procSharedEdges, globalShared, sharedEdgeI);
 422             }
 423         }
 424
 425         // Now our globalShared should have some edges with -1 as edge label
 426         // These were only used once so are not proper shared edges.
 427         // Remove them.
 428         {
 429             EdgeMap<label> oldSharedEdges(globalShared);
 430
 431             globalShared.clear();
 432
 433             forAllConstIter(EdgeMap<label>, oldSharedEdges, iter)
 434             {
 435                 if (iter() != -1)
 436                 {
 437                     globalShared.insert(iter.key(), iter());
 438                 }
 439             }
 440             if (debug)
 441             {
 442                 Pout<< "globalMeshData::calcSharedEdges : Filtered "
 443                     << oldSharedEdges.size()
 444                     << " down to " << globalShared.size() << endl;
 445             }
 446         }
 447
 448
 449         // Send back to slaves.
 450         if (Pstream::parRun())
 451         {
 452             for
 453             (
 454                 int slave=Pstream::firstSlave();
 455                 slave<=Pstream::lastSlave();
 456                 slave++
 457             )
 458             {
 459                 // Receive the edges using shared points from the slave.
 460                 OPstream toSlave(Pstream::blocking, slave);
 461                 toSlave << globalShared;
 462             }
 463         }
 464     }
 465     else
 466     {
 467         // Send local edges to master
 468         {
 469             OPstream toMaster(Pstream::blocking, Pstream::masterNo());
 470
 471             toMaster << localShared;
 472         }
 473         // Receive merged edges from master.
 474         {
 475             IPstream fromMaster(Pstream::blocking, Pstream::masterNo());
 476
 477             fromMaster >> globalShared;
 478         }
 479     }
 480
 481     // Now use the global shared edges list (globalShared) to classify my local
 482     // ones (localShared)
 483
 484     nGlobalEdges_ = globalShared.size();
 485
 486     DynamicList<label> dynSharedEdgeLabels(globalShared.size());
 487     DynamicList<label> dynSharedEdgeAddr(globalShared.size());
 488
 489     forAllConstIter(EdgeMap<labelList>, localShared, iter)
 490     {
 491         const edge& e = iter.key();
 492
 493         EdgeMap<label>::const_iterator edgeFnd = globalShared.find(e);
 494
 495         if (edgeFnd != globalShared.end())
 496         {
 497             // My local edge is indeed a shared one. Go through all local edge
 498             // labels with this point combination.
 499             const labelList& edgeLabels = iter();
 500
 501             forAll(edgeLabels, i)
 502             {
 503                 // Store label of local mesh edge
 504                 dynSharedEdgeLabels.append(edgeLabels[i]);
 505
 506                 // Store label of shared edge
 507                 dynSharedEdgeAddr.append(edgeFnd());
 508             }
 509         }
 510     }
 511
 512     sharedEdgeLabelsPtr_.reset(new labelList());
 513     labelList& sharedEdgeLabels = sharedEdgeLabelsPtr_();
 514     sharedEdgeLabels.transfer(dynSharedEdgeLabels);
 515
 516     sharedEdgeAddrPtr_.reset(new labelList());
 517     labelList& sharedEdgeAddr = sharedEdgeAddrPtr_();
 518     sharedEdgeAddr.transfer(dynSharedEdgeAddr);
 519
 520     if (debug)
 521     {
 522         Pout<< "globalMeshData : nGlobalEdges_:" << nGlobalEdges_ << nl
 523             << "globalMeshData : sharedEdgeLabels:" << sharedEdgeLabels.size()
 524             << nl
 525             << "globalMeshData : sharedEdgeAddr:" << sharedEdgeAddr.size()
 526             << endl;
 527     }
 528 }
 529
 530
 531 void Foam::globalMeshData::calcGlobalPointSlaves() const
 532 {
 533     if (debug)
 534     {
 535         Pout<< "globalMeshData::calcGlobalPointSlaves() :"
 536             << " calculating coupled master to slave point addressing."
 537             << endl;
 538     }
 539
 540     // Calculate connected points for master points.
 541     globalPoints globalData(mesh_, coupledPatch(), true, true);
 542
 543     globalPointSlavesPtr_.reset
 544     (
 545         new labelListList
 546         (
 547             globalData.pointPoints().xfer()
 548         )
 549     );
 550     globalPointTransformedSlavesPtr_.reset
 551     (
 552         new labelListList
 553         (
 554             globalData.transformedPointPoints().xfer()
 555         )
 556     );
 557
 558     globalPointSlavesMapPtr_.reset
 559     (
 560         new mapDistribute
 561         (
 562             globalData.map().xfer()
 563         )
 564     );
 565 }
 566
 567
 568 void Foam::globalMeshData::calcPointConnectivity
 569 (
 570     List<labelPairList>& allPointConnectivity
 571 ) const
 572 {
 573     const globalIndexAndTransform& transforms = globalTransforms();
 574     const labelListList& slaves = globalPointSlaves();
 575     const labelListList& transformedSlaves = globalPointTransformedSlaves();
 576
 577
 578     // Create field with my local data
 579     labelPairList myData(globalPointSlavesMap().constructSize());
 580     forAll(slaves, pointI)
 581     {
 582         myData[pointI] = globalIndexAndTransform::encode
 583         (
 584             Pstream::myProcNo(),
 585             pointI,
 586             transforms.nullTransformIndex()
 587         );
 588     }
 589     // Send over.
 590     globalPointSlavesMap().distribute(myData);
 591
 592
 593     // String of connected points with their transform
 594     allPointConnectivity.setSize(globalPointSlavesMap().constructSize());
 595     forAll(slaves, pointI)
 596     {
 597         // Reconstruct string of connected points
 598         const labelList& pSlaves = slaves[pointI];
 599         const labelList& pTransformSlaves = transformedSlaves[pointI];
 600         labelPairList& pConnectivity = allPointConnectivity[pointI];
 601         pConnectivity.setSize(1+pSlaves.size()+pTransformSlaves.size());
 602         label connI = 0;
 603
 604         // Add myself
 605         pConnectivity[connI++] = myData[pointI];
 606         // Add untransformed points
 607         forAll(pSlaves, i)
 608         {
 609             pConnectivity[connI++] = myData[pSlaves[i]];
 610         }
 611         // Add transformed points.
 612         forAll(pTransformSlaves, i)
 613         {
 614             // Get transform from index
 615             label transformI = globalPointSlavesMap().whichTransform
 616             (
 617                 pTransformSlaves[i]
 618             );
 619             // Add transform to connectivity
 620             const labelPair& n = myData[pTransformSlaves[i]];
 621             label procI = globalIndexAndTransform::processor(n);
 622             label index = globalIndexAndTransform::index(n);
 623             pConnectivity[connI++] = globalIndexAndTransform::encode
 624             (
 625                 procI,
 626                 index,
 627                 transformI
 628             );
 629         }
 630
 631         // Put back in slots
 632         forAll(pSlaves, i)
 633         {
 634             allPointConnectivity[pSlaves[i]] = pConnectivity;
 635         }
 636         forAll(pTransformSlaves, i)
 637         {
 638             allPointConnectivity[pTransformSlaves[i]] = pConnectivity;
 639         }
 640     }
 641     globalPointSlavesMap().reverseDistribute
 642     (
 643         allPointConnectivity.size(),
 644         allPointConnectivity
 645     );
 646 }
 647
 648
 649 void Foam::globalMeshData::calcGlobalPointEdges
 650 (
 651     labelListList& globalPointEdges,
 652     List<labelPairList>& globalPointPoints
 653 ) const
 654 {
 655     const edgeList& edges = coupledPatch().edges();
 656     const labelListList& pointEdges = coupledPatch().pointEdges();
 657     const globalIndex& globalEdgeNumbers = globalEdgeNumbering();
 658     const labelListList& slaves = globalPointSlaves();
 659     const labelListList& transformedSlaves = globalPointTransformedSlaves();
 660
 661     // Create local version
 662     globalPointEdges.setSize(globalPointSlavesMap().constructSize());
 663     globalPointPoints.setSize(globalPointSlavesMap().constructSize());
 664     forAll(pointEdges, pointI)
 665     {
 666         const labelList& pEdges = pointEdges[pointI];
 667         labelList& globalPEdges = globalPointEdges[pointI];
 668         globalPEdges.setSize(pEdges.size());
 669         forAll(pEdges, i)
 670         {
 671             globalPEdges[i] = globalEdgeNumbers.toGlobal(pEdges[i]);
 672         }
 673
 674         labelPairList& globalPPoints = globalPointPoints[pointI];
 675         globalPPoints.setSize(pEdges.size());
 676         forAll(pEdges, i)
 677         {
 678             label otherPointI = edges[pEdges[i]].otherVertex(pointI);
 679             globalPPoints[i] = globalIndexAndTransform::encode
 680             (
 681                 Pstream::myProcNo(),
 682                 otherPointI,
 683                 globalTransforms().nullTransformIndex()
 684             );
 685         }
 686     }
 687
 688     // Pull slave data to master. Dummy transform.
 689     globalPointSlavesMap().distribute(globalPointEdges);
 690     globalPointSlavesMap().distribute(globalPointPoints);
 691     // Add all pointEdges
 692     forAll(slaves, pointI)
 693     {
 694         const labelList& pSlaves = slaves[pointI];
 695         const labelList& pTransformSlaves = transformedSlaves[pointI];
 696
 697         label n = 0;
 698         forAll(pSlaves, i)
 699         {
 700             n += globalPointEdges[pSlaves[i]].size();
 701         }
 702         forAll(pTransformSlaves, i)
 703         {
 704             n += globalPointEdges[pTransformSlaves[i]].size();
 705         }
 706
 707         // Add all the point edges of the slaves to those of the (master) point
 708         {
 709             labelList& globalPEdges = globalPointEdges[pointI];
 710             label sz = globalPEdges.size();
 711             globalPEdges.setSize(sz+n);
 712             forAll(pSlaves, i)
 713             {
 714                 const labelList& otherData = globalPointEdges[pSlaves[i]];
 715                 forAll(otherData, j)
 716                 {
 717                     globalPEdges[sz++] = otherData[j];
 718                 }
 719             }
 720             forAll(pTransformSlaves, i)
 721             {
 722                 const labelList& otherData =
 723                     globalPointEdges[pTransformSlaves[i]];
 724                 forAll(otherData, j)
 725                 {
 726                     globalPEdges[sz++] = otherData[j];
 727                 }
 728             }
 729
 730             // Put back in slots
 731             forAll(pSlaves, i)
 732             {
 733                 globalPointEdges[pSlaves[i]] = globalPEdges;
 734             }
 735             forAll(pTransformSlaves, i)
 736             {
 737                 globalPointEdges[pTransformSlaves[i]] = globalPEdges;
 738             }
 739         }
 740
 741
 742         // Same for corresponding pointPoints
 743         {
 744             labelPairList& globalPPoints = globalPointPoints[pointI];
 745             label sz = globalPPoints.size();
 746             globalPPoints.setSize(sz + n);
 747
 748             // Add untransformed points
 749             forAll(pSlaves, i)
 750             {
 751                 const labelPairList& otherData = globalPointPoints[pSlaves[i]];
 752                 forAll(otherData, j)
 753                 {
 754                     globalPPoints[sz++] = otherData[j];
 755                 }
 756             }
 757             // Add transformed points.
 758             forAll(pTransformSlaves, i)
 759             {
 760                 // Get transform from index
 761                 label transformI = globalPointSlavesMap().whichTransform
 762                 (
 763                     pTransformSlaves[i]
 764                 );
 765
 766                 const labelPairList& otherData =
 767                     globalPointPoints[pTransformSlaves[i]];
 768                 forAll(otherData, j)
 769                 {
 770                     // Add transform to connectivity
 771                     const labelPair& n = otherData[j];
 772                     label procI = globalIndexAndTransform::processor(n);
 773                     label index = globalIndexAndTransform::index(n);
 774                     globalPPoints[sz++] = globalIndexAndTransform::encode
 775                     (
 776                         procI,
 777                         index,
 778                         transformI
 779                     );
 780                 }
 781             }
 782
 783             // Put back in slots
 784             forAll(pSlaves, i)
 785             {
 786                 globalPointPoints[pSlaves[i]] = globalPPoints;
 787             }
 788             forAll(pTransformSlaves, i)
 789             {
 790                 globalPointPoints[pTransformSlaves[i]] = globalPPoints;
 791             }
 792         }
 793     }
 794     // Push back
 795     globalPointSlavesMap().reverseDistribute
 796     (
 797         globalPointEdges.size(),
 798         globalPointEdges
 799     );
 800     // Push back
 801     globalPointSlavesMap().reverseDistribute
 802     (
 803         globalPointPoints.size(),
 804         globalPointPoints
 805     );
 806 }
 807
 808
 809 // Find transformation to take remotePoint to localPoint. Use info to find
 810 // the transforms.
 811 Foam::label Foam::globalMeshData::findTransform
 812 (
 813     const labelPairList& info,
 814     const labelPair& remotePoint,
 815     const label localPoint
 816 ) const
 817 {
 818     const label remoteProcI = globalIndexAndTransform::processor(remotePoint);
 819     const label remoteIndex = globalIndexAndTransform::index(remotePoint);
 820
 821     label remoteTransformI = -1;
 822     label localTransformI = -1;
 823     forAll(info, i)
 824     {
 825         label procI = globalIndexAndTransform::processor(info[i]);
 826         label pointI = globalIndexAndTransform::index(info[i]);
 827         label transformI = globalIndexAndTransform::transformIndex(info[i]);
 828
 829         if (procI == Pstream::myProcNo() && pointI == localPoint)
 830         {
 831             localTransformI = transformI;
 832             //Pout<< "For local :" << localPoint
 833             //    << " found transform:" << localTransformI
 834             //    << endl;
 835         }
 836         if (procI == remoteProcI && pointI == remoteIndex)
 837         {
 838             remoteTransformI = transformI;
 839             //Pout<< "For remote:" << remotePoint
 840             //    << " found transform:" << remoteTransformI
 841             //    << " at index:" << i
 842             //    << endl;
 843         }
 844     }
 845
 846     if (remoteTransformI == -1 || localTransformI == -1)
 847     {
 848         FatalErrorIn("globalMeshData::findTransform(..)")
 849             << "Problem. Cannot find " << remotePoint
 850             << " or " << localPoint << " in " << info
 851             << endl
 852             << "remoteTransformI:" << remoteTransformI << endl
 853             << "localTransformI:" << localTransformI
 854             << abort(FatalError);
 855     }
 856
 857     return globalTransforms().subtractTransformIndex
 858     (
 859         remoteTransformI,
 860         localTransformI
 861     );
 862 }
 863
 864
 865 void Foam::globalMeshData::calcGlobalEdgeSlaves() const
 866 {
 867     if (debug)
 868     {
 869         Pout<< "globalMeshData::calcGlobalEdgeSlaves() :"
 870             << " calculating coupled master to slave edge addressing." << endl;
 871     }
 872
 873     const edgeList& edges = coupledPatch().edges();
 874     const globalIndex& globalEdgeNumbers = globalEdgeNumbering();
 875
 876
 877     // The whole problem with deducting edge-connectivity from
 878     // point-connectivity is that one of the the endpoints might be
 879     // a local master but the other endpoint might not. So we first
 880     // need to make sure that all points know about connectivity and
 881     // the transformations.
 882
 883
 884     // 1. collect point connectivity - basically recreating globalPoints ouput.
 885     // All points will now have a string of points. The transforms are
 886     // in respect to the master.
 887     List<labelPairList> allPointConnectivity;
 888     calcPointConnectivity(allPointConnectivity);
 889
 890
 891     // 2. Get all pointEdges and pointPoints
 892     // Coupled point to global coupled edges and corresponding endpoint.
 893     labelListList globalPointEdges;
 894     List<labelPairList> globalPointPoints;
 895     calcGlobalPointEdges(globalPointEdges, globalPointPoints);
 896
 897
 898     // 3. Now all points have
 899     //      - all the connected points with original transform
 900     //      - all the connected global edges
 901
 902     // Now all we need to do is go through all the edges and check
 903     // both endpoints. If there is a edge between the two which is
 904     // produced by transforming both points in the same way it is a shared
 905     // edge.
 906
 907     // Collect strings of connected edges.
 908     List<labelPairList> allEdgeConnectivity(edges.size());
 909
 910     forAll(edges, edgeI)
 911     {
 912         const edge& e = edges[edgeI];
 913         const labelList& pEdges0 = globalPointEdges[e[0]];
 914         const labelPairList& pPoints0 = globalPointPoints[e[0]];
 915         const labelList& pEdges1 = globalPointEdges[e[1]];
 916         const labelPairList& pPoints1 = globalPointPoints[e[1]];
 917
 918         // Most edges will be size 2
 919         DynamicList<labelPair> eEdges(2);
 920         // Append myself.
 921         eEdges.append
 922         (
 923             globalIndexAndTransform::encode
 924             (
 925                 Pstream::myProcNo(),
 926                 edgeI,
 927                 globalTransforms().nullTransformIndex()
 928             )
 929         );
 930
 931         forAll(pEdges0, i)
 932         {
 933             forAll(pEdges1, j)
 934             {
 935                 if
 936                 (
 937                     pEdges0[i] == pEdges1[j]
 938                  && pEdges0[i] != globalEdgeNumbers.toGlobal(edgeI)
 939                 )
 940                 {
 941                     // Found a shared edge. Now check if the endpoints
 942                     // go through the same transformation.
 943                     // Local: e[0]    remote:pPoints1[j]
 944                     // Local: e[1]    remote:pPoints0[i]
 945
 946
 947                     // Find difference in transforms to go from point on remote
 948                     // edge (pPoints1[j]) to this point.
 949
 950                     label transform0 = findTransform
 951                     (
 952                         allPointConnectivity[e[0]],
 953                         pPoints1[j],
 954                         e[0]
 955                     );
 956                     label transform1 = findTransform
 957                     (
 958                         allPointConnectivity[e[1]],
 959                         pPoints0[i],
 960                         e[1]
 961                     );
 962
 963                     if (transform0 == transform1)
 964                     {
 965                         label procI = globalEdgeNumbers.whichProcID(pEdges0[i]);
 966                         eEdges.append
 967                         (
 968                             globalIndexAndTransform::encode
 969                             (
 970                                 procI,
 971                                 globalEdgeNumbers.toLocal(procI, pEdges0[i]),
 972                                 transform0
 973                             )
 974                         );
 975                     }
 976                 }
 977             }
 978         }
 979
 980         allEdgeConnectivity[edgeI].transfer(eEdges);
 981         sort(allEdgeConnectivity[edgeI], globalIndexAndTransform::less());
 982     }
 983
 984     // We now have - in allEdgeConnectivity - a list of edges which are shared
 985     // between multiple processors. Filter into non-transformed and transformed
 986     // connections.
 987
 988     globalEdgeSlavesPtr_.reset(new labelListList(edges.size()));
 989     labelListList& globalEdgeSlaves = globalEdgeSlavesPtr_();
 990     List<labelPairList> transformedEdges(edges.size());
 991     forAll(allEdgeConnectivity, edgeI)
 992     {
 993         const labelPairList& edgeInfo = allEdgeConnectivity[edgeI];
 994         if (edgeInfo.size() >= 2)
 995         {
 996             const labelPair& masterInfo = edgeInfo[0];
 997
 998             // Check if master edge (= first element (since sorted)) is me.
 999             if
1000             (
1001                 (
1002                     globalIndexAndTransform::processor(masterInfo)
1003                  == Pstream::myProcNo()
1004                 )
1005              && (globalIndexAndTransform::index(masterInfo) == edgeI)
1006             )
1007             {
1008                 // Sort into transformed and untransformed
1009                 labelList& eEdges = globalEdgeSlaves[edgeI];
1010                 eEdges.setSize(edgeInfo.size()-1);
1011
1012                 labelPairList& trafoEEdges = transformedEdges[edgeI];
1013                 trafoEEdges.setSize(edgeInfo.size()-1);
1014
1015                 label nonTransformI = 0;
1016                 label transformI = 0;
1017
1018                 for (label i = 1; i < edgeInfo.size(); i++)
1019                 {
1020                     const labelPair& info = edgeInfo[i];
1021                     label procI = globalIndexAndTransform::processor(info);
1022                     label index = globalIndexAndTransform::index(info);
1023                     label transform = globalIndexAndTransform::transformIndex
1024                     (
1025                         info
1026                     );
1027
1028                     if (transform == globalTransforms().nullTransformIndex())
1029                     {
1030                         eEdges[nonTransformI++] = globalEdgeNumbers.toGlobal
1031                         (
1032                             procI,
1033                             index
1034                         );
1035                     }
1036                     else
1037                     {
1038                         trafoEEdges[transformI++] = info;
1039                     }
1040                 }
1041
1042                 eEdges.setSize(nonTransformI);
1043                 trafoEEdges.setSize(transformI);
1044             }
1045         }
1046     }
1047
1048
1049     // Construct map
1050     globalEdgeTransformedSlavesPtr_.reset(new labelListList());
1051
1052     List<Map<label> > compactMap(Pstream::nProcs());
1053     globalEdgeSlavesMapPtr_.reset
1054     (
1055         new mapDistribute
1056         (
1057             globalEdgeNumbers,
1058             globalEdgeSlaves,
1059
1060             globalTransforms(),
1061             transformedEdges,
1062             globalEdgeTransformedSlavesPtr_(),
1063
1064             compactMap
1065         )
1066     );
1067
1068
1069     if (debug)
1070     {
1071         Pout<< "globalMeshData::calcGlobalEdgeSlaves() :"
1072             << " coupled edges:" << edges.size()
1073             << " additional coupled edges:"
1074             << globalEdgeSlavesMapPtr_().constructSize() - edges.size()
1075             << endl;
1076     }
1077 }
1078
1079
1080 void Foam::globalMeshData::calcGlobalEdgeOrientation() const
1081 {
1082     if (debug)
1083     {
1084         Pout<< "globalMeshData::calcGlobalEdgeOrientation() :"
1085             << " calculating edge orientation w.r.t. master edge." << endl;
1086     }
1087
1088     const globalIndex& globalPoints = globalPointNumbering();
1089
1090     // 1. Determine master point
1091     labelList masterPoint;
1092     {
1093         const mapDistribute& map = globalPointSlavesMap();
1094
1095         masterPoint.setSize(map.constructSize());
1096         masterPoint = labelMax;
1097
1098         for (label pointI = 0; pointI < coupledPatch().nPoints(); pointI++)
1099         {
1100             masterPoint[pointI] = globalPoints.toGlobal(pointI);
1101         }
1102         syncData
1103         (
1104             masterPoint,
1105             globalPointSlaves(),
1106             globalPointTransformedSlaves(),
1107             map,
1108             minEqOp<label>()
1109         );
1110     }
1111
1112     // Now all points should know who is master by comparing their global
1113     // pointID with the masterPointID. We now can use this information
1114     // to find the orientation of the master edge.
1115
1116     {
1117         const mapDistribute& map = globalEdgeSlavesMap();
1118         const labelListList& slaves = globalEdgeSlaves();
1119         const labelListList& transformedSlaves = globalEdgeTransformedSlaves();
1120
1121         // Distribute orientation of master edge (in masterPoint numbering)
1122         labelPairList masterEdgeVerts(map.constructSize());
1123         masterEdgeVerts = labelPair(labelMax, labelMax);
1124
1125         for (label edgeI = 0; edgeI < coupledPatch().nEdges(); edgeI++)
1126         {
1127             if
1128             (
1129                 (
1130                     slaves[edgeI].size()
1131                   + transformedSlaves[edgeI].size()
1132                 )
1133               > 0
1134             )
1135             {
1136                 // I am master. Fill in my masterPoint equivalent.
1137
1138                 const edge& e = coupledPatch().edges()[edgeI];
1139                 masterEdgeVerts[edgeI] = labelPair
1140                 (
1141                     masterPoint[e[0]],
1142                     masterPoint[e[1]]
1143                 );
1144             }
1145         }
1146         syncData
1147         (
1148             masterEdgeVerts,
1149             slaves,
1150             transformedSlaves,
1151             map,
1152             minEqOp<labelPair>()
1153         );
1154
1155         // Now check my edges on how they relate to the master's edgeVerts
1156         globalEdgeOrientationPtr_.reset
1157         (
1158             new PackedBoolList(coupledPatch().nEdges())
1159         );
1160         PackedBoolList& globalEdgeOrientation = globalEdgeOrientationPtr_();
1161
1162         forAll(coupledPatch().edges(), edgeI)
1163         {
1164             const edge& e = coupledPatch().edges()[edgeI];
1165             const labelPair masterE
1166             (
1167                 masterPoint[e[0]],
1168                 masterPoint[e[1]]
1169             );
1170
1171             label stat = labelPair::compare
1172             (
1173                 masterE,
1174                 masterEdgeVerts[edgeI]
1175             );
1176             if (stat == 0)
1177             {
1178                 FatalErrorIn
1179                 (
1180                     "globalMeshData::calcGlobalEdgeOrientation() const"
1181                 )   << "problem : my edge:" << e
1182                     << " in master points:" << masterE
1183                     << " v.s. masterEdgeVerts:" << masterEdgeVerts[edgeI]
1184                     << exit(FatalError);
1185             }
1186             else
1187             {
1188                 globalEdgeOrientation[edgeI] = (stat == 1);
1189             }
1190         }
1191     }
1192
1193     if (debug)
1194     {
1195         Pout<< "globalMeshData::calcGlobalEdgeOrientation() :"
1196             << " finished calculating edge orientation."
1197             << endl;
1198     }
1199 }
1200
1201
1202 // Calculate uncoupled boundary faces (without calculating
1203 // primitiveMesh::pointFaces())
1204 void Foam::globalMeshData::calcPointBoundaryFaces
1205 (
1206     labelListList& pointBoundaryFaces
1207 ) const
1208 {
1209     const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
1210     const Map<label>& meshPointMap = coupledPatch().meshPointMap();
1211
1212     // 1. Count
1213
1214     labelList nPointFaces(coupledPatch().nPoints(), 0);
1215
1216     forAll(bMesh, patchI)
1217     {
1218         const polyPatch& pp = bMesh[patchI];
1219
1220         if (!pp.coupled())
1221         {
1222             forAll(pp, i)
1223             {
1224                 const face& f = pp[i];
1225
1226                 forAll(f, fp)
1227                 {
1228                     Map<label>::const_iterator iter = meshPointMap.find
1229                     (
1230                         f[fp]
1231                     );
1232                     if (iter != meshPointMap.end())
1233                     {
1234                         nPointFaces[iter()]++;
1235                     }
1236                 }
1237             }
1238         }
1239     }
1240
1241
1242     // 2. Size
1243
1244     pointBoundaryFaces.setSize(coupledPatch().nPoints());
1245     forAll(nPointFaces, pointI)
1246     {
1247         pointBoundaryFaces[pointI].setSize(nPointFaces[pointI]);
1248     }
1249     nPointFaces = 0;
1250
1251
1252     // 3. Fill
1253
1254     forAll(bMesh, patchI)
1255     {
1256         const polyPatch& pp = bMesh[patchI];
1257
1258         if (!pp.coupled())
1259         {
1260             forAll(pp, i)
1261             {
1262                 const face& f = pp[i];
1263                 forAll(f, fp)
1264                 {
1265                     Map<label>::const_iterator iter = meshPointMap.find
1266                     (
1267                         f[fp]
1268                     );
1269                     if (iter != meshPointMap.end())
1270                     {
1271                         label bFaceI =
1272                              pp.start() + i - mesh_.nInternalFaces();
1273                         pointBoundaryFaces[iter()][nPointFaces[iter()]++] =
1274                             bFaceI;
1275                     }
1276                 }
1277             }
1278         }
1279     }
1280 }
1281
1282
1283 void Foam::globalMeshData::calcGlobalPointBoundaryFaces() const
1284 {
1285     if (debug)
1286     {
1287         Pout<< "globalMeshData::calcGlobalPointBoundaryFaces() :"
1288             << " calculating coupled point to boundary face addressing."
1289             << endl;
1290     }
1291
1292     // Construct local point to (uncoupled)boundaryfaces.
1293     labelListList pointBoundaryFaces;
1294     calcPointBoundaryFaces(pointBoundaryFaces);
1295
1296
1297     // Global indices for boundary faces
1298     globalBoundaryFaceNumberingPtr_.reset
1299     (
1300         new globalIndex(mesh_.nFaces()-mesh_.nInternalFaces())
1301     );
1302     globalIndex& globalIndices = globalBoundaryFaceNumberingPtr_();
1303
1304
1305     // Convert local boundary faces to global numbering
1306     globalPointBoundaryFacesPtr_.reset
1307     (
1308         new labelListList(globalPointSlavesMap().constructSize())
1309     );
1310     labelListList& globalPointBoundaryFaces = globalPointBoundaryFacesPtr_();
1311
1312     forAll(pointBoundaryFaces, pointI)
1313     {
1314         const labelList& bFaces = pointBoundaryFaces[pointI];
1315         labelList& globalFaces = globalPointBoundaryFaces[pointI];
1316         globalFaces.setSize(bFaces.size());
1317         forAll(bFaces, i)
1318         {
1319             globalFaces[i] = globalIndices.toGlobal(bFaces[i]);
1320         }
1321     }
1322
1323
1324     // Pull slave pointBoundaryFaces to master
1325     globalPointSlavesMap().distribute
1326     (
1327         globalPointBoundaryFaces,
1328         true    // put data on transformed points into correct slots
1329     );
1330
1331
1332     // Merge slave labels into master globalPointBoundaryFaces.
1333     // Split into untransformed and transformed values.
1334     const labelListList& pointSlaves = globalPointSlaves();
1335     const labelListList& pointTransformSlaves =
1336         globalPointTransformedSlaves();
1337
1338
1339     // Any faces coming in through transformation
1340     List<labelPairList> transformedFaces(pointSlaves.size());
1341
1342
1343     forAll(pointSlaves, pointI)
1344     {
1345         const labelList& slaves = pointSlaves[pointI];
1346         const labelList& transformedSlaves = pointTransformSlaves[pointI];
1347
1348         if (slaves.size() > 0)
1349         {
1350             labelList& myBFaces = globalPointBoundaryFaces[pointI];
1351             label sz = myBFaces.size();
1352
1353             // Count
1354             label n = 0;
1355             forAll(slaves, i)
1356             {
1357                 n += globalPointBoundaryFaces[slaves[i]].size();
1358             }
1359             // Fill
1360             myBFaces.setSize(sz+n);
1361             n = sz;
1362             forAll(slaves, i)
1363             {
1364                 const labelList& slaveBFaces =
1365                     globalPointBoundaryFaces[slaves[i]];
1366
1367                 // Add all slaveBFaces. Note that need to check for
1368                 // uniqueness only in case of cyclics.
1369
1370                 forAll(slaveBFaces, j)
1371                 {
1372                     label slave = slaveBFaces[j];
1373                     if (findIndex(SubList<label>(myBFaces, sz), slave) == -1)
1374                     {
1375                         myBFaces[n++] = slave;
1376                     }
1377                 }
1378             }
1379             myBFaces.setSize(n);
1380         }
1381
1382
1383         if (transformedSlaves.size() > 0)
1384         {
1385             const labelList& untrafoFaces = globalPointBoundaryFaces[pointI];
1386
1387             labelPairList& myBFaces = transformedFaces[pointI];
1388             label sz = myBFaces.size();
1389
1390             // Count
1391             label n = 0;
1392             forAll(transformedSlaves, i)
1393             {
1394                 n += globalPointBoundaryFaces[transformedSlaves[i]].size();
1395             }
1396             // Fill
1397             myBFaces.setSize(sz+n);
1398             n = sz;
1399             forAll(transformedSlaves, i)
1400             {
1401                 label transformI = globalPointSlavesMap().whichTransform
1402                 (
1403                     transformedSlaves[i]
1404                 );
1405
1406                 const labelList& slaveBFaces =
1407                     globalPointBoundaryFaces[transformedSlaves[i]];
1408
1409                 forAll(slaveBFaces, j)
1410                 {
1411                     label slave = slaveBFaces[j];
1412                     // Check that same face not already present untransformed
1413                     if (findIndex(untrafoFaces, slave)== -1)
1414                     {
1415                         label procI = globalIndices.whichProcID(slave);
1416                         label faceI = globalIndices.toLocal(procI, slave);
1417
1418                         myBFaces[n++] = globalIndexAndTransform::encode
1419                         (
1420                             procI,
1421                             faceI,
1422                             transformI
1423                         );
1424                     }
1425                 }
1426             }
1427             myBFaces.setSize(n);
1428         }
1429
1430
1431         if (slaves.size() + transformedSlaves.size() == 0)
1432         {
1433              globalPointBoundaryFaces[pointI].clear();
1434         }
1435     }
1436
1437     // Construct a map to get the face data directly
1438     List<Map<label> > compactMap(Pstream::nProcs());
1439
1440     globalPointTransformedBoundaryFacesPtr_.reset
1441     (
1442         new labelListList(transformedFaces.size())
1443     );
1444
1445     globalPointBoundaryFacesMapPtr_.reset
1446     (
1447         new mapDistribute
1448         (
1449             globalIndices,
1450             globalPointBoundaryFaces,
1451
1452             globalTransforms(),
1453             transformedFaces,
1454             globalPointTransformedBoundaryFacesPtr_(),
1455
1456             compactMap
1457         )
1458     );
1459     globalPointBoundaryFaces.setSize(coupledPatch().nPoints());
1460     globalPointTransformedBoundaryFacesPtr_().setSize(coupledPatch().nPoints());
1461
1462     if (debug)
1463     {
1464         Pout<< "globalMeshData::calcGlobalPointBoundaryFaces() :"
1465             << " coupled points:" << coupledPatch().nPoints()
1466             << " local boundary faces:" <<  globalIndices.localSize()
1467             << " additional coupled faces:"
1468             <<  globalPointBoundaryFacesMapPtr_().constructSize()
1469               - globalIndices.localSize()
1470             << endl;
1471     }
1472 }
1473
1474
1475 void Foam::globalMeshData::calcGlobalPointBoundaryCells() const
1476 {
1477     if (debug)
1478     {
1479         Pout<< "globalMeshData::calcGlobalPointBoundaryCells() :"
1480             << " calculating coupled point to boundary cell addressing."
1481             << endl;
1482     }
1483
1484     // Create map of boundary cells and point-cell addressing
1485     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1486
1487     label bCellI = 0;
1488     Map<label> meshCellMap(4*coupledPatch().nPoints());
1489     DynamicList<label> cellMap(meshCellMap.size());
1490
1491     // Create addressing for point to boundary cells (local)
1492     labelListList pointBoundaryCells(coupledPatch().nPoints());
1493
1494     forAll(coupledPatch().meshPoints(), pointI)
1495     {
1496         label meshPointI = coupledPatch().meshPoints()[pointI];
1497         const labelList& pCells = mesh_.pointCells(meshPointI);
1498
1499         labelList& bCells = pointBoundaryCells[pointI];
1500         bCells.setSize(pCells.size());
1501
1502         forAll(pCells, i)
1503         {
1504             label cellI = pCells[i];
1505             Map<label>::iterator fnd = meshCellMap.find(cellI);
1506
1507             if (fnd != meshCellMap.end())
1508             {
1509                 bCells[i] = fnd();
1510             }
1511             else
1512             {
1513                 meshCellMap.insert(cellI, bCellI);
1514                 cellMap.append(cellI);
1515                 bCells[i] = bCellI;
1516                 bCellI++;
1517             }
1518         }
1519     }
1520
1521
1522     boundaryCellsPtr_.reset(new labelList());
1523     labelList& boundaryCells = boundaryCellsPtr_();
1524     boundaryCells.transfer(cellMap.shrink());
1525
1526
1527     // Convert point-cells to global (boundary)cell numbers
1528     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1529
1530     globalBoundaryCellNumberingPtr_.reset
1531     (
1532         new globalIndex(boundaryCells.size())
1533     );
1534     globalIndex& globalIndices = globalBoundaryCellNumberingPtr_();
1535
1536
1537     globalPointBoundaryCellsPtr_.reset
1538     (
1539         new labelListList(globalPointSlavesMap().constructSize())
1540     );
1541     labelListList& globalPointBoundaryCells = globalPointBoundaryCellsPtr_();
1542
1543     forAll(pointBoundaryCells, pointI)
1544     {
1545         const labelList& pCells = pointBoundaryCells[pointI];
1546         labelList& globalCells = globalPointBoundaryCells[pointI];
1547         globalCells.setSize(pCells.size());
1548         forAll(pCells, i)
1549         {
1550             globalCells[i] = globalIndices.toGlobal(pCells[i]);
1551         }
1552     }
1553
1554
1555     // Pull slave pointBoundaryCells to master
1556     globalPointSlavesMap().distribute
1557     (
1558         globalPointBoundaryCells,
1559         true    // put data on transformed points into correct slots
1560     );
1561
1562
1563     // Merge slave labels into master globalPointBoundaryCells
1564     const labelListList& pointSlaves = globalPointSlaves();
1565     const labelListList& pointTransformSlaves =
1566         globalPointTransformedSlaves();
1567
1568     List<labelPairList> transformedCells(pointSlaves.size());
1569
1570
1571     forAll(pointSlaves, pointI)
1572     {
1573         const labelList& slaves = pointSlaves[pointI];
1574         const labelList& transformedSlaves = pointTransformSlaves[pointI];
1575
1576         if (slaves.size() > 0)
1577         {
1578             labelList& myBCells = globalPointBoundaryCells[pointI];
1579             label sz = myBCells.size();
1580
1581             // Count
1582             label n = 0;
1583             forAll(slaves, i)
1584             {
1585                 n += globalPointBoundaryCells[slaves[i]].size();
1586             }
1587             // Fill
1588             myBCells.setSize(sz+n);
1589             n = sz;
1590             forAll(slaves, i)
1591             {
1592                 const labelList& slaveBCells =
1593                     globalPointBoundaryCells[slaves[i]];
1594
1595                 // Add all slaveBCells. Note that need to check for
1596                 // uniqueness only in case of cyclics.
1597
1598                 forAll(slaveBCells, j)
1599                 {
1600                     label slave = slaveBCells[j];
1601                     if (findIndex(SubList<label>(myBCells, sz), slave) == -1)
1602                     {
1603                         myBCells[n++] = slave;
1604                     }
1605                 }
1606             }
1607             myBCells.setSize(n);
1608         }
1609
1610
1611         if (transformedSlaves.size() > 0)
1612         {
1613             const labelList& untrafoCells = globalPointBoundaryCells[pointI];
1614
1615             labelPairList& myBCells = transformedCells[pointI];
1616             label sz = myBCells.size();
1617
1618             // Count
1619             label n = 0;
1620             forAll(transformedSlaves, i)
1621             {
1622                 n += globalPointBoundaryCells[transformedSlaves[i]].size();
1623             }
1624             // Fill
1625             myBCells.setSize(sz+n);
1626             n = sz;
1627             forAll(transformedSlaves, i)
1628             {
1629                 label transformI = globalPointSlavesMap().whichTransform
1630                 (
1631                     transformedSlaves[i]
1632                 );
1633
1634                 const labelList& slaveBCells =
1635                     globalPointBoundaryCells[transformedSlaves[i]];
1636
1637                 forAll(slaveBCells, j)
1638                 {
1639                     label slave = slaveBCells[j];
1640
1641                     // Check that same cell not already present untransformed
1642                     if (findIndex(untrafoCells, slave)== -1)
1643                     {
1644                         label procI = globalIndices.whichProcID(slave);
1645                         label cellI = globalIndices.toLocal(procI, slave);
1646                         myBCells[n++] = globalIndexAndTransform::encode
1647                         (
1648                             procI,
1649                             cellI,
1650                             transformI
1651                         );
1652                     }
1653                 }
1654             }
1655             myBCells.setSize(n);
1656         }
1657
1658         if (slaves.size() + transformedSlaves.size() == 0)
1659         {
1660              globalPointBoundaryCells[pointI].clear();
1661         }
1662     }
1663
1664     // Construct a map to get the cell data directly
1665     List<Map<label> > compactMap(Pstream::nProcs());
1666
1667     globalPointTransformedBoundaryCellsPtr_.reset
1668     (
1669         new labelListList(transformedCells.size())
1670     );
1671
1672     globalPointBoundaryCellsMapPtr_.reset
1673     (
1674         new mapDistribute
1675         (
1676             globalIndices,
1677             globalPointBoundaryCells,
1678
1679             globalTransforms(),
1680             transformedCells,
1681             globalPointTransformedBoundaryCellsPtr_(),
1682
1683             compactMap
1684         )
1685     );
1686     globalPointBoundaryCells.setSize(coupledPatch().nPoints());
1687     globalPointTransformedBoundaryCellsPtr_().setSize(coupledPatch().nPoints());
1688
1689     if (debug)
1690     {
1691         Pout<< "globalMeshData::calcGlobalPointBoundaryCells() :"
1692             << " coupled points:" << coupledPatch().nPoints()
1693             << " local boundary cells:" <<  globalIndices.localSize()
1694             << " additional coupled cells:"
1695             <<  globalPointBoundaryCellsMapPtr_().constructSize()
1696               - globalIndices.localSize()
1697             << endl;
1698     }
1699 }
1700
1701
1702 void Foam::globalMeshData::calcGlobalCoPointSlaves() const
1703 {
1704     if (debug)
1705     {
1706         Pout<< "globalMeshData::calcGlobalCoPointSlaves() :"
1707             << " calculating coupled master to collocated"
1708             << " slave point addressing." << endl;
1709     }
1710
1711     // Calculate connected points for master points.
1712     globalPoints globalData(mesh_, coupledPatch(), true, false);
1713
1714     globalCoPointSlavesPtr_.reset
1715     (
1716         new labelListList
1717         (
1718             globalData.pointPoints().xfer()
1719         )
1720     );
1721     globalCoPointSlavesMapPtr_.reset
1722     (
1723         new mapDistribute
1724         (
1725             globalData.map().xfer()
1726         )
1727     );
1728
1729     if (debug)
1730     {
1731         Pout<< "globalMeshData::calcGlobalCoPointSlaves() :"
1732             << " finished calculating coupled master to collocated"
1733             << " slave point addressing." << endl;
1734     }
1735 }
1736
1737
1738 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
1739
1740 // Construct from polyMesh
1741 Foam::globalMeshData::globalMeshData(const polyMesh& mesh)
1742 :
1743     processorTopology(mesh.boundaryMesh()),
1744     mesh_(mesh),
1745     nTotalPoints_(-1),
1746     nTotalFaces_(-1),
1747     nTotalCells_(-1),
1748     processorPatches_(0),
1749     processorPatchIndices_(0),
1750     processorPatchNeighbours_(0),
1751     nGlobalPoints_(-1),
1752     sharedPointLabelsPtr_(NULL),
1753     sharedPointAddrPtr_(NULL),
1754     sharedPointGlobalLabelsPtr_(NULL),
1755     nGlobalEdges_(-1),
1756     sharedEdgeLabelsPtr_(NULL),
1757     sharedEdgeAddrPtr_(NULL)
1758 {
1759     updateMesh();
1760 }
1761
1762
1763 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
1764
1765 Foam::globalMeshData::~globalMeshData()
1766 {
1767     clearOut();
1768 }
1769
1770
1771 void Foam::globalMeshData::clearOut()
1772 {
1773     // Point
1774     nGlobalPoints_ = -1;
1775     sharedPointLabelsPtr_.clear();
1776     sharedPointAddrPtr_.clear();
1777     sharedPointGlobalLabelsPtr_.clear();
1778
1779     // Edge
1780     nGlobalEdges_ = -1;
1781     sharedEdgeLabelsPtr_.clear();
1782     sharedEdgeAddrPtr_.clear();
1783
1784     // Coupled patch
1785     coupledPatchPtr_.clear();
1786     coupledPatchMeshEdgesPtr_.clear();
1787     coupledPatchMeshEdgeMapPtr_.clear();
1788     globalTransformsPtr_.clear();
1789
1790     // Point
1791     globalPointNumberingPtr_.clear();
1792     globalPointSlavesPtr_.clear();
1793     globalPointTransformedSlavesPtr_.clear();
1794     globalPointSlavesMapPtr_.clear();
1795     // Edge
1796     globalEdgeNumberingPtr_.clear();
1797     globalEdgeSlavesPtr_.clear();
1798     globalEdgeTransformedSlavesPtr_.clear();
1799     globalEdgeOrientationPtr_.clear();
1800     globalEdgeSlavesMapPtr_.clear();
1801
1802     // Face
1803     globalBoundaryFaceNumberingPtr_.clear();
1804     globalPointBoundaryFacesPtr_.clear();
1805     globalPointTransformedBoundaryFacesPtr_.clear();
1806     globalPointBoundaryFacesMapPtr_.clear();
1807
1808     // Cell
1809     boundaryCellsPtr_.clear();
1810     globalBoundaryCellNumberingPtr_.clear();
1811     globalPointBoundaryCellsPtr_.clear();
1812     globalPointTransformedBoundaryCellsPtr_.clear();
1813     globalPointBoundaryCellsMapPtr_.clear();
1814
1815     // Other: collocated points
1816     globalCoPointSlavesPtr_.clear();
1817     globalCoPointSlavesMapPtr_.clear();
1818 }
1819
1820
1821 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1822
1823 // Return shared point global labels.
1824 const Foam::labelList& Foam::globalMeshData::sharedPointGlobalLabels() const
1825 {
1826     if (!sharedPointGlobalLabelsPtr_.valid())
1827     {
1828         sharedPointGlobalLabelsPtr_.reset
1829         (
1830             new labelList(sharedPointLabels().size())
1831         );
1832         labelList& sharedPointGlobalLabels = sharedPointGlobalLabelsPtr_();
1833
1834         IOobject addrHeader
1835         (
1836             "pointProcAddressing",
1837             mesh_.facesInstance()/mesh_.meshSubDir,
1838             mesh_,
1839             IOobject::MUST_READ
1840         );
1841
1842         if (addrHeader.headerOk())
1843         {
1844             // There is a pointProcAddressing file so use it to get labels
1845             // on the original mesh
1846             Pout<< "globalMeshData::sharedPointGlobalLabels : "
1847                 << "Reading pointProcAddressing" << endl;
1848
1849             labelIOList pointProcAddressing(addrHeader);
1850
1851             const labelList& pointLabels = sharedPointLabels();
1852
1853             forAll(pointLabels, i)
1854             {
1855                 // Get my mesh point
1856                 label pointI = pointLabels[i];
1857
1858                 // Map to mesh point of original mesh
1859                 sharedPointGlobalLabels[i] = pointProcAddressing[pointI];
1860             }
1861         }
1862         else
1863         {
1864             Pout<< "globalMeshData::sharedPointGlobalLabels :"
1865                 << " Setting pointProcAddressing to -1" << endl;
1866
1867             sharedPointGlobalLabels = -1;
1868         }
1869     }
1870     return sharedPointGlobalLabelsPtr_();
1871 }
1872
1873
1874 // Collect coordinates of shared points. (does parallel communication!)
1875 Foam::pointField Foam::globalMeshData::sharedPoints() const
1876 {
1877     // Get all processors to send their shared points to master.
1878     // (not very efficient)
1879
1880     pointField sharedPoints(nGlobalPoints());
1881     const labelList& pointAddr = sharedPointAddr();
1882     const labelList& pointLabels = sharedPointLabels();
1883
1884     if (Pstream::master())
1885     {
1886         // Master:
1887         // insert my own data first
1888         forAll(pointLabels, i)
1889         {
1890             label sharedPointI = pointAddr[i];
1891
1892             sharedPoints[sharedPointI] = mesh_.points()[pointLabels[i]];
1893         }
1894
1895         // Receive data from slaves and insert
1896         for
1897         (
1898             int slave=Pstream::firstSlave();
1899             slave<=Pstream::lastSlave();
1900             slave++
1901         )
1902         {
1903             IPstream fromSlave(Pstream::blocking, slave);
1904
1905             labelList nbrSharedPointAddr;
1906             pointField nbrSharedPoints;
1907             fromSlave >> nbrSharedPointAddr >> nbrSharedPoints;
1908
1909             forAll(nbrSharedPointAddr, i)
1910             {
1911                 label sharedPointI = nbrSharedPointAddr[i];
1912
1913                 sharedPoints[sharedPointI] = nbrSharedPoints[i];
1914             }
1915         }
1916
1917         // Send back
1918         for
1919         (
1920             int slave=Pstream::firstSlave();
1921             slave<=Pstream::lastSlave();
1922             slave++
1923         )
1924         {
1925             OPstream toSlave
1926             (
1927                 Pstream::blocking,
1928                 slave,
1929                 sharedPoints.size()*sizeof(vector::zero)
1930             );
1931             toSlave << sharedPoints;
1932         }
1933     }
1934     else
1935     {
1936         // Slave:
1937         // send points
1938         {
1939             OPstream toMaster(Pstream::blocking, Pstream::masterNo());
1940
1941             toMaster
1942                 << pointAddr
1943                 << UIndirectList<point>(mesh_.points(), pointLabels)();
1944         }
1945
1946         // Receive sharedPoints
1947         {
1948             IPstream fromMaster(Pstream::blocking, Pstream::masterNo());
1949             fromMaster >> sharedPoints;
1950         }
1951     }
1952
1953     return sharedPoints;
1954 }
1955
1956
1957 // Collect coordinates of shared points. (does parallel communication!)
1958 Foam::pointField Foam::globalMeshData::geometricSharedPoints() const
1959 {
1960     // Get coords of my shared points
1961     pointField sharedPoints(mesh_.points(), sharedPointLabels());
1962
1963     // Append from all processors
1964     combineReduce(sharedPoints, plusEqOp<pointField>());
1965
1966     // Merge tolerance
1967     scalar tolDim = matchTol_ * mesh_.bounds().mag();
1968
1969     // And see how many are unique
1970     labelList pMap;
1971     pointField mergedPoints;
1972
1973     Foam::mergePoints
1974     (
1975         sharedPoints,   // coordinates to merge
1976         tolDim,         // tolerance
1977         false,          // verbosity
1978         pMap,
1979         mergedPoints
1980     );
1981
1982     return mergedPoints;
1983 }
1984
1985
1986 Foam::label Foam::globalMeshData::nGlobalPoints() const
1987 {
1988     if (nGlobalPoints_ == -1)
1989     {
1990         calcSharedPoints();
1991     }
1992     return nGlobalPoints_;
1993 }
1994
1995
1996 const Foam::labelList& Foam::globalMeshData::sharedPointLabels() const
1997 {
1998     if (!sharedPointLabelsPtr_.valid())
1999     {
2000         calcSharedPoints();
2001     }
2002     return sharedPointLabelsPtr_();
2003 }
2004
2005
2006 const Foam::labelList& Foam::globalMeshData::sharedPointAddr() const
2007 {
2008     if (!sharedPointAddrPtr_.valid())
2009     {
2010         calcSharedPoints();
2011     }
2012     return sharedPointAddrPtr_();
2013 }
2014
2015
2016 Foam::label Foam::globalMeshData::nGlobalEdges() const
2017 {
2018     if (nGlobalEdges_ == -1)
2019     {
2020         calcSharedEdges();
2021     }
2022     return nGlobalEdges_;
2023 }
2024
2025
2026 const Foam::labelList& Foam::globalMeshData::sharedEdgeLabels() const
2027 {
2028     if (!sharedEdgeLabelsPtr_.valid())
2029     {
2030         calcSharedEdges();
2031     }
2032     return sharedEdgeLabelsPtr_();
2033 }
2034
2035
2036 const Foam::labelList& Foam::globalMeshData::sharedEdgeAddr() const
2037 {
2038     if (!sharedEdgeAddrPtr_.valid())
2039     {
2040         calcSharedEdges();
2041     }
2042     return sharedEdgeAddrPtr_();
2043 }
2044
2045
2046 const Foam::indirectPrimitivePatch& Foam::globalMeshData::coupledPatch() const
2047 {
2048     if (!coupledPatchPtr_.valid())
2049     {
2050         const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
2051
2052         label nCoupled = 0;
2053
2054         forAll(bMesh, patchI)
2055         {
2056             const polyPatch& pp = bMesh[patchI];
2057
2058             if (pp.coupled())
2059             {
2060                 nCoupled += pp.size();
2061             }
2062         }
2063         labelList coupledFaces(nCoupled);
2064         nCoupled = 0;
2065
2066         forAll(bMesh, patchI)
2067         {
2068             const polyPatch& pp = bMesh[patchI];
2069
2070             if (pp.coupled())
2071             {
2072                 label faceI = pp.start();
2073
2074                 forAll(pp, i)
2075                 {
2076                     coupledFaces[nCoupled++] = faceI++;
2077                 }
2078             }
2079         }
2080
2081         coupledPatchPtr_.reset
2082         (
2083             new indirectPrimitivePatch
2084             (
2085                 IndirectList<face>
2086                 (
2087                     mesh_.faces(),
2088                     coupledFaces
2089                 ),
2090                 mesh_.points()
2091             )
2092         );
2093
2094         if (debug)
2095         {
2096             Pout<< "globalMeshData::coupledPatch() :"
2097                 << " constructed  coupled faces patch:"
2098                 << "  faces:" << coupledPatchPtr_().size()
2099                 << "  points:" << coupledPatchPtr_().nPoints()
2100                 << endl;
2101         }
2102     }
2103     return coupledPatchPtr_();
2104 }
2105
2106
2107 const Foam::labelList& Foam::globalMeshData::coupledPatchMeshEdges() const
2108 {
2109     if (!coupledPatchMeshEdgesPtr_.valid())
2110     {
2111         coupledPatchMeshEdgesPtr_.reset
2112         (
2113             new labelList
2114             (
2115                 coupledPatch().meshEdges
2116                 (
2117                     mesh_.edges(),
2118                     mesh_.pointEdges()
2119                 )
2120             )
2121         );
2122     }
2123     return coupledPatchMeshEdgesPtr_();
2124 }
2125
2126
2127 const Foam::Map<Foam::label>& Foam::globalMeshData::coupledPatchMeshEdgeMap()
2128 const
2129 {
2130     if (!coupledPatchMeshEdgeMapPtr_.valid())
2131     {
2132         const labelList& me = coupledPatchMeshEdges();
2133
2134         coupledPatchMeshEdgeMapPtr_.reset(new Map<label>(2*me.size()));
2135         Map<label>& em = coupledPatchMeshEdgeMapPtr_();
2136
2137         forAll(me, i)
2138         {
2139             em.insert(me[i], i);
2140         }
2141     }
2142     return coupledPatchMeshEdgeMapPtr_();
2143 }
2144
2145
2146 const Foam::globalIndex& Foam::globalMeshData::globalPointNumbering() const
2147 {
2148     if (!globalPointNumberingPtr_.valid())
2149     {
2150         globalPointNumberingPtr_.reset
2151         (
2152             new globalIndex(coupledPatch().nPoints())
2153         );
2154     }
2155     return globalPointNumberingPtr_();
2156 }
2157
2158
2159 const Foam::globalIndexAndTransform&
2160 Foam::globalMeshData::globalTransforms() const
2161 {
2162     if (!globalTransformsPtr_.valid())
2163     {
2164         globalTransformsPtr_.reset(new globalIndexAndTransform(mesh_));
2165     }
2166     return globalTransformsPtr_();
2167 }
2168
2169
2170 const Foam::labelListList& Foam::globalMeshData::globalPointSlaves() const
2171 {
2172     if (!globalPointSlavesPtr_.valid())
2173     {
2174         calcGlobalPointSlaves();
2175     }
2176     return globalPointSlavesPtr_();
2177 }
2178
2179
2180 const Foam::labelListList& Foam::globalMeshData::globalPointTransformedSlaves()
2181 const
2182 {
2183     if (!globalPointTransformedSlavesPtr_.valid())
2184     {
2185         calcGlobalPointSlaves();
2186     }
2187     return globalPointTransformedSlavesPtr_();
2188 }
2189
2190
2191 const Foam::mapDistribute& Foam::globalMeshData::globalPointSlavesMap() const
2192 {
2193     if (!globalPointSlavesMapPtr_.valid())
2194     {
2195         calcGlobalPointSlaves();
2196     }
2197     return globalPointSlavesMapPtr_();
2198 }
2199
2200
2201 const Foam::globalIndex& Foam::globalMeshData::globalEdgeNumbering() const
2202 {
2203     if (!globalEdgeNumberingPtr_.valid())
2204     {
2205         globalEdgeNumberingPtr_.reset
2206         (
2207             new globalIndex(coupledPatch().nEdges())
2208         );
2209     }
2210     return globalEdgeNumberingPtr_();
2211 }
2212
2213
2214 const Foam::labelListList& Foam::globalMeshData::globalEdgeSlaves() const
2215 {
2216     if (!globalEdgeSlavesPtr_.valid())
2217     {
2218         calcGlobalEdgeSlaves();
2219     }
2220     return globalEdgeSlavesPtr_();
2221 }
2222
2223
2224 const Foam::labelListList& Foam::globalMeshData::globalEdgeTransformedSlaves()
2225 const
2226 {
2227     if (!globalEdgeTransformedSlavesPtr_.valid())
2228     {
2229         calcGlobalEdgeSlaves();
2230     }
2231     return globalEdgeTransformedSlavesPtr_();
2232 }
2233
2234
2235 const Foam::PackedBoolList& Foam::globalMeshData::globalEdgeOrientation() const
2236 {
2237     if (!globalEdgeOrientationPtr_.valid())
2238     {
2239         calcGlobalEdgeOrientation();
2240     }
2241     return globalEdgeOrientationPtr_();
2242 }
2243
2244
2245 const Foam::mapDistribute& Foam::globalMeshData::globalEdgeSlavesMap() const
2246 {
2247     if (!globalEdgeSlavesMapPtr_.valid())
2248     {
2249         calcGlobalEdgeSlaves();
2250     }
2251     return globalEdgeSlavesMapPtr_();
2252 }
2253
2254
2255 const Foam::globalIndex& Foam::globalMeshData::globalBoundaryFaceNumbering()
2256 const
2257 {
2258     if (!globalBoundaryFaceNumberingPtr_.valid())
2259     {
2260         calcGlobalPointBoundaryFaces();
2261     }
2262     return globalBoundaryFaceNumberingPtr_();
2263 }
2264
2265
2266 const Foam::labelListList& Foam::globalMeshData::globalPointBoundaryFaces()
2267 const
2268 {
2269     if (!globalPointBoundaryFacesPtr_.valid())
2270     {
2271         calcGlobalPointBoundaryFaces();
2272     }
2273     return globalPointBoundaryFacesPtr_();
2274 }
2275
2276
2277 const Foam::labelListList&
2278 Foam::globalMeshData::globalPointTransformedBoundaryFaces() const
2279 {
2280     if (!globalPointTransformedBoundaryFacesPtr_.valid())
2281     {
2282         calcGlobalPointBoundaryFaces();
2283     }
2284     return globalPointTransformedBoundaryFacesPtr_();
2285 }
2286
2287
2288 const Foam::mapDistribute& Foam::globalMeshData::globalPointBoundaryFacesMap()
2289 const
2290 {
2291     if (!globalPointBoundaryFacesMapPtr_.valid())
2292     {
2293         calcGlobalPointBoundaryFaces();
2294     }
2295     return globalPointBoundaryFacesMapPtr_();
2296 }
2297
2298
2299 const Foam::labelList& Foam::globalMeshData::boundaryCells() const
2300 {
2301     if (!boundaryCellsPtr_.valid())
2302     {
2303         calcGlobalPointBoundaryCells();
2304     }
2305     return boundaryCellsPtr_();
2306 }
2307
2308
2309 const Foam::globalIndex& Foam::globalMeshData::globalBoundaryCellNumbering()
2310 const
2311 {
2312     if (!globalBoundaryCellNumberingPtr_.valid())
2313     {
2314         calcGlobalPointBoundaryCells();
2315     }
2316     return globalBoundaryCellNumberingPtr_();
2317 }
2318
2319
2320 const Foam::labelListList& Foam::globalMeshData::globalPointBoundaryCells()
2321 const
2322 {
2323     if (!globalPointBoundaryCellsPtr_.valid())
2324     {
2325         calcGlobalPointBoundaryCells();
2326     }
2327     return globalPointBoundaryCellsPtr_();
2328 }
2329
2330
2331 const Foam::labelListList&
2332 Foam::globalMeshData::globalPointTransformedBoundaryCells() const
2333 {
2334     if (!globalPointTransformedBoundaryCellsPtr_.valid())
2335     {
2336         calcGlobalPointBoundaryCells();
2337     }
2338     return globalPointTransformedBoundaryCellsPtr_();
2339 }
2340
2341
2342 const Foam::mapDistribute& Foam::globalMeshData::globalPointBoundaryCellsMap()
2343 const
2344 {
2345     if (!globalPointBoundaryCellsMapPtr_.valid())
2346     {
2347         calcGlobalPointBoundaryCells();
2348     }
2349     return globalPointBoundaryCellsMapPtr_();
2350 }
2351
2352
2353 const Foam::labelListList& Foam::globalMeshData::globalCoPointSlaves() const
2354 {
2355     if (!globalCoPointSlavesPtr_.valid())
2356     {
2357         calcGlobalCoPointSlaves();
2358     }
2359     return globalCoPointSlavesPtr_();
2360 }
2361
2362
2363 const Foam::mapDistribute& Foam::globalMeshData::globalCoPointSlavesMap() const
2364 {
2365     if (!globalCoPointSlavesMapPtr_.valid())
2366     {
2367         calcGlobalCoPointSlaves();
2368     }
2369     return globalCoPointSlavesMapPtr_();
2370 }
2371
2372
2373 Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
2374 (
2375     labelList& pointToGlobal,
2376     labelList& uniquePoints
2377 ) const
2378 {
2379     const indirectPrimitivePatch& cpp = coupledPatch();
2380     const globalIndex& globalCoupledPoints = globalPointNumbering();
2381     // Use collocated only
2382     const labelListList& pointSlaves = globalCoPointSlaves();
2383     const mapDistribute& pointSlavesMap = globalCoPointSlavesMap();
2384
2385
2386     // Points are either
2387     // - master with slaves
2388     // - slave with a master
2389     // - other (since e.g. non-collocated cyclics not connected)
2390
2391     labelList masterGlobalPoint(cpp.nPoints(), -1);
2392     forAll(masterGlobalPoint, pointI)
2393     {
2394         const labelList& slavePoints = pointSlaves[pointI];
2395         if (slavePoints.size() > 0)
2396         {
2397             masterGlobalPoint[pointI] = globalCoupledPoints.toGlobal(pointI);
2398         }
2399     }
2400
2401     // Sync by taking max
2402     syncData
2403     (
2404         masterGlobalPoint,
2405         pointSlaves,
2406         labelListList(cpp.nPoints()),   // no transforms
2407         pointSlavesMap,
2408         maxEqOp<label>()
2409     );
2410
2411
2412     // 1. Count number of masters on my processor.
2413     label nMaster = 0;
2414     PackedBoolList isMaster(mesh_.nPoints(), 1);
2415     forAll(pointSlaves, pointI)
2416     {
2417         if (masterGlobalPoint[pointI] == -1)
2418         {
2419             // unconnected point (e.g. from separated cyclic)
2420             nMaster++;
2421         }
2422         else if
2423         (
2424             masterGlobalPoint[pointI]
2425          == globalCoupledPoints.toGlobal(pointI)
2426         )
2427         {
2428             // connected master
2429             nMaster++;
2430         }
2431         else
2432         {
2433             // connected slave point
2434             isMaster[cpp.meshPoints()[pointI]] = 0;
2435         }
2436     }
2437
2438     label myUniquePoints = mesh_.nPoints() - cpp.nPoints() + nMaster;
2439
2440     //Pout<< "Points :" << nl
2441     //    << "    mesh             : " << mesh_.nPoints() << nl
2442     //    << "    of which coupled : " << cpp.nPoints() << nl
2443     //    << "    of which master  : " << nMaster << nl
2444     //    << endl;
2445
2446
2447     // 2. Create global indexing for unique points.
2448     autoPtr<globalIndex> globalPointsPtr(new globalIndex(myUniquePoints));
2449
2450
2451     // 3. Assign global point numbers. Keep slaves unset.
2452     pointToGlobal.setSize(mesh_.nPoints());
2453     pointToGlobal = -1;
2454     uniquePoints.setSize(myUniquePoints);
2455     nMaster = 0;
2456
2457     forAll(isMaster, meshPointI)
2458     {
2459         if (isMaster[meshPointI])
2460         {
2461             pointToGlobal[meshPointI] = globalPointsPtr().toGlobal(nMaster);
2462             uniquePoints[nMaster] = meshPointI;
2463             nMaster++;
2464         }
2465     }
2466
2467
2468     // 4. Push global index for coupled points to slaves.
2469     {
2470         labelList masterToGlobal(pointSlavesMap.constructSize(), -1);
2471
2472         forAll(pointSlaves, pointI)
2473         {
2474             const labelList& slaves = pointSlaves[pointI];
2475
2476             if (slaves.size() > 0)
2477             {
2478                 // Duplicate master globalpoint into slave slots
2479                 label meshPointI = cpp.meshPoints()[pointI];
2480                 masterToGlobal[pointI] = pointToGlobal[meshPointI];
2481                 forAll(slaves, i)
2482                 {
2483                     masterToGlobal[slaves[i]] = masterToGlobal[pointI];
2484                 }
2485             }
2486         }
2487
2488         // Send back
2489         pointSlavesMap.reverseDistribute(cpp.nPoints(), masterToGlobal);
2490
2491         // On slave copy master index into overal map.
2492         forAll(pointSlaves, pointI)
2493         {
2494             label meshPointI = cpp.meshPoints()[pointI];
2495
2496             if (!isMaster[meshPointI])
2497             {
2498                 pointToGlobal[meshPointI] = masterToGlobal[pointI];
2499             }
2500         }
2501     }
2502
2503     return globalPointsPtr;
2504 }
2505
2506
2507 Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
2508 (
2509     const labelList& meshPoints,
2510     const Map<label>& meshPointMap,
2511     labelList& pointToGlobal,
2512     labelList& uniqueMeshPoints
2513 ) const
2514 {
2515     const indirectPrimitivePatch& cpp = coupledPatch();
2516     const labelListList& pointSlaves = globalCoPointSlaves();
2517     const mapDistribute& pointSlavesMap = globalCoPointSlavesMap();
2518
2519
2520     // The patch points come in two variants:
2521     // - not on a coupled patch so guaranteed unique
2522     // - on a coupled patch
2523     // If the point is on a coupled patch the problem is that the
2524     // master-slave structure (globalPointSlaves etc.) assigns one of the
2525     // coupled points to be the master but this master point is not
2526     // necessarily on the patch itself! (it might just be connected to the
2527     // patch point via coupled patches).
2528
2529
2530     // Determine mapping:
2531     // - from patch point to coupled point (or -1)
2532     // - from coupled point to global patch point
2533     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2534
2535     globalIndex globalPPoints(meshPoints.size());
2536
2537     labelList patchToCoupled(meshPoints.size(), -1);
2538     label nCoupled = 0;
2539     labelList coupledToGlobalPatch(pointSlavesMap.constructSize(), -1);
2540
2541     // Note: loop over patch since usually smaller
2542     forAll(meshPoints, patchPointI)
2543     {
2544         label meshPointI = meshPoints[patchPointI];
2545
2546         Map<label>::const_iterator iter = cpp.meshPointMap().find(meshPointI);
2547
2548         if (iter != cpp.meshPointMap().end())
2549         {
2550             patchToCoupled[patchPointI] = iter();
2551             coupledToGlobalPatch[iter()] = globalPPoints.toGlobal(patchPointI);
2552             nCoupled++;
2553         }
2554     }
2555
2556     //Pout<< "Patch:" << nl
2557     //    << "    points:" << meshPoints.size() << nl
2558     //    << "    of which on coupled patch:" << nCoupled << endl;
2559
2560
2561     // Determine master of connected points
2562     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2563     // Problem is that the coupled master might not be on the patch. So
2564     // work out the best patch-point master from all connected points.
2565     // - if the coupled master is on the patch it becomes the patch-point master
2566     // - else the slave with the lowest numbered patch point label
2567
2568     // Get all data on master
2569     pointSlavesMap.distribute(coupledToGlobalPatch);
2570     forAll(pointSlaves, coupledPointI)
2571     {
2572         const labelList& slaves = pointSlaves[coupledPointI];
2573
2574         if (slaves.size() > 0)
2575         {
2576             // I am master. What is the best candidate for patch-point master
2577             label masterI = labelMax;
2578             if (coupledToGlobalPatch[coupledPointI] != -1)
2579             {
2580                 // I am master and on the coupled patch. Use me.
2581                 masterI = coupledToGlobalPatch[coupledPointI];
2582             }
2583             else
2584             {
2585                 // Get min of slaves as master.
2586                 forAll(slaves, i)
2587                 {
2588                     label slavePp = coupledToGlobalPatch[slaves[i]];
2589                     if (slavePp != -1 && slavePp < masterI)
2590                     {
2591                         masterI = slavePp;
2592                     }
2593                 }
2594             }
2595
2596             if (masterI != labelMax)
2597             {
2598                 // Push back
2599                 coupledToGlobalPatch[coupledPointI] = masterI;
2600                 forAll(slaves, i)
2601                 {
2602                     coupledToGlobalPatch[slaves[i]] = masterI;
2603                 }
2604             }
2605         }
2606     }
2607     pointSlavesMap.reverseDistribute(cpp.nPoints(), coupledToGlobalPatch);
2608
2609
2610     // Generate compact numbering for master points
2611     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2612     // Now coupledToGlobalPatch is the globalIndex of the master point.
2613     // Now every processor can check whether they hold it and generate a
2614     // compact numbering.
2615
2616     label nMasters = 0;
2617     forAll(meshPoints, patchPointI)
2618     {
2619         if (patchToCoupled[patchPointI] == -1)
2620         {
2621             nMasters++;
2622         }
2623         else
2624         {
2625             label coupledPointI = patchToCoupled[patchPointI];
2626             if
2627             (
2628                 globalPPoints.toGlobal(patchPointI)
2629              == coupledToGlobalPatch[coupledPointI]
2630             )
2631             {
2632                 // I am the master
2633                 nMasters++;
2634             }
2635         }
2636     }
2637
2638     autoPtr<globalIndex> globalPointsPtr(new globalIndex(nMasters));
2639
2640     //Pout<< "Patch:" << nl
2641     //    << "    points:" << meshPoints.size() << nl
2642     //    << "    of which on coupled patch:" << nCoupled << nl
2643     //    << "    of which master:" << nMasters << endl;
2644
2645
2646
2647     // Push back compact numbering for master point onto slaves
2648     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2649
2650     pointToGlobal.setSize(meshPoints.size());
2651     pointToGlobal = -1;
2652     uniqueMeshPoints.setSize(nMasters);
2653
2654     // Sync master in global point numbering so all know the master point.
2655     // Initialise globalMaster to be -1 except at a globalMaster.
2656     labelList globalMaster(cpp.nPoints(), -1);
2657
2658     nMasters = 0;
2659     forAll(meshPoints, patchPointI)
2660     {
2661         if (patchToCoupled[patchPointI] == -1)
2662         {
2663             uniqueMeshPoints[nMasters++] = meshPoints[patchPointI];
2664         }
2665         else
2666         {
2667             label coupledPointI = patchToCoupled[patchPointI];
2668             if
2669             (
2670                 globalPPoints.toGlobal(patchPointI)
2671              == coupledToGlobalPatch[coupledPointI]
2672             )
2673             {
2674                 globalMaster[coupledPointI] =
2675                     globalPointsPtr().toGlobal(nMasters);
2676                 uniqueMeshPoints[nMasters++] = meshPoints[patchPointI];
2677             }
2678         }
2679     }
2680
2681
2682     // Sync by taking max
2683     syncData
2684     (
2685         globalMaster,
2686         pointSlaves,
2687         labelListList(cpp.nPoints()),   // no transforms
2688         pointSlavesMap,
2689         maxEqOp<label>()
2690     );
2691
2692
2693     // Now everyone has the master point in globalPointsPtr numbering. Fill
2694     // in the pointToGlobal map.
2695     nMasters = 0;
2696     forAll(meshPoints, patchPointI)
2697     {
2698         if (patchToCoupled[patchPointI] == -1)
2699         {
2700             pointToGlobal[patchPointI] = globalPointsPtr().toGlobal(nMasters++);
2701         }
2702         else
2703         {
2704             label coupledPointI = patchToCoupled[patchPointI];
2705             pointToGlobal[patchPointI] = globalMaster[coupledPointI];
2706
2707             if
2708             (
2709                 globalPPoints.toGlobal(patchPointI)
2710              == coupledToGlobalPatch[coupledPointI]
2711             )
2712             {
2713                 nMasters++;
2714             }
2715         }
2716     }
2717
2718     return globalPointsPtr;
2719 }
2720
2721
2722 void Foam::globalMeshData::movePoints(const pointField& newPoints)
2723 {
2724     // Topology does not change and we don't store any geometry so nothing
2725     // needs to be done.
2726     // Only global transformations might change but this is not really
2727     // supported.
2728 }
2729
2730
2731 // Update all data after morph
2732 void Foam::globalMeshData::updateMesh()
2733 {
2734     // Clear out old data
2735     clearOut();
2736
2737     // Do processor patch addressing
2738     initProcAddr();
2739
2740     scalar tolDim = matchTol_ * mesh_.bounds().mag();
2741
2742     if (debug)
2743     {
2744         Pout<< "globalMeshData : merge dist:" << tolDim << endl;
2745     }
2746
2747     // Total number of faces.
2748     nTotalFaces_ = returnReduce(mesh_.nFaces(), sumOp<label>());
2749
2750     if (debug)
2751     {
2752         Pout<< "globalMeshData : nTotalFaces_:" << nTotalFaces_ << endl;
2753     }
2754
2755     nTotalCells_ = returnReduce(mesh_.nCells(), sumOp<label>());
2756
2757     if (debug)
2758     {
2759         Pout<< "globalMeshData : nTotalCells_:" << nTotalCells_ << endl;
2760     }
2761
2762     nTotalPoints_ = returnReduce(mesh_.nPoints(), sumOp<label>());
2763
2764     if (debug)
2765     {
2766         Pout<< "globalMeshData : nTotalPoints_:" << nTotalPoints_ << endl;
2767     }
2768 }
2769
2770
2771 // ************************************************************************* //