src/foam/meshes/polyMesh/globalMeshData/globalMeshData.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "globalMeshData.H"
  27 #include "foamTime.H"
  28 #include "Pstream.H"
  29 #include "PstreamCombineReduceOps.H"
  30 #include "processorPolyPatch.H"
  31 #include "demandDrivenData.H"
  32 #include "globalPoints.H"
  33 //#include "geomGlobalPoints.H"
  34 #include "labelIOList.H"
  35 #include "PackedList.H"
  36 #include "mergePoints.H"
  37 #include "matchPoints.H"
  38 #include "OFstream.H"
  39
  40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  41
  42 defineTypeNameAndDebug(Foam::globalMeshData, 0);
  43
  44 // Geometric matching tolerance. Factor of mesh bounding box.
  45 const Foam::debug::tolerancesSwitch
  46 Foam::globalMeshData::matchTol_
  47 (
  48     "globalMeshDataMatchTol",
  49     1e-8
  50 );
  51
  52
  53 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  54
  55 // Collect processor patch addressing.
  56 void Foam::globalMeshData::initProcAddr()
  57 {
  58     processorPatchIndices_.setSize(mesh_.boundaryMesh().size());
  59     processorPatchIndices_ = -1;
  60
  61     processorPatchNeighbours_.setSize(mesh_.boundaryMesh().size());
  62     processorPatchNeighbours_ = -1;
  63
  64     // Construct processor patch indexing. processorPatchNeighbours_ only
  65     // set if running in parallel!
  66     processorPatches_.setSize(mesh_.boundaryMesh().size());
  67
  68     label nNeighbours = 0;
  69
  70     forAll (mesh_.boundaryMesh(), patchi)
  71     {
  72         if (isA<processorPolyPatch>(mesh_.boundaryMesh()[patchi]))
  73         {
  74             processorPatches_[nNeighbours] = patchi;
  75             processorPatchIndices_[patchi] = nNeighbours++;
  76         }
  77     }
  78     processorPatches_.setSize(nNeighbours);
  79
  80     if (Pstream::parRun())
  81     {
  82         // Send indices of my processor patches to my neighbours
  83         forAll (processorPatches_, i)
  84         {
  85             label patchi = processorPatches_[i];
  86
  87             OPstream toNeighbour
  88             (
  89                 Pstream::blocking,
  90                 refCast<const processorPolyPatch>
  91                 (
  92                     mesh_.boundaryMesh()[patchi]
  93                 ).neighbProcNo()
  94             );
  95
  96             toNeighbour << processorPatchIndices_[patchi];
  97         }
  98
  99         forAll(processorPatches_, i)
 100         {
 101             label patchi = processorPatches_[i];
 102
 103             IPstream fromNeighbour
 104             (
 105                 Pstream::blocking,
 106                 refCast<const processorPolyPatch>
 107                 (
 108                     mesh_.boundaryMesh()[patchi]
 109                 ).neighbProcNo()
 110             );
 111
 112             fromNeighbour >> processorPatchNeighbours_[patchi];
 113         }
 114     }
 115 }
 116
 117
 118 // Given information about locally used edges allocate global shared edges.
 119 void Foam::globalMeshData::countSharedEdges
 120 (
 121     const HashTable<labelList, edge, Hash<edge> >& procSharedEdges,
 122     HashTable<label, edge, Hash<edge> >& globalShared,
 123     label& sharedEdgeI
 124 )
 125 {
 126     // Count occurrences of procSharedEdges in global shared edges table.
 127     for
 128     (
 129         HashTable<labelList, edge, Hash<edge> >::const_iterator iter =
 130             procSharedEdges.begin();
 131         iter != procSharedEdges.end();
 132         ++iter
 133     )
 134     {
 135         const edge& e = iter.key();
 136
 137         HashTable<label, edge, Hash<edge> >::iterator globalFnd =
 138             globalShared.find(e);
 139
 140         if (globalFnd == globalShared.end())
 141         {
 142             // First time occurrence of this edge. Check how many we are adding
 143             if (iter().size() == 1)
 144             {
 145                 // Only one edge. Mark with special value.
 146                 globalShared.insert(e, -1);
 147             }
 148             else
 149             {
 150                 // Edge used more than once (even by local shared edges alone)
 151                 // so allocate proper shared edge label.
 152                 globalShared.insert(e, sharedEdgeI++);
 153             }
 154         }
 155         else
 156         {
 157             if (globalFnd() == -1)
 158             {
 159                 // Second time occurence of this edge. Assign proper
 160                 // edge label.
 161                 globalFnd() = sharedEdgeI++;
 162             }
 163         }
 164     }
 165 }
 166
 167
 168 // Shared edges are shared between multiple processors. By their nature both
 169 // of their endpoints are shared points. (but not all edges using two shared
 170 // points are shared edges! There might e.g. be an edge between two unrelated
 171 // clusters of shared points)
 172 void Foam::globalMeshData::calcSharedEdges() const
 173 {
 174     if (nGlobalEdges_ != -1 || sharedEdgeLabelsPtr_ || sharedEdgeAddrPtr_)
 175     {
 176         FatalErrorIn("globalMeshData::calcSharedEdges()")
 177             << "Shared edge addressing already done" << abort(FatalError);
 178     }
 179
 180
 181     const labelList& sharedPtAddr = sharedPointAddr();
 182     const labelList& sharedPtLabels = sharedPointLabels();
 183
 184     // Since don't want to construct pointEdges for whole mesh create
 185     // Map for all shared points.
 186     Map<label> meshToShared(2*sharedPtLabels.size());
 187     forAll(sharedPtLabels, i)
 188     {
 189         meshToShared.insert(sharedPtLabels[i], i);
 190     }
 191
 192     // Find edges using shared points. Store correspondence to local edge
 193     // numbering. Note that multiple local edges can have the same shared
 194     // points! (for cyclics or separated processor patches)
 195     HashTable<labelList, edge, Hash<edge> > localShared
 196     (
 197         2*sharedPtAddr.size()
 198     );
 199
 200     const edgeList& edges = mesh_.edges();
 201
 202     forAll(edges, edgeI)
 203     {
 204         const edge& e = edges[edgeI];
 205
 206         Map<label>::const_iterator e0Fnd = meshToShared.find(e[0]);
 207
 208         if (e0Fnd != meshToShared.end())
 209         {
 210             Map<label>::const_iterator e1Fnd = meshToShared.find(e[1]);
 211
 212             if (e1Fnd != meshToShared.end())
 213             {
 214                 // Found edge which uses shared points. Probably shared.
 215
 216                 // Construct the edge in shared points (or rather global
 217                 //  indices of the shared points)
 218                 edge sharedEdge
 219                 (
 220                     sharedPtAddr[e0Fnd()],
 221                     sharedPtAddr[e1Fnd()]
 222                 );
 223
 224                 HashTable<labelList, edge, Hash<edge> >::iterator iter =
 225                     localShared.find(sharedEdge);
 226
 227                 if (iter == localShared.end())
 228                 {
 229                     // First occurrence of this point combination. Store.
 230                     localShared.insert(sharedEdge, labelList(1, edgeI));
 231                 }
 232                 else
 233                 {
 234                     // Add this edge to list of edge labels.
 235                     labelList& edgeLabels = iter();
 236
 237                     label sz = edgeLabels.size();
 238                     edgeLabels.setSize(sz+1);
 239                     edgeLabels[sz] = edgeI;
 240                 }
 241             }
 242         }
 243     }
 244
 245
 246     // Now we have a table on every processors which gives its edges which use
 247     // shared points. Send this all to the master and have it allocate
 248     // global edge numbers for it. But only allocate a global edge number for
 249     // edge if it is used more than once!
 250     // Note that we are now sending the whole localShared to the master whereas
 251     // we only need the local count (i.e. the number of times a global edge is
 252     // used). But then this only gets done once so not too bothered about the
 253     // extra global communication.
 254
 255     HashTable<label, edge, Hash<edge> > globalShared(nGlobalPoints());
 256
 257     if (Pstream::master())
 258     {
 259         label sharedEdgeI = 0;
 260
 261         // Merge my shared edges into the global list
 262         if (debug)
 263         {
 264             Pout<< "globalMeshData::calcSharedEdges : Merging in from proc0 : "
 265                 << localShared.size() << endl;
 266         }
 267         countSharedEdges(localShared, globalShared, sharedEdgeI);
 268
 269         // Receive data from slaves and insert
 270         if (Pstream::parRun())
 271         {
 272             for
 273             (
 274                 int slave=Pstream::firstSlave();
 275                 slave<=Pstream::lastSlave();
 276                 slave++
 277             )
 278             {
 279                 // Receive the edges using shared points from the slave.
 280                 IPstream fromSlave(Pstream::blocking, slave);
 281                 HashTable<labelList, edge, Hash<edge> > procSharedEdges
 282                 (
 283                     fromSlave
 284                 );
 285
 286                 if (debug)
 287                 {
 288                     Pout<< "globalMeshData::calcSharedEdges : "
 289                         << "Merging in from proc"
 290                         << Foam::name(slave) << " : " << procSharedEdges.size()
 291                         << endl;
 292                 }
 293                 countSharedEdges(procSharedEdges, globalShared, sharedEdgeI);
 294             }
 295         }
 296
 297         // Now our globalShared should have some edges with -1 as edge label
 298         // These were only used once so are not proper shared edges.
 299         // Remove them.
 300         {
 301             HashTable<label, edge, Hash<edge> > oldSharedEdges(globalShared);
 302
 303             globalShared.clear();
 304
 305             for
 306             (
 307                 HashTable<label, edge, Hash<edge> >::const_iterator iter =
 308                     oldSharedEdges.begin();
 309                 iter != oldSharedEdges.end();
 310                 ++iter
 311             )
 312             {
 313                 if (iter() != -1)
 314                 {
 315                     globalShared.insert(iter.key(), iter());
 316                 }
 317             }
 318             if (debug)
 319             {
 320                 Pout<< "globalMeshData::calcSharedEdges : Filtered "
 321                     << oldSharedEdges.size()
 322                     << " down to " << globalShared.size() << endl;
 323             }
 324         }
 325
 326
 327         // Send back to slaves.
 328         if (Pstream::parRun())
 329         {
 330             for
 331             (
 332                 int slave=Pstream::firstSlave();
 333                 slave<=Pstream::lastSlave();
 334                 slave++
 335             )
 336             {
 337                 // Receive the edges using shared points from the slave.
 338                 OPstream toSlave(Pstream::blocking, slave);
 339                 toSlave << globalShared;
 340             }
 341         }
 342     }
 343     else
 344     {
 345         // Send local edges to master
 346         {
 347             OPstream toMaster(Pstream::blocking, Pstream::masterNo());
 348
 349             toMaster << localShared;
 350         }
 351         // Receive merged edges from master.
 352         {
 353             IPstream fromMaster(Pstream::blocking, Pstream::masterNo());
 354
 355             fromMaster >> globalShared;
 356         }
 357     }
 358
 359     // Now use the global shared edges list (globalShared) to classify my local
 360     // ones (localShared)
 361
 362     nGlobalEdges_ = globalShared.size();
 363
 364     DynamicList<label> dynSharedEdgeLabels(globalShared.size());
 365     DynamicList<label> dynSharedEdgeAddr(globalShared.size());
 366
 367     for
 368     (
 369         HashTable<labelList, edge, Hash<edge> >::const_iterator iter =
 370             localShared.begin();
 371         iter != localShared.end();
 372         ++iter
 373     )
 374     {
 375         const edge& e = iter.key();
 376
 377         HashTable<label, edge, Hash<edge> >::const_iterator edgeFnd =
 378             globalShared.find(e);
 379
 380         if (edgeFnd != globalShared.end())
 381         {
 382             // My local edge is indeed a shared one. Go through all local edge
 383             // labels with this point combination.
 384             const labelList& edgeLabels = iter();
 385
 386             forAll(edgeLabels, i)
 387             {
 388                 // Store label of local mesh edge
 389                 dynSharedEdgeLabels.append(edgeLabels[i]);
 390
 391                 // Store label of shared edge
 392                 dynSharedEdgeAddr.append(edgeFnd());
 393             }
 394         }
 395     }
 396
 397     sharedEdgeLabelsPtr_ = new labelList();
 398     labelList& sharedEdgeLabels = *sharedEdgeLabelsPtr_;
 399     sharedEdgeLabels.transfer(dynSharedEdgeLabels);
 400
 401     sharedEdgeAddrPtr_ = new labelList();
 402     labelList& sharedEdgeAddr = *sharedEdgeAddrPtr_;
 403     sharedEdgeAddr.transfer(dynSharedEdgeAddr);
 404
 405     if (debug)
 406     {
 407         Pout<< "globalMeshData : nGlobalEdges_:" << nGlobalEdges_ << nl
 408             << "globalMeshData : sharedEdgeLabels:" << sharedEdgeLabels.size()
 409             << nl
 410             << "globalMeshData : sharedEdgeAddr:" << sharedEdgeAddr.size()
 411             << endl;
 412     }
 413 }
 414
 415
 416 // Helper function to count coincident faces. This part used to be
 417 // in updateMesh but I've had to move it to a separate function
 418 // because of aliasing optimisation errors in icc9.1 on the
 419 // Itanium.
 420 Foam::label Foam::globalMeshData::countCoincidentFaces
 421 (
 422     const scalar tolDim,
 423     const vectorField& separationDist
 424 )
 425 {
 426     label nCoincident = 0;
 427
 428     forAll(separationDist, faceI)
 429     {
 430         if (mag(separationDist[faceI]) < tolDim)
 431         {
 432             // Faces coincide
 433             nCoincident++;
 434         }
 435     }
 436     return nCoincident;
 437 }
 438
 439
 440 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 441
 442 // Construct from polyMesh
 443 Foam::globalMeshData::globalMeshData(const polyMesh& mesh)
 444 :
 445     processorTopology(mesh.boundaryMesh()),
 446     mesh_(mesh),
 447     bb_(vector::zero, vector::zero),
 448     nTotalPoints_(-1),
 449     nTotalFaces_(-1),
 450     nTotalCells_(-1),
 451     processorPatches_(0),
 452     processorPatchIndices_(0),
 453     processorPatchNeighbours_(0),
 454     nGlobalPoints_(-1),
 455     sharedPointLabels_(0),
 456     sharedPointAddr_(0),
 457     sharedPointGlobalLabelsPtr_(NULL),
 458     nGlobalEdges_(-1),
 459     sharedEdgeLabelsPtr_(NULL),
 460     sharedEdgeAddrPtr_(NULL)
 461 {
 462     updateMesh();
 463 }
 464
 465
 466 // Read constructor given IOobject and a polyMesh reference
 467 Foam::globalMeshData::globalMeshData(const IOobject& io, const polyMesh& mesh)
 468 :
 469     processorTopology(mesh.boundaryMesh()),
 470     mesh_(mesh),
 471     bb_(mesh.points()),
 472     nTotalPoints_(-1),
 473     nTotalFaces_(-1),
 474     nTotalCells_(-1),
 475     processorPatches_(0),
 476     processorPatchIndices_(0),
 477     processorPatchNeighbours_(0),
 478     nGlobalPoints_(-1),
 479     sharedPointLabels_(0),
 480     sharedPointAddr_(0),
 481     sharedPointGlobalLabelsPtr_(NULL),
 482     nGlobalEdges_(-1),
 483     sharedEdgeLabelsPtr_(NULL),
 484     sharedEdgeAddrPtr_(NULL)
 485 {
 486     initProcAddr();
 487
 488     IOdictionary dict(io);
 489
 490     dict.lookup("nTotalPoints") >> nTotalPoints_;
 491     dict.lookup("nTotalFaces") >> nTotalFaces_;
 492     dict.lookup("nTotalCells") >> nTotalCells_;
 493     dict.lookup("nGlobalPoints") >> nGlobalPoints_;
 494     dict.lookup("sharedPointLabels") >> sharedPointLabels_;
 495     dict.lookup("sharedPointAddr") >> sharedPointAddr_;
 496     labelList sharedPointGlobalLabels(dict.lookup("sharedPointGlobalLabels"));
 497
 498     sharedPointGlobalLabelsPtr_ = new labelList(sharedPointGlobalLabels);
 499 }
 500
 501
 502 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 503
 504 Foam::globalMeshData::~globalMeshData()
 505 {
 506     clearOut();
 507 }
 508
 509
 510 void Foam::globalMeshData::clearOut()
 511 {
 512     deleteDemandDrivenData(sharedPointGlobalLabelsPtr_);
 513     // Edge
 514     nGlobalPoints_ = -1;
 515     deleteDemandDrivenData(sharedEdgeLabelsPtr_);
 516     deleteDemandDrivenData(sharedEdgeAddrPtr_);
 517 }
 518
 519
 520 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 521
 522 // Return shared point global labels.
 523 const Foam::labelList& Foam::globalMeshData::sharedPointGlobalLabels() const
 524 {
 525     if (!sharedPointGlobalLabelsPtr_)
 526     {
 527         sharedPointGlobalLabelsPtr_ = new labelList(sharedPointLabels_.size());
 528         labelList& sharedPointGlobalLabels = *sharedPointGlobalLabelsPtr_;
 529
 530         IOobject addrHeader
 531         (
 532             "pointProcAddressing",
 533             mesh_.facesInstance()/mesh_.meshSubDir,
 534             mesh_,
 535             IOobject::MUST_READ
 536         );
 537
 538         if (addrHeader.headerOk())
 539         {
 540             // There is a pointProcAddressing file so use it to get labels
 541             // on the original mesh
 542             Pout<< "globalMeshData::sharedPointGlobalLabels : "
 543                 << "Reading pointProcAddressing" << endl;
 544
 545             labelIOList pointProcAddressing(addrHeader);
 546
 547             forAll(sharedPointLabels_, i)
 548             {
 549                 // Get my mesh point
 550                 label pointI = sharedPointLabels_[i];
 551
 552                 // Map to mesh point of original mesh
 553                 sharedPointGlobalLabels[i] = pointProcAddressing[pointI];
 554             }
 555         }
 556         else
 557         {
 558             Pout<< "globalMeshData::sharedPointGlobalLabels :"
 559                 << " Setting pointProcAddressing to -1" << endl;
 560
 561             sharedPointGlobalLabels = -1;
 562         }
 563     }
 564     return *sharedPointGlobalLabelsPtr_;
 565 }
 566
 567
 568 // Collect coordinates of shared points. (does parallel communication!)
 569 Foam::pointField Foam::globalMeshData::sharedPoints() const
 570 {
 571     // Get all processors to send their shared points to master.
 572     // (not very efficient)
 573
 574     pointField sharedPoints(nGlobalPoints_);
 575
 576     if (Pstream::master())
 577     {
 578         // Master:
 579         // insert my own data first
 580         forAll(sharedPointLabels_, i)
 581         {
 582             label sharedPointI = sharedPointAddr_[i];
 583
 584             sharedPoints[sharedPointI] = mesh_.points()[sharedPointLabels_[i]];
 585         }
 586
 587         // Receive data from slaves and insert
 588         for
 589         (
 590             int slave=Pstream::firstSlave();
 591             slave<=Pstream::lastSlave();
 592             slave++
 593         )
 594         {
 595             IPstream fromSlave(Pstream::blocking, slave);
 596
 597             labelList nbrSharedPointAddr;
 598             pointField nbrSharedPoints;
 599             fromSlave >> nbrSharedPointAddr >> nbrSharedPoints;
 600
 601             forAll(nbrSharedPointAddr, i)
 602             {
 603                 label sharedPointI = nbrSharedPointAddr[i];
 604
 605                 sharedPoints[sharedPointI] = nbrSharedPoints[i];
 606             }
 607         }
 608
 609         // Send back
 610         for
 611         (
 612             int slave=Pstream::firstSlave();
 613             slave<=Pstream::lastSlave();
 614             slave++
 615         )
 616         {
 617             OPstream toSlave
 618             (
 619                 Pstream::blocking,
 620                 slave,
 621                 sharedPoints.size()*sizeof(vector::zero)
 622             );
 623             toSlave << sharedPoints;
 624         }
 625     }
 626     else
 627     {
 628         // Slave:
 629         // send points
 630         {
 631             OPstream toMaster(Pstream::blocking, Pstream::masterNo());
 632
 633             toMaster
 634                 << sharedPointAddr_
 635                 << UIndirectList<point>(mesh_.points(), sharedPointLabels_)();
 636         }
 637
 638         // Receive sharedPoints
 639         {
 640             IPstream fromMaster(Pstream::blocking, Pstream::masterNo());
 641             fromMaster >> sharedPoints;
 642         }
 643     }
 644
 645     return sharedPoints;
 646 }
 647
 648
 649 // Collect coordinates of shared points. (does parallel communication!)
 650 Foam::pointField Foam::globalMeshData::geometricSharedPoints() const
 651 {
 652     // Get coords of my shared points
 653     pointField sharedPoints(sharedPointLabels_.size());
 654
 655     forAll(sharedPointLabels_, i)
 656     {
 657         label meshPointI = sharedPointLabels_[i];
 658
 659         sharedPoints[i] = mesh_.points()[meshPointI];
 660     }
 661
 662     // Append from all processors
 663     combineReduce(sharedPoints, plusEqOp<pointField>());
 664
 665     // Merge tolerance
 666     scalar tolDim = matchTol_()*bb_.mag();
 667
 668     // And see how many are unique
 669     labelList pMap;
 670     pointField mergedPoints;
 671
 672     mergePoints
 673     (
 674         sharedPoints,   // coordinates to merge
 675         tolDim,         // tolerance
 676         false,          // verbosity
 677         pMap,
 678         mergedPoints
 679     );
 680
 681     return mergedPoints;
 682 }
 683
 684
 685 Foam::label Foam::globalMeshData::nGlobalEdges() const
 686 {
 687     if (nGlobalEdges_ == -1)
 688     {
 689         calcSharedEdges();
 690     }
 691     return nGlobalEdges_;
 692 }
 693
 694
 695 const Foam::labelList& Foam::globalMeshData::sharedEdgeLabels() const
 696 {
 697     if (!sharedEdgeLabelsPtr_)
 698     {
 699         calcSharedEdges();
 700     }
 701     return *sharedEdgeLabelsPtr_;
 702 }
 703
 704
 705 const Foam::labelList& Foam::globalMeshData::sharedEdgeAddr() const
 706 {
 707     if (!sharedEdgeAddrPtr_)
 708     {
 709         calcSharedEdges();
 710     }
 711     return *sharedEdgeAddrPtr_;
 712 }
 713
 714
 715 void Foam::globalMeshData::movePoints(const pointField& newPoints)
 716 {
 717     // Topology does not change and we don't store any geometry so nothing
 718     // needs to be done.
 719 }
 720
 721
 722 // Update all data after morph
 723 void Foam::globalMeshData::updateMesh()
 724 {
 725     // Clear out old data
 726     clearOut();
 727
 728     // Do processor patch addressing
 729     initProcAddr();
 730
 731     // Note: boundBox does reduce
 732     bb_ = boundBox(mesh_.points());
 733
 734     scalar tolDim = matchTol_()*bb_.mag();
 735
 736     if (debug)
 737     {
 738         Pout<< "globalMeshData : bb_:" << bb_
 739             << " merge dist:" << tolDim << endl;
 740     }
 741
 742
 743     // Option 1. Topological
 744     {
 745         // Calculate all shared points. This does all the hard work.
 746         globalPoints parallelPoints(mesh_);
 747
 748         // Copy data out.
 749         nGlobalPoints_ = parallelPoints.nGlobalPoints();
 750         sharedPointLabels_ = parallelPoints.sharedPointLabels();
 751         sharedPointAddr_ = parallelPoints.sharedPointAddr();
 752     }
 753
 754     //// Option 2. Geometric
 755     //{
 756     //    // Calculate all shared points. This does all the hard work.
 757     //    geomGlobalPoints parallelPoints(mesh_, tolDim);
 758     //
 759     //    // Copy data out.
 760     //    nGlobalPoints_ = parallelPoints.nGlobalPoints();
 761     //    sharedPointLabels_ = parallelPoints.sharedPointLabels();
 762     //    sharedPointAddr_ = parallelPoints.sharedPointAddr();
 763     //
 764     //    nGlobalEdges_ = parallelPoints.nGlobalEdges();
 765     //    sharedEdgeLabels_ = parallelPoints.sharedEdgeLabels();
 766     //    sharedEdgeAddr_ = parallelPoints.sharedEdgeAddr();
 767     //}
 768
 769     // Total number of faces. Start off from all faces. Remove coincident
 770     // processor faces (on highest numbered processor) before summing.
 771     nTotalFaces_ = mesh_.nFaces();
 772
 773     // Do not count processor-patch faces that are coincident.
 774     forAll(processorPatches_, i)
 775     {
 776         label patchI = processorPatches_[i];
 777
 778         const processorPolyPatch& procPatch =
 779             refCast<const processorPolyPatch>(mesh_.boundaryMesh()[patchI]);
 780
 781         if (Pstream::myProcNo() > procPatch.neighbProcNo())
 782         {
 783             // Uncount my faces. Handle cyclics separately.
 784
 785             if (procPatch.separated())
 786             {
 787                 const vectorField& separationDist = procPatch.separation();
 788
 789                 nTotalFaces_ -= countCoincidentFaces(tolDim, separationDist);
 790             }
 791             else
 792             {
 793                 // Normal, unseparated processor patch. Remove duplicates.
 794                 nTotalFaces_ -= procPatch.size();
 795             }
 796         }
 797     }
 798     reduce(nTotalFaces_, sumOp<label>());
 799
 800     if (debug)
 801     {
 802         Pout<< "globalMeshData : nTotalFaces_:" << nTotalFaces_ << endl;
 803     }
 804
 805
 806     nTotalCells_ = mesh_.nCells();
 807     reduce(nTotalCells_, sumOp<label>());
 808
 809     if (debug)
 810     {
 811         Pout<< "globalMeshData : nTotalCells_:" << nTotalCells_ << endl;
 812     }
 813
 814     nTotalPoints_ = mesh_.nPoints();
 815
 816     // Correct points for duplicate ones. We have
 817     // - points shared between 2 processor patches only. Count only on
 818     //   lower numbered processor. Make sure to count only once since points
 819     //   can be on multiple patches on the same processor.
 820     // - globally shared points.
 821
 822     if (Pstream::parRun())
 823     {
 824         const label UNSET = 0;      // not set
 825         const label SHARED = 1;     // globally shared
 826         const label VISITED = 2;    // corrected for
 827
 828         // Mark globally shared points
 829         PackedList<2> pointStatus(mesh_.nPoints(), UNSET);
 830
 831         forAll(sharedPointLabels_, i)
 832         {
 833             label meshPointI = sharedPointLabels_[i];
 834
 835             pointStatus.set(meshPointI, SHARED);
 836         }
 837
 838         // Send patch local points
 839         forAll(processorPatches_, i)
 840         {
 841             label patchI = processorPatches_[i];
 842
 843             const processorPolyPatch& procPatch =
 844                 refCast<const processorPolyPatch>
 845                 (
 846                     mesh_.boundaryMesh()[patchI]
 847                 );
 848
 849             OPstream toNeighbour(Pstream::blocking, procPatch.neighbProcNo());
 850
 851             toNeighbour << procPatch.localPoints();
 852         }
 853
 854         // Receive patch local points and uncount if coincident and not shared
 855         forAll(processorPatches_, i)
 856         {
 857             label patchI = processorPatches_[i];
 858
 859             const processorPolyPatch& procPatch =
 860                 refCast<const processorPolyPatch>
 861                 (
 862                     mesh_.boundaryMesh()[patchI]
 863                 );
 864
 865             IPstream fromNeighbour
 866             (
 867                 Pstream::blocking,
 868                 procPatch.neighbProcNo()
 869             );
 870
 871             pointField nbrPoints(fromNeighbour);
 872
 873             if (Pstream::myProcNo() > procPatch.neighbProcNo())
 874             {
 875                 labelList pMap;
 876                 matchPoints
 877                 (
 878                     procPatch.localPoints(),
 879                     nbrPoints,
 880                     scalarField(procPatch.nPoints(), tolDim),   // tolerance
 881                     false,      // verbosity
 882                     pMap        // map from my points to nbrPoints
 883                 );
 884
 885                 forAll(pMap, patchPointI)
 886                 {
 887                     label meshPointI = procPatch.meshPoints()[patchPointI];
 888
 889                     label stat = pointStatus.get(meshPointI);
 890
 891                     if (stat == UNSET)
 892                     {
 893                         // Mark point as visited so if point is on multipl
 894                         // e processor patches it only gets uncounted once
 895                         pointStatus.set(meshPointI, VISITED);
 896
 897                         if (pMap[patchPointI] != -1)
 898                         {
 899                             // Points share same coordinate so uncount.
 900                             nTotalPoints_--;
 901                         }
 902                     }
 903                 }
 904             }
 905         }
 906         // Sum all points
 907         reduce(nTotalPoints_, sumOp<label>());
 908     }
 909
 910     // nTotalPoints has not been corrected yet for shared points. For these
 911     // just collect all their coordinates and count unique ones.
 912
 913     label mySharedPoints = sharedPointLabels_.size();
 914     reduce(mySharedPoints, sumOp<label>());
 915
 916     // Collect and merge shared points (does parallel communication)
 917     pointField geomSharedPoints(geometricSharedPoints());
 918     label nGeomSharedPoints = geomSharedPoints.size();
 919
 920     // Shared points merged down to mergedPoints size.
 921     nTotalPoints_ -= mySharedPoints - nGeomSharedPoints;
 922
 923     if (debug)
 924     {
 925         Pout<< "globalMeshData : nTotalPoints_:" << nTotalPoints_ << endl;
 926     }
 927
 928     //
 929     // Now we have all info we wanted.
 930     // Do some checking (if debug is set)
 931     //
 932
 933     if (debug)
 934     {
 935         if (Pstream::master())
 936         {
 937             // We have the geometricSharedPoints already so write them.
 938             // Ideally would like to write the networks of connected points as
 939             // well but this is harder. (Todo)
 940             Pout<< "globalMeshData : writing geometrically separated shared"
 941                 << " points to geomSharedPoints.obj" << endl;
 942
 943             OFstream str("geomSharedPoints.obj");
 944
 945             forAll(geomSharedPoints, i)
 946             {
 947                 const point& pt = geomSharedPoints[i];
 948
 949                 str << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z()
 950                     << nl;
 951             }
 952         }
 953     }
 954 }
 955
 956
 957 // Write data
 958 bool Foam::globalMeshData::write() const
 959 {
 960     IOdictionary dict
 961     (
 962         IOobject
 963         (
 964             "parallelData",
 965             mesh_.facesInstance(),
 966             mesh_.meshSubDir,
 967             mesh_
 968         )
 969     );
 970
 971     dict.add("nTotalPoints", nTotalPoints());
 972     dict.add("nTotalFaces", nTotalFaces());
 973     dict.add("nTotalCells", nTotalCells());
 974
 975     dict.add("nGlobalPoints", nGlobalPoints());
 976     dict.add("sharedPointLabels", sharedPointLabels());
 977     dict.add("sharedPointAddr", sharedPointAddr());
 978     dict.add("sharedPointGlobalLabels", sharedPointGlobalLabels());
 979
 980     return dict.writeObject
 981     (
 982         IOstream::ASCII,
 983         IOstream::currentVersion,
 984         IOstream::UNCOMPRESSED
 985     );
 986 }
 987
 988
 989 // * * * * * * * * * * * * * * * Ostream Operators * * * * * * * * * * * * * //
 990
 991 Foam::Ostream& Foam::operator<<(Ostream& os, const globalMeshData& p)
 992 {
 993     os  << "nTotalPoints " << p.nTotalPoints() << token::END_STATEMENT << nl
 994         << "nTotalFaces " << p.nTotalFaces() << token::END_STATEMENT << nl
 995         << "nTotalCells " << p.nTotalCells() << token::END_STATEMENT << nl
 996         << "nGlobalPoints " << p.nGlobalPoints() << token::END_STATEMENT << nl
 997         << "sharedPointLabels " << p.sharedPointLabels()
 998         << token::END_STATEMENT << nl
 999         << "sharedPointAddr " << p.sharedPointAddr()
1000         << token::END_STATEMENT << endl;
1001
1002     return os;
1003 }
1004
1005
1006 // ************************************************************************* //