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