src/mesh/cfMesh/meshLibrary/utilities/octrees/meshOctree/meshOctreeAddressing/meshOctreeAddressingParallelAddressing.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
   6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh 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     cfMesh 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 cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "meshOctreeAddressing.H"
  29 #include "meshOctree.H"
  30 #include "demandDrivenData.H"
  31 #include "helperFunctions.H"
  32
  33 #include <map>
  34
  35 # ifdef USE_OMP
  36 #include <omp.h>
  37 # endif
  38
  39 //#define DEBUGAddressing
  40
  41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  42
  43 namespace Foam
  44 {
  45
  46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  47
  48 void meshOctreeAddressing::calcGlobalPointLabels() const
  49 {
  50     if( !Pstream::parRun() )
  51         FatalErrorIn("void meshOctreeAddressing::calcGlobalPointLabels() const")
  52             << "Cannot calculate global labels! Exitting" << exit(FatalError);
  53
  54     const VRWGraph& nodeLabels = this->nodeLabels();
  55     const FRWGraph<label, 8>& nodeLeaves = this->nodeLeaves();
  56
  57     //- allocate containers
  58     globalPointLabelPtr_ = new labelLongList(nodeLeaves.size(), -1);
  59     labelLongList& globalPointLabel = *globalPointLabelPtr_;
  60
  61     globalPointToLocalPtr_ = new Map<label>();
  62     Map<label>& globalToLocal = *globalPointToLocalPtr_;
  63
  64     pointProcsPtr_ = new VRWGraph(nodeLeaves.size());
  65     VRWGraph& pointProcs = *pointProcsPtr_;
  66
  67     //- find the number of points local to each processor
  68     //- The point is taken to be local if it belongs to one processor, only,
  69     //- or to the leaf with the smallest label
  70     labelList nLocalPoints(Pstream::nProcs(), 0);
  71
  72     forAll(nodeLeaves, pointI)
  73     {
  74         DynList<label> procs;
  75         procs.append(Pstream::myProcNo());
  76         FixedList<bool, 8> validLeaf(true);
  77
  78         for(label nlI=0;nlI<8;++nlI)
  79         {
  80             const label leafI = nodeLeaves(pointI, nlI);
  81
  82             if( leafI < 0 )
  83             {
  84                 validLeaf[nlI] = false;
  85                 continue;
  86             }
  87
  88             for(label i=0;i<nlI;++i)
  89                 if( nodeLeaves(pointI, i) == nodeLeaves(pointI, nlI) )
  90                 {
  91                     validLeaf[nlI] = false;
  92                     validLeaf[i] = false;
  93                 }
  94
  95             procs.appendIfNotIn(octree_.returnLeaf(leafI).procNo());
  96         }
  97
  98         label minLeaf(octree_.numberOfLeaves());
  99         bool found(false);
 100
 101         for(label nlI=0;nlI<8;++nlI)
 102         {
 103             if( !validLeaf[nlI] )
 104                 continue;
 105
 106             const label leafI = nodeLeaves(pointI, nlI);
 107
 108             minLeaf = Foam::min(leafI, minLeaf);
 109             found = true;
 110         }
 111
 112         if( found && octree_.returnLeaf(minLeaf).procNo() == Pstream::myProcNo() )
 113         {
 114             if( procs.size() > 1 )
 115                 pointProcs.setRow(pointI, procs);
 116
 117             globalPointLabel[pointI] = -2;
 118             ++nLocalPoints[Pstream::myProcNo()];
 119         }
 120     }
 121
 122     //- exchange data with other processors
 123     Pstream::gatherList(nLocalPoints);
 124     Pstream::scatterList(nLocalPoints);
 125
 126     //- find the starting point label
 127     label startPoint(0);
 128     for(label procI=0;procI<Pstream::myProcNo();++procI)
 129         startPoint += nLocalPoints[procI];
 130
 131     //- assign global labels to local points
 132     forAll(globalPointLabel, pointI)
 133     {
 134         if( globalPointLabel[pointI] == -2 )
 135         {
 136             globalPointLabel[pointI] = startPoint++;
 137
 138             if( pointProcs.sizeOfRow(pointI) != 0 )
 139                 globalToLocal.insert(globalPointLabel[pointI], pointI);
 140         }
 141     }
 142
 143     //- distribute the labels to other processors
 144     //- it is done by sending the global leaf label and the node labels
 145     //- to processors which contain the leaves as part of buffer layers
 146     //- it is performed in reduce-like manner
 147     const labelLongList& globalLeafLabel = this->globalLeafLabel();
 148     const Map<label>& globalToLocalLeaf = this->globalToLocalLeafAddressing();
 149     const VRWGraph& leafAtProcs = this->leafAtProcs();
 150     const labelList& neiProcs = octree_.neiProcs();
 151
 152     DynList<label> below, above;
 153     forAll(neiProcs, i)
 154     {
 155         if( neiProcs[i] < Pstream::myProcNo() )
 156         {
 157             above.append(neiProcs[i]);
 158         }
 159         else if( neiProcs[i] > Pstream::myProcNo() )
 160         {
 161             below.append(neiProcs[i]);
 162         }
 163     }
 164
 165     VRWGraph procLeaves;
 166     procLeaves.reverseAddressing(Pstream::nProcs(), leafAtProcs);
 167
 168     //- scatter the data from the processors above to the processors below
 169     //- receive the data from the processors above
 170     forAll(above, aboveI)
 171     {
 172         //- receive the data
 173         labelList receivedLabels;
 174         IPstream fromOtherProc(Pstream::blocking, above[aboveI]);
 175         fromOtherProc >> receivedLabels;
 176
 177         label counter(0);
 178         while( counter < receivedLabels.size() )
 179         {
 180             const label leafI = globalToLocalLeaf[receivedLabels[counter++]];
 181
 182             if( nodeLabels.sizeOfRow(leafI) == 0 )
 183                 FatalErrorIn
 184                 (
 185                    "void meshOctreeAddressing::"
 186                     "calcGlobalPointLabels() const"
 187                 ) << "1. Leaf " << leafI << " is not used in the mesh!"
 188                  << " Exitting.." << abort(FatalError);
 189
 190             for(label i=0;i<8;++i)
 191             {
 192                 const label nI = nodeLabels(leafI, i);
 193
 194                 const label globalLabel = receivedLabels[counter++];
 195
 196                 const label nProcs = receivedLabels[counter++];
 197                 for(label ppI=0;ppI<nProcs;++ppI)
 198                     pointProcs.appendIfNotIn(nI, receivedLabels[counter++]);
 199
 200                 if( globalLabel < 0 )
 201                     continue;
 202
 203                 label& gpl = globalPointLabel[nI];
 204
 205                 if( (gpl != -1) && (gpl != globalLabel) )
 206                     FatalErrorIn
 207                     (
 208                         "void meshOctreeAddressing::"
 209                         "calcGlobalPointLabels() const"
 210                     ) << "Wrong global label for point " << nI
 211                       << " Exitting.." << abort(FatalError);
 212
 213                 gpl = globalLabel;
 214                 globalToLocal.insert(globalLabel, nI);
 215             }
 216         }
 217     }
 218
 219     //- send the data to the processors below
 220     forAll(below, belowI)
 221     {
 222         const label neiProc = below[belowI];
 223
 224         labelLongList dts;
 225         forAllRow(procLeaves, neiProc, i)
 226         {
 227             const label leafI = procLeaves(neiProc, i);
 228
 229             if( nodeLabels.sizeOfRow(leafI) == 0 )
 230                 continue;
 231
 232             dts.append(globalLeafLabel[leafI]);
 233             for(label nI=0;nI<8;++nI)
 234             {
 235                 const label nodeI = nodeLabels(leafI, nI);
 236                 dts.append(globalPointLabel[nodeI]);
 237
 238                 //- add the current processor
 239                 pointProcs.appendIfNotIn(nodeI, Pstream::myProcNo());
 240
 241                 dts.append(pointProcs.sizeOfRow(nodeI));
 242                 forAllRow(pointProcs, nodeI, ppI)
 243                     dts.append(pointProcs(nodeI, ppI));
 244             }
 245         }
 246
 247         //- send the data
 248         OPstream toOtherProc(Pstream::blocking, neiProc, dts.byteSize());
 249         toOtherProc << dts;
 250     }
 251
 252     //- gather the data from the processors below to the processors above
 253     //- receive the data from the processors below
 254     forAllReverse(below, belowI)
 255     {
 256         //- receive the data
 257         labelList receivedLabels;
 258         IPstream fromOtherProc(Pstream::blocking, below[belowI]);
 259         fromOtherProc >> receivedLabels;
 260
 261         label counter(0);
 262         while( counter < receivedLabels.size() )
 263         {
 264             const label leafI = globalToLocalLeaf[receivedLabels[counter++]];
 265
 266             if( nodeLabels.sizeOfRow(leafI) == 0 )
 267                 FatalErrorIn
 268                 (
 269                    "void meshOctreeAddressing::"
 270                     "calcGlobalPointLabels() const"
 271                 ) << "2. Leaf " << leafI << " is not used in the mesh!"
 272                  << " Exitting.." << abort(FatalError);
 273
 274             for(label i=0;i<8;++i)
 275             {
 276                 const label nI = nodeLabels(leafI, i);
 277
 278                 const label globalLabel = receivedLabels[counter++];
 279
 280                 const label nProcs = receivedLabels[counter++];
 281                 for(label ppI=0;ppI<nProcs;++ppI)
 282                     pointProcs.appendIfNotIn(nI, receivedLabels[counter++]);
 283
 284                 if( globalLabel < 0 )
 285                     continue;
 286
 287                 label & gpl = globalPointLabel[nI];
 288
 289                 if( (gpl != -1) && (gpl != globalLabel) )
 290                     FatalErrorIn
 291                     (
 292                         "void meshOctreeAddressing::"
 293                         "calcGlobalPointLabels() const"
 294                     ) << "Wrong global label for point " << nI
 295                       << " Exitting.." << abort(FatalError);
 296
 297                 gpl = globalLabel;
 298                 globalToLocal.insert(globalLabel, nI);
 299             }
 300         }
 301     }
 302
 303     //- send the data to the processors below
 304     forAllReverse(above, aboveI)
 305     {
 306         const label neiProc = above[aboveI];
 307
 308         labelLongList dts;
 309         forAllRow(procLeaves, neiProc, i)
 310         {
 311             const label leafI = procLeaves(neiProc, i);
 312
 313             if( nodeLabels.sizeOfRow(leafI) == 0 )
 314                 continue;
 315
 316             dts.append(globalLeafLabel[leafI]);
 317             for(label nI=0;nI<8;++nI)
 318             {
 319                 const label nodeI = nodeLabels(leafI, nI);
 320                 dts.append(globalPointLabel[nodeI]);
 321
 322                 //- add the current processor
 323                 pointProcs.appendIfNotIn(nodeI, Pstream::myProcNo());
 324
 325                 dts.append(pointProcs.sizeOfRow(nodeI));
 326                 forAllRow(pointProcs, nodeI, ppI)
 327                     dts.append(pointProcs(nodeI, ppI));
 328             }
 329         }
 330
 331         //- send the data
 332         OPstream toOtherProc(Pstream::blocking, neiProc, dts.byteSize());
 333         toOtherProc << dts;
 334     }
 335 }
 336
 337 void meshOctreeAddressing::calcGlobalFaceLabels() const
 338 {
 339     if( !Pstream::parRun() )
 340         FatalErrorIn("void meshOctreeAddressing::calcGlobalFaceLabels() const")
 341             << "Cannot calculate global labels! Exitting" << exit(FatalError);
 342
 343     FatalError << "Not implemented" << exit(FatalError);
 344 }
 345
 346 void meshOctreeAddressing::calcGlobalLeafLabels() const
 347 {
 348     if( !Pstream::parRun() )
 349         FatalErrorIn("void meshOctreeAddressing::calcGlobalLeafLabels() const")
 350             << "Cannot calculate global labels! Exitting" << exit(FatalError);
 351
 352     //- allocate the memory
 353     globalLeafLabelPtr_ = new labelLongList(octree_.numberOfLeaves(), -1);
 354     labelLongList& globalLeafLabel = *globalLeafLabelPtr_;
 355
 356     globalLeafToLocalPtr_ = new Map<label>();
 357     Map<label>& globalToLocal = *globalLeafToLocalPtr_;
 358
 359     leafAtProcsPtr_ = new VRWGraph(octree_.numberOfLeaves());
 360     VRWGraph& leafAtProcs = *leafAtProcsPtr_;
 361
 362     //- find the number of leaves local to each processor
 363     labelList nLeavesAtProc(Pstream::nProcs(), 0);
 364
 365     label nLeaves(0);
 366
 367     # ifdef USE_OMP
 368     # pragma omp parallel for schedule(static) reduction(+:nLeaves)
 369     # endif
 370     for(label leafI=0;leafI<octree_.numberOfLeaves();++leafI)
 371     {
 372         const meshOctreeCubeBasic& oc = octree_.returnLeaf(leafI);
 373
 374         if( oc.procNo() == Pstream::myProcNo() )
 375             ++nLeaves;
 376     }
 377
 378     nLeavesAtProc[Pstream::myProcNo()] = nLeaves;
 379
 380     //- exchange the data with other processors
 381     Pstream::gatherList(nLeavesAtProc);
 382     Pstream::scatterList(nLeavesAtProc);
 383
 384     //- find the starting labels for
 385     nLeaves = 0;
 386     for(label procI=0;procI<Pstream::myProcNo();++procI)
 387         nLeaves += nLeavesAtProc[procI];
 388
 389     //- set the global labels to local leaves
 390     labelLongList otherProcLeaves;
 391     for(label leafI=0;leafI<octree_.numberOfLeaves();++leafI)
 392     {
 393         const meshOctreeCubeBasic& oc = octree_.returnLeaf(leafI);
 394
 395         if( oc.procNo() == Pstream::myProcNo() )
 396         {
 397             globalLeafLabel[leafI] = nLeaves++;
 398         }
 399         else
 400         {
 401             otherProcLeaves.append(leafI);
 402             leafAtProcs.append(leafI, Pstream::myProcNo());
 403             leafAtProcs.append(leafI, oc.procNo());
 404         }
 405     }
 406
 407     //- the rest of the code is needed in case an additional layer of
 408     //- of octree leaves belonging to other processors is added in order to
 409     //- simplify the procedure for generation of mesh templates
 410
 411     //- allocate the map for exchanging of data
 412     std::map<label, LongList<meshOctreeCubeBasic> > exchangeData;
 413     const labelList& neiProcs = octree_.neiProcs();
 414     forAll(neiProcs, procI)
 415     {
 416         const std::pair<label, LongList<meshOctreeCubeBasic> > pp
 417         (
 418             neiProcs[procI],
 419             LongList<meshOctreeCubeBasic>()
 420         );
 421
 422         exchangeData.insert(pp);
 423     }
 424
 425     //- Here we have to combine the information from all processors
 426     //- it is started such that all processors send the leaves to the processor
 427     //- that contains them locally
 428
 429     //- fill the map with data
 430     forAll(otherProcLeaves, i)
 431     {
 432         const meshOctreeCubeBasic& oc = octree_.returnLeaf(otherProcLeaves[i]);
 433         meshOctreeCubeBasic coc(oc);
 434         coc.setProcNo(Pstream::myProcNo());
 435         exchangeData[oc.procNo()].append(coc);
 436     }
 437
 438     //- exchange the data with other processors
 439     LongList<meshOctreeCubeBasic> rLeaves;
 440     help::exchangeMap(exchangeData, rLeaves, Pstream::scheduled);
 441
 442     //- update the local data
 443     forAll(rLeaves, i)
 444     {
 445         const label cLabel = octree_.findLeafLabelForPosition(rLeaves[i]);
 446
 447         globalToLocal.insert(globalLeafLabel[cLabel], cLabel);
 448         leafAtProcs.appendIfNotIn(cLabel, Pstream::myProcNo());
 449         leafAtProcs.appendIfNotIn(cLabel, rLeaves[i].procNo());
 450     }
 451
 452     //- now the global leaf labels shall be sent from the processors
 453     //- that own the leaves to the processors that also contain them
 454     std::map<label, labelLongList> exchangeLabels;
 455     std::map<label, LongList<meshOctreeCubeBasic> >::iterator it;
 456     for(it=exchangeData.begin();it!=exchangeData.end();++it)
 457     {
 458         it->second.clear();
 459         exchangeLabels.insert(std::make_pair(it->first, labelLongList()));
 460     }
 461
 462     //- fill in the data
 463     forAll(leafAtProcs, leafI)
 464     {
 465         if( octree_.returnLeaf(leafI).procNo() == Pstream::myProcNo() )
 466         {
 467             forAllRow(leafAtProcs, leafI, i)
 468             {
 469                 const label procI = leafAtProcs(leafI, i);
 470
 471                 if( procI == Pstream::myProcNo() )
 472                     continue;
 473
 474                 exchangeData[procI].append(octree_.returnLeaf(leafI));
 475                 exchangeLabels[procI].append(globalLeafLabel[leafI]);
 476             }
 477         }
 478     }
 479
 480     //- exchange the data
 481     rLeaves.clear();
 482     help::exchangeMap(exchangeData, rLeaves, Pstream::scheduled);
 483     labelLongList rLabels;
 484     help::exchangeMap(exchangeLabels, rLabels, Pstream::scheduled);
 485
 486     if( rLeaves.size() != rLabels.size() )
 487         FatalErrorIn("void meshOctreeAddressing::calcGlobalLeafLabels() const")
 488             << "Invalid list size!" << abort(FatalError);
 489
 490     //- set the labels to the leaves originating from other processors
 491     forAll(rLeaves, i)
 492     {
 493         const label cLabel = octree_.findLeafLabelForPosition(rLeaves[i]);
 494
 495         globalLeafLabel[cLabel] = rLabels[i];
 496         globalToLocal.insert(rLabels[i], cLabel);
 497     }
 498
 499     //- update leafAtProcs for all processors
 500     exchangeLabels.clear();
 501     forAll(neiProcs, procI)
 502         exchangeLabels.insert(std::make_pair(neiProcs[procI], labelLongList()));
 503
 504     forAllConstIter(Map<label>, globalToLocal, iter)
 505     {
 506         const label leafI = iter();
 507
 508         if( octree_.returnLeaf(leafI).procNo() != Pstream::myProcNo() )
 509             continue;
 510
 511         forAllRow(leafAtProcs, leafI, i)
 512         {
 513             const label procI = leafAtProcs(leafI, i);
 514
 515             if( procI == Pstream::myProcNo() )
 516                 continue;
 517
 518             labelLongList& dts = exchangeLabels[procI];
 519             dts.append(iter.key());
 520
 521             dts.append(leafAtProcs.sizeOfRow(leafI));
 522             forAllRow(leafAtProcs, leafI, j)
 523                 dts.append(leafAtProcs(leafI, j));
 524         }
 525     }
 526
 527     //- exchange the data
 528     rLabels.clear();
 529     help::exchangeMap(exchangeLabels, rLabels, Pstream::scheduled);
 530
 531     //- update the local data
 532     label counter(0);
 533     while( counter < rLabels.size() )
 534     {
 535         const label gLabel = rLabels[counter++];
 536
 537         if( !globalToLocal.found(gLabel) )
 538             FatalError << "Cannot find global label " << gLabel
 539                 << exit(FatalError);
 540
 541         const label leafI = globalToLocal[gLabel];
 542
 543         const label numberOfProcs = rLabels[counter++];
 544         for(label i=0;i<numberOfProcs;++i)
 545             leafAtProcs.append(leafI, rLabels[counter++]);
 546     }
 547
 548     # ifdef DEBUGAddressing
 549     returnReduce(1, sumOp<label>());
 550     const List<direction>& boxType = this->boxType();
 551     const VRWGraph& nodeLabels = this->nodeLabels();
 552     for(label i=0;i<Pstream::nProcs();++i)
 553     {
 554         if( i == Pstream::myProcNo() )
 555         {
 556             forAll(globalLeafLabel, leafI)
 557             {
 558                 Pout << "Leaf " << leafI << "of type " << label(boxType[leafI])
 559                     << " has global label " << globalLeafLabel[leafI]
 560                     << " and coordinates " << octree_.returnLeaf(leafI)
 561                     << " and located at procs " << leafAtProcs[leafI]
 562                     << " node labels " << nodeLabels[leafI] << endl;
 563
 564                 if( octree_.returnLeaf(leafI).procNo() == Pstream::myProcNo() )
 565                     continue;
 566                 if( globalToLocal[globalLeafLabel[leafI]] != leafI )
 567                     FatalError << "Crap!!" << abort(FatalError);
 568             }
 569         }
 570
 571         returnReduce(1, sumOp<label>());
 572     }
 573
 574     //- check if the leaf at procs is ok
 575     exchangeData.clear();
 576     for(label i=0;i<Pstream::nProcs();++i)
 577     {
 578         if( i == Pstream::myProcNo() )
 579             continue;
 580
 581         exchangeData[i];
 582     }
 583
 584     for(label leafI=0;leafI<octree_.numberOfLeaves();++leafI)
 585         for(label i=0;i<Pstream::nProcs();++i)
 586         {
 587             if( i == Pstream::myProcNo() )
 588                 continue;
 589
 590             exchangeData[i].append(octree_.returnLeaf(leafI));
 591         }
 592
 593     std::map<label, List<meshOctreeCubeBasic> > rMap;
 594     help::exchangeMap(exchangeData, rMap);
 595
 596     for(std::map<label, List<meshOctreeCubeBasic> >::const_iterator it=rMap.begin();it!=rMap.end();++it)
 597     {
 598         const List<meshOctreeCubeBasic>& data = it->second;
 599
 600         forAll(data, i)
 601         {
 602             const label leafI = octree_.findLeafLabelForPosition(data[i]);
 603
 604             if( leafI < 0 )
 605                 continue;
 606
 607             if( !leafAtProcs.contains(leafI, it->first) )
 608                 FatalError << "Problem!!" << leafI
 609                     << " does not contain processor " << it->first
 610                     << " contains procs " << leafAtProcs[leafI]
 611                     << abort(FatalError);
 612         }
 613     }
 614
 615     returnReduce(1, sumOp<label>());
 616     # endif
 617 }
 618
 619 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 620
 621 } // End namespace Foam
 622
 623 // ************************************************************************* //