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) 1991-2008 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "globalMeshData.H"
  28 #include "Time.H"
  29 #include "Pstream.H"
  30 #include "PstreamCombineReduceOps.H"
  31 #include "processorPolyPatch.H"
  32 #include "demandDrivenData.H"
  33 #include "globalPoints.H"
  34 //#include "geomGlobalPoints.H"
  35 #include "labelIOList.H"
  36 #include "PackedList.H"
  37 #include "mergePoints.H"
  38 #include "matchPoints.H"
  39 #include "OFstream.H"
  40
  41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  42
  43 namespace Foam
  44 {
  45     defineTypeNameAndDebug(globalMeshData, 0);
  46 }
  47
  48 // Geometric matching tolerance. Factor of mesh bounding box.
  49 const Foam::scalar Foam::globalMeshData::matchTol_ = 1E-8;
  50
  51
  52 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  53
  54 // Collect processor patch addressing.
  55 void Foam::globalMeshData::initProcAddr()
  56 {
  57     processorPatchIndices_.setSize(mesh_.boundaryMesh().size());
  58     processorPatchIndices_ = -1;
  59
  60     processorPatchNeighbours_.setSize(mesh_.boundaryMesh().size());
  61     processorPatchNeighbours_ = -1;
  62
  63     // Construct processor patch indexing. processorPatchNeighbours_ only
  64     // set if running in parallel!
  65     processorPatches_.setSize(mesh_.boundaryMesh().size());
  66
  67     label nNeighbours = 0;
  68
  69     forAll (mesh_.boundaryMesh(), patchi)
  70     {
  71         if (isA<processorPolyPatch>(mesh_.boundaryMesh()[patchi]))
  72         {
  73             processorPatches_[nNeighbours] = patchi;
  74             processorPatchIndices_[patchi] = nNeighbours++;
  75         }
  76     }
  77     processorPatches_.setSize(nNeighbours);
  78
  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 indices
 217                 // 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     dynSharedEdgeLabels.shrink();
 397     sharedEdgeLabelsPtr_ = new labelList();
 398     labelList& sharedEdgeLabels = *sharedEdgeLabelsPtr_;
 399     sharedEdgeLabels.transfer(dynSharedEdgeLabels);
 400     dynSharedEdgeLabels.clear();
 401
 402     dynSharedEdgeAddr.shrink();
 403     sharedEdgeAddrPtr_ = new labelList();
 404     labelList& sharedEdgeAddr = *sharedEdgeAddrPtr_;
 405     sharedEdgeAddr.transfer(dynSharedEdgeAddr);
 406     dynSharedEdgeAddr.clear();
 407
 408     if (debug)
 409     {
 410         Pout<< "globalMeshData : nGlobalEdges_:" << nGlobalEdges_ << nl
 411             << "globalMeshData : sharedEdgeLabels:" << sharedEdgeLabels.size()
 412             << nl
 413             << "globalMeshData : sharedEdgeAddr:" << sharedEdgeAddr.size()
 414             << endl;
 415     }
 416 }
 417
 418
 419 // Helper function to count coincident faces. This part used to be
 420 // in updateMesh but I've had to move it to a separate function
 421 // because of aliasing optimisation errors in icc9.1 on the
 422 // Itanium.
 423 Foam::label Foam::globalMeshData::countCoincidentFaces
 424 (
 425     const scalar tolDim,
 426     const vectorField& separationDist
 427 )
 428 {
 429     label nCoincident = 0;
 430
 431     forAll(separationDist, faceI)
 432     {
 433         if (mag(separationDist[faceI]) < tolDim)
 434         {
 435             // Faces coincide
 436             nCoincident++;
 437         }
 438     }
 439     return nCoincident;
 440 }
 441
 442
 443 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 444
 445 // Construct from polyMesh
 446 Foam::globalMeshData::globalMeshData(const polyMesh& mesh)
 447 :
 448     processorTopology(mesh.boundaryMesh()),
 449     mesh_(mesh),
 450     bb_(vector::zero, vector::zero),
 451     nTotalPoints_(-1),
 452     nTotalFaces_(-1),
 453     nTotalCells_(-1),
 454     processorPatches_(0),
 455     processorPatchIndices_(0),
 456     processorPatchNeighbours_(0),
 457     nGlobalPoints_(-1),
 458     sharedPointLabels_(0),
 459     sharedPointAddr_(0),
 460     sharedPointGlobalLabelsPtr_(NULL),
 461     nGlobalEdges_(-1),
 462     sharedEdgeLabelsPtr_(NULL),
 463     sharedEdgeAddrPtr_(NULL)
 464 {
 465     updateMesh();
 466 }
 467
 468
 469 // Read constructor given IOobject and a polyMesh reference
 470 Foam::globalMeshData::globalMeshData(const IOobject& io, const polyMesh& mesh)
 471 :
 472     processorTopology(mesh.boundaryMesh()),
 473     mesh_(mesh),
 474     bb_(mesh.points()),
 475     nTotalPoints_(-1),
 476     nTotalFaces_(-1),
 477     nTotalCells_(-1),
 478     processorPatches_(0),
 479     processorPatchIndices_(0),
 480     processorPatchNeighbours_(0),
 481     nGlobalPoints_(-1),
 482     sharedPointLabels_(0),
 483     sharedPointAddr_(0),
 484     sharedPointGlobalLabelsPtr_(NULL),
 485     nGlobalEdges_(-1),
 486     sharedEdgeLabelsPtr_(NULL),
 487     sharedEdgeAddrPtr_(NULL)
 488 {
 489     initProcAddr();
 490
 491     IOdictionary dict(io);
 492
 493     dict.lookup("nTotalPoints") >> nTotalPoints_;
 494     dict.lookup("nTotalFaces") >> nTotalFaces_;
 495     dict.lookup("nTotalCells") >> nTotalCells_;
 496     dict.lookup("nGlobalPoints") >> nGlobalPoints_;
 497     dict.lookup("sharedPointLabels") >> sharedPointLabels_;
 498     dict.lookup("sharedPointAddr") >> sharedPointAddr_;
 499     labelList sharedPointGlobalLabels(dict.lookup("sharedPointGlobalLabels"));
 500
 501     sharedPointGlobalLabelsPtr_ = new labelList(sharedPointGlobalLabels);
 502 }
 503
 504
 505 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 506
 507 Foam::globalMeshData::~globalMeshData()
 508 {
 509     clearOut();
 510 }
 511
 512
 513 void Foam::globalMeshData::clearOut()
 514 {
 515     deleteDemandDrivenData(sharedPointGlobalLabelsPtr_);
 516     // Edge
 517     nGlobalPoints_ = -1;
 518     deleteDemandDrivenData(sharedEdgeLabelsPtr_);
 519     deleteDemandDrivenData(sharedEdgeAddrPtr_);
 520 }
 521
 522
 523 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 524
 525 // Return shared point global labels.
 526 const Foam::labelList& Foam::globalMeshData::sharedPointGlobalLabels() const
 527 {
 528     if (!sharedPointGlobalLabelsPtr_)
 529     {
 530         sharedPointGlobalLabelsPtr_ = new labelList(sharedPointLabels_.size());
 531         labelList& sharedPointGlobalLabels = *sharedPointGlobalLabelsPtr_;
 532
 533         IOobject addrHeader
 534         (
 535             "pointProcAddressing",
 536             mesh_.facesInstance()/mesh_.meshSubDir,
 537             mesh_,
 538             IOobject::MUST_READ
 539         );
 540
 541         if (addrHeader.headerOk())
 542         {
 543             // There is a pointProcAddressing file so use it to get labels
 544             // on the original mesh
 545             Pout<< "globalMeshData::sharedPointGlobalLabels : "
 546                 << "Reading pointProcAddressing" << endl;
 547
 548             labelIOList pointProcAddressing(addrHeader);
 549
 550             forAll(sharedPointLabels_, i)
 551             {
 552                 // Get my mesh point
 553                 label pointI = sharedPointLabels_[i];
 554
 555                 // Map to mesh point of original mesh
 556                 sharedPointGlobalLabels[i] = pointProcAddressing[pointI];
 557             }
 558         }
 559         else
 560         {
 561             Pout<< "globalMeshData::sharedPointGlobalLabels :"
 562                 << " Setting pointProcAddressing to -1" << endl;
 563
 564             sharedPointGlobalLabels = -1;
 565         }
 566     }
 567     return *sharedPointGlobalLabelsPtr_;
 568 }
 569
 570
 571 // Collect coordinates of shared points. (does parallel communication!)
 572 Foam::pointField Foam::globalMeshData::sharedPoints() const
 573 {
 574     // Get all processors to send their shared points to master.
 575     // (not very efficient)
 576
 577     pointField sharedPoints(nGlobalPoints_);
 578
 579     if (Pstream::master())
 580     {
 581         // Master:
 582         // insert my own data first
 583         forAll(sharedPointLabels_, i)
 584         {
 585             label sharedPointI = sharedPointAddr_[i];
 586
 587             sharedPoints[sharedPointI] = mesh_.points()[sharedPointLabels_[i]];
 588         }
 589
 590         // Receive data from slaves and insert
 591         for
 592         (
 593             int slave=Pstream::firstSlave();
 594             slave<=Pstream::lastSlave();
 595             slave++
 596         )
 597         {
 598             IPstream fromSlave(Pstream::blocking, slave);
 599
 600             labelList nbrSharedPointAddr;
 601             pointField nbrSharedPoints;
 602             fromSlave >> nbrSharedPointAddr >> nbrSharedPoints;
 603
 604             forAll(nbrSharedPointAddr, i)
 605             {
 606                 label sharedPointI = nbrSharedPointAddr[i];
 607
 608                 sharedPoints[sharedPointI] = nbrSharedPoints[i];
 609             }
 610         }
 611
 612         // Send back
 613         for
 614         (
 615             int slave=Pstream::firstSlave();
 616             slave<=Pstream::lastSlave();
 617             slave++
 618         )
 619         {
 620             OPstream toSlave
 621             (
 622                 Pstream::blocking,
 623                 slave,
 624                 sharedPoints.size()*sizeof(vector::zero)
 625             );
 626             toSlave << sharedPoints;
 627         }
 628     }
 629     else
 630     {
 631         // Slave:
 632         // send points
 633         {
 634             OPstream toMaster(Pstream::blocking, Pstream::masterNo());
 635
 636             toMaster
 637                 << sharedPointAddr_
 638                 << IndirectList<point>(mesh_.points(), sharedPointLabels_)();
 639         }
 640
 641         // Receive sharedPoints
 642         {
 643             IPstream fromMaster(Pstream::blocking, Pstream::masterNo());
 644             fromMaster >> sharedPoints;
 645         }
 646     }
 647
 648     return sharedPoints;
 649 }
 650
 651
 652 // Collect coordinates of shared points. (does parallel communication!)
 653 Foam::pointField Foam::globalMeshData::geometricSharedPoints() const
 654 {
 655     // Get coords of my shared points
 656     pointField sharedPoints(sharedPointLabels_.size());
 657
 658     forAll(sharedPointLabels_, i)
 659     {
 660         label meshPointI = sharedPointLabels_[i];
 661
 662         sharedPoints[i] = mesh_.points()[meshPointI];
 663     }
 664
 665     // Append from all processors
 666     combineReduce(sharedPoints, plusEqOp<pointField>());
 667
 668     // Merge tolerance
 669     scalar tolDim = matchTol_*mag(bb_.max() - bb_.min());
 670
 671     // And see how many are unique
 672     labelList pMap;
 673     pointField mergedPoints;
 674
 675     mergePoints
 676     (
 677         sharedPoints,   // coordinates to merge
 678         tolDim,         // tolerance
 679         false,          // verbosity
 680         pMap,
 681         mergedPoints
 682     );
 683
 684     return mergedPoints;
 685 }
 686
 687
 688 Foam::label Foam::globalMeshData::nGlobalEdges() const
 689 {
 690     if (nGlobalEdges_ == -1)
 691     {
 692         calcSharedEdges();
 693     }
 694     return nGlobalEdges_;
 695 }
 696
 697
 698 const Foam::labelList& Foam::globalMeshData::sharedEdgeLabels() const
 699 {
 700     if (!sharedEdgeLabelsPtr_)
 701     {
 702         calcSharedEdges();
 703     }
 704     return *sharedEdgeLabelsPtr_;
 705 }
 706
 707
 708 const Foam::labelList& Foam::globalMeshData::sharedEdgeAddr() const
 709 {
 710     if (!sharedEdgeAddrPtr_)
 711     {
 712         calcSharedEdges();
 713     }
 714     return *sharedEdgeAddrPtr_;
 715 }
 716
 717
 718 void Foam::globalMeshData::movePoints(const pointField& newPoints)
 719 {
 720     // Topology does not change and we don't store any geometry so nothing
 721     // needs to be done.
 722 }
 723
 724
 725 // Update all data after morph
 726 void Foam::globalMeshData::updateMesh()
 727 {
 728     // Clear out old data
 729     clearOut();
 730
 731     // Do processor patch addressing
 732     initProcAddr();
 733
 734     // Bounding box (does communication)
 735     bb_ = boundBox(mesh_.points(), true);
 736
 737     scalar tolDim = matchTol_*mag(bb_.max() - bb_.min());
 738
 739     if (debug)
 740     {
 741         Pout<< "globalMeshData : bb_:" << bb_
 742             << " merge dist:" << tolDim << endl;
 743     }
 744
 745
 746
 747     // Option 1. Topological
 748     {
 749         // Calculate all shared points. This does all the hard work.
 750         globalPoints parallelPoints(mesh_);
 751
 752         // Copy data out.
 753         nGlobalPoints_ = parallelPoints.nGlobalPoints();
 754         sharedPointLabels_ = parallelPoints.sharedPointLabels();
 755         sharedPointAddr_ = parallelPoints.sharedPointAddr();
 756     }
 757     //// Option 2. Geometric
 758     //{
 759     //    // Calculate all shared points. This does all the hard work.
 760     //    geomGlobalPoints parallelPoints(mesh_, tolDim);
 761     //
 762     //    // Copy data out.
 763     //    nGlobalPoints_ = parallelPoints.nGlobalPoints();
 764     //    sharedPointLabels_ = parallelPoints.sharedPointLabels();
 765     //    sharedPointAddr_ = parallelPoints.sharedPointAddr();
 766     //
 767     //    nGlobalEdges_ = parallelPoints.nGlobalEdges();
 768     //    sharedEdgeLabels_ = parallelPoints.sharedEdgeLabels();
 769     //    sharedEdgeAddr_ = parallelPoints.sharedEdgeAddr();
 770     //}
 771
 772     // Total number of faces. Start off from all faces. Remove coincident
 773     // processor faces (on highest numbered processor) before summing.
 774     nTotalFaces_ = mesh_.nFaces();
 775
 776     // Do not count processorpatch faces that are coincident.
 777     forAll(processorPatches_, i)
 778     {
 779         label patchI = processorPatches_[i];
 780
 781         const processorPolyPatch& procPatch =
 782             refCast<const processorPolyPatch>(mesh_.boundaryMesh()[patchI]);
 783
 784         if (Pstream::myProcNo() > procPatch.neighbProcNo())
 785         {
 786             // Uncount my faces. Handle cyclics separately.
 787
 788             if (procPatch.separated())
 789             {
 790                 const vectorField& separationDist = procPatch.separation();
 791
 792                 nTotalFaces_ -= countCoincidentFaces(tolDim, separationDist);
 793             }
 794             else
 795             {
 796                 // Normal, unseparated processor patch. Remove duplicates.
 797                 nTotalFaces_ -= procPatch.size();
 798             }
 799         }
 800     }
 801     reduce(nTotalFaces_, sumOp<label>());
 802
 803     if (debug)
 804     {
 805         Pout<< "globalMeshData : nTotalFaces_:" << nTotalFaces_ << endl;
 806     }
 807
 808
 809     nTotalCells_ = mesh_.nCells();
 810     reduce(nTotalCells_, sumOp<label>());
 811
 812     if (debug)
 813     {
 814         Pout<< "globalMeshData : nTotalCells_:" << nTotalCells_ << endl;
 815     }
 816
 817     nTotalPoints_ = mesh_.nPoints();
 818
 819     // Correct points for duplicate ones. We have
 820     // - points shared between 2 processor patches only. Count only on
 821     //   lower numbered processor. Make sure to count only once since points
 822     //   can be on multiple patches on the same processor.
 823     // - globally shared points.
 824
 825     if (Pstream::parRun())
 826     {
 827         const label UNSET = 0;      // not set
 828         const label SHARED = 1;     // globally shared
 829         const label VISITED = 2;    // corrected for
 830
 831         // Mark globally shared points
 832         PackedList<2> pointStatus(mesh_.nPoints(), UNSET);
 833
 834         forAll(sharedPointLabels_, i)
 835         {
 836             label meshPointI = sharedPointLabels_[i];
 837
 838             pointStatus.set(meshPointI, SHARED);
 839         }
 840
 841         // Send patch local points
 842         forAll(processorPatches_, i)
 843         {
 844             label patchI = processorPatches_[i];
 845
 846             const processorPolyPatch& procPatch =
 847                 refCast<const processorPolyPatch>(mesh_.boundaryMesh()[patchI]);
 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>(mesh_.boundaryMesh()[patchI]);
 861
 862             IPstream fromNeighbour(Pstream::blocking, procPatch.neighbProcNo());
 863
 864             pointField nbrPoints(fromNeighbour);
 865
 866             if (Pstream::myProcNo() > procPatch.neighbProcNo())
 867             {
 868                 labelList pMap;
 869                 matchPoints
 870                 (
 871                     procPatch.localPoints(),
 872                     nbrPoints,
 873                     scalarField(procPatch.nPoints(), tolDim),   // tolerance
 874                     false,      // verbosity
 875                     pMap        // map from my points to nbrPoints
 876                 );
 877
 878                 forAll(pMap, patchPointI)
 879                 {
 880                     label meshPointI = procPatch.meshPoints()[patchPointI];
 881
 882                     label stat = pointStatus.get(meshPointI);
 883
 884                     if (stat == UNSET)
 885                     {
 886                         // Mark point as visited so if point is on multiple proc
 887                         // patches it only gets uncounted once.
 888                         pointStatus.set(meshPointI, VISITED);
 889
 890                         if (pMap[patchPointI] != -1)
 891                         {
 892                             // Points share same coordinate so uncount.
 893                             nTotalPoints_--;
 894                         }
 895                     }
 896                 }
 897             }
 898         }
 899         // Sum all points
 900         reduce(nTotalPoints_, sumOp<label>());
 901     }
 902
 903     // nTotalPoints has not been corrected yet for shared points. For these
 904     // just collect all their coordinates and count unique ones.
 905
 906     label mySharedPoints = sharedPointLabels_.size();
 907     reduce(mySharedPoints, sumOp<label>());
 908
 909     // Collect and merge shared points (does parallel communication)
 910     pointField geomSharedPoints(geometricSharedPoints());
 911     label nGeomSharedPoints = geomSharedPoints.size();
 912
 913     // Shared points merged down to mergedPoints size.
 914     nTotalPoints_ -= mySharedPoints - nGeomSharedPoints;
 915
 916     if (debug)
 917     {
 918         Pout<< "globalMeshData : nTotalPoints_:" << nTotalPoints_ << endl;
 919     }
 920
 921     //
 922     // Now we have all info we wanted.
 923     // Do some checking (if debug is set)
 924     //
 925
 926     if (debug)
 927     {
 928         if (Pstream::master())
 929         {
 930             // We have the geometricSharedPoints already so write them.
 931             // Ideally would like to write the networks of connected points as
 932             // well but this is harder. (Todo)
 933             Pout<< "globalMeshData : writing geometrically separated shared"
 934                 << " points to geomSharedPoints.obj" << endl;
 935
 936             OFstream str("geomSharedPoints.obj");
 937
 938             forAll(geomSharedPoints, i)
 939             {
 940                 const point& pt = geomSharedPoints[i];
 941
 942                 str << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z()
 943                     << nl;
 944             }
 945         }
 946     }
 947 }
 948
 949
 950 // Write data
 951 bool Foam::globalMeshData::write() const
 952 {
 953     IOdictionary dict
 954     (
 955         IOobject
 956         (
 957             "parallelData",
 958             mesh_.facesInstance(),
 959             mesh_.meshSubDir,
 960             mesh_
 961         )
 962     );
 963
 964     dict.add("nTotalPoints", nTotalPoints());
 965     dict.add("nTotalFaces", nTotalFaces());
 966     dict.add("nTotalCells", nTotalCells());
 967
 968     dict.add("nGlobalPoints", nGlobalPoints());
 969     dict.add("sharedPointLabels", sharedPointLabels());
 970     dict.add("sharedPointAddr", sharedPointAddr());
 971     dict.add("sharedPointGlobalLabels", sharedPointGlobalLabels());
 972
 973     return dict.writeObject
 974     (
 975         IOstream::ASCII,
 976         IOstream::currentVersion,
 977         IOstream::UNCOMPRESSED
 978     );
 979 }
 980
 981
 982 // * * * * * * * * * * * * * * * Ostream Operators * * * * * * * * * * * * * //
 983
 984 Foam::Ostream& Foam::operator<<(Ostream& os, const globalMeshData& p)
 985 {
 986     os  << "nTotalPoints " << p.nTotalPoints() << token::END_STATEMENT << nl
 987         << "nTotalFaces " << p.nTotalFaces() << token::END_STATEMENT << nl
 988         << "nTotalCells " << p.nTotalCells() << token::END_STATEMENT << nl
 989         << "nGlobalPoints " << p.nGlobalPoints() << token::END_STATEMENT << nl
 990         << "sharedPointLabels " << p.sharedPointLabels()
 991         << token::END_STATEMENT << nl
 992         << "sharedPointAddr " << p.sharedPointAddr()
 993         << token::END_STATEMENT << endl;
 994
 995     return os;
 996 }
 997
 998
 999 // ************************************************************************* //