src/mesh/cfMesh/utilities/surfaceTools/meshSurfaceEngine/meshSurfaceEngineCalculateBoundaryNodesAndFaces.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 "meshSurfaceEngine.H"
  29 #include "demandDrivenData.H"
  30 #include "boolList.H"
  31 #include "helperFunctions.H"
  32 #include "helperFunctionsPar.H"
  33 #include "VRWGraphSMPModifier.H"
  34 #include "labelledPoint.H"
  35 #include "HashSet.H"
  36
  37 #include <map>
  38 #include <set>
  39
  40 # ifdef USE_OMP
  41 #include <omp.h>
  42 # endif
  43
  44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  45
  46 namespace Foam
  47 {
  48
  49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  50
  51 void meshSurfaceEngine::calculateBoundaryFaces() const
  52 {
  53     if( mesh_.boundaries().size() != 0 )
  54     {
  55         const faceListPMG& faces = mesh_.faces();
  56         const PtrList<boundaryPatch>& boundaries = mesh_.boundaries();
  57
  58         label nBoundaryFaces(0);
  59         if( activePatch_ < 0 )
  60         {
  61             //- take all patches
  62             forAll(boundaries, patchI)
  63                 nBoundaryFaces += boundaries[patchI].patchSize();
  64
  65             boundaryFacesPtr_ =
  66                 new faceList::subList
  67                 (
  68                     faces,
  69                     nBoundaryFaces,
  70                     boundaries[0].patchStart()
  71                 );
  72         }
  73         else if( activePatch_ < boundaries.size() )
  74         {
  75             nBoundaryFaces = boundaries[activePatch_].patchSize();
  76
  77             boundaryFacesPtr_ =
  78                 new faceList::subList
  79                 (
  80                     faces,
  81                     nBoundaryFaces,
  82                     boundaries[activePatch_].patchStart()
  83                 );
  84         }
  85         else
  86         {
  87             FatalErrorIn
  88             (
  89                 "void meshSurfaceEngine::calculateBoundaryFaces() const"
  90             ) << "Cannot select boundary faces. Invalid patch index "
  91               << activePatch_ << exit(FatalError);
  92         }
  93
  94         reduce(nBoundaryFaces, sumOp<label>());
  95         Info << "Found " << nBoundaryFaces << " boundary faces " << endl;
  96     }
  97     else
  98     {
  99         FatalErrorIn
 100         (
 101             "void meshSurfaceEngine::calculateBoundaryFaces() const"
 102         ) << "Boundary faces are not at the end of the face list!"
 103             << exit(FatalError);
 104     }
 105 }
 106
 107 void meshSurfaceEngine::calculateBoundaryOwners() const
 108 {
 109     const labelList& owner = mesh_.owner();
 110
 111     const faceList::subList& boundaryFaces = this->boundaryFaces();
 112
 113     if( !boundaryFaceOwnersPtr_ )
 114     boundaryFaceOwnersPtr_ = new labelList(boundaryFaces.size());
 115
 116     labelList& owners = *boundaryFaceOwnersPtr_;
 117
 118     const label start = mesh_.boundaries()[0].patchStart();
 119
 120     # ifdef USE_OMP
 121     # pragma omp parallel for schedule(static, 1)
 122     # endif
 123     forAll(boundaryFaces, fI)
 124     owners[fI] = owner[start+fI];
 125 }
 126
 127 void meshSurfaceEngine::calculateBoundaryNodes() const
 128 {
 129     //- mark boundary points
 130     label pointI(0);
 131     if( !bppPtr_ )
 132         bppPtr_ = new labelList(mesh_.points().size(), -1);
 133     labelList& bp = *bppPtr_;
 134
 135     const faceList::subList& boundaryFaces = this->boundaryFaces();
 136
 137     boolList isBndPoint(bp.size(), false);
 138
 139     # ifdef USE_OMP
 140     const label nThreads = 3 * omp_get_num_procs();
 141     # pragma omp parallel for num_threads(nThreads) schedule(static, 1)
 142     # endif
 143     forAll(boundaryFaces, bfI)
 144     {
 145         const face& bf = boundaryFaces[bfI];
 146
 147         forAll(bf, pI)
 148             isBndPoint[bf[pI]] = true;
 149     }
 150
 151     forAll(isBndPoint, pI)
 152     {
 153         if( isBndPoint[pI] )
 154             bp[pI] = pointI++;
 155     }
 156
 157     if( Pstream::parRun() )
 158     {
 159         const faceListPMG& faces = mesh_.faces();
 160         const PtrList<processorBoundaryPatch>& procBoundaries =
 161             mesh_.procBoundaries();
 162
 163         //- exchange information with processors
 164         //- this is need sometimes to find all nodes at the boundary
 165         bool found;
 166         do
 167         {
 168             found = false;
 169
 170             //- send bnd nodes to other processor
 171             forAll(procBoundaries, patchI)
 172             {
 173                 const label start = procBoundaries[patchI].patchStart();
 174                 const label end = start + procBoundaries[patchI].patchSize();
 175
 176                 labelLongList dts;
 177                 labelHashSet addedPoint;
 178                 for(label faceI=start;faceI<end;++faceI)
 179                 {
 180                     const face& f = faces[faceI];
 181                     forAll(f, pI)
 182                         if( (bp[f[pI]] != -1) && !addedPoint.found(f[pI]) )
 183                         {
 184                             addedPoint.insert(f[pI]);
 185
 186                             //- data is sent as follows
 187                             //- 1. local face label in patch
 188                             //- 2. local node in face
 189                             dts.append(faceI-start);
 190                             dts.append((f.size() - pI) % f.size());
 191                         }
 192                 }
 193
 194                 OPstream toOtherProc
 195                 (
 196                     Pstream::blocking,
 197                     procBoundaries[patchI].neiProcNo(),
 198                     dts.byteSize()
 199                 );
 200                 toOtherProc << dts;
 201             }
 202
 203             //- receive data and update positions if needed
 204             forAll(procBoundaries, patchI)
 205             {
 206                 const label start = procBoundaries[patchI].patchStart();
 207                 labelList receiveData;
 208                 IPstream fromOtherProc
 209                 (
 210                     Pstream::blocking,
 211                     procBoundaries[patchI].neiProcNo()
 212                 );
 213                 fromOtherProc >> receiveData;
 214
 215                 label counter(0);
 216                 while( counter < receiveData.size() )
 217                 {
 218                     const label fI = receiveData[counter++];
 219                     const label pI = receiveData[counter++];
 220
 221                     if( bp[faces[start+fI][pI]] == -1 )
 222                     {
 223                         bp[faces[start+fI][pI]] = pointI++;
 224                         found = true;
 225                     }
 226                 }
 227             }
 228
 229             reduce(found, maxOp<bool>());
 230
 231         } while( found );
 232     }
 233
 234     if( !boundaryPointsPtr_ )
 235         boundaryPointsPtr_ = new labelList();
 236
 237     labelList& boundaryPoints = *boundaryPointsPtr_;
 238     boundaryPoints.setSize(pointI);
 239
 240     //- fill the boundaryPoints list
 241     # ifdef USE_OMP
 242     # pragma omp parallel for num_threads(nThreads) schedule(static, 1)
 243     # endif
 244     forAll(bp, bpI)
 245     {
 246         if( bp[bpI] != -1 )
 247             boundaryPoints[bp[bpI]] = bpI;
 248     }
 249 }
 250
 251 void meshSurfaceEngine::calculateBoundaryFacePatches() const
 252 {
 253     const faceList::subList& bFaces = this->boundaryFaces();
 254     boundaryFacePatchPtr_ = new labelList(bFaces.size());
 255     labelList& facePatch = *boundaryFacePatchPtr_;
 256
 257     label faceI(0);
 258     const PtrList<boundaryPatch>& boundaries = mesh_.boundaries();
 259     forAll(boundaries, patchI)
 260     {
 261         const label nFaces = boundaries[patchI].patchSize();
 262         for(label patchFaceI=0;patchFaceI<nFaces;++patchFaceI)
 263         {
 264             facePatch[faceI] = patchI;
 265             ++faceI;
 266         }
 267     }
 268 }
 269
 270 void meshSurfaceEngine::calculatePointFaces() const
 271 {
 272     //- fill pointFacesAddr
 273     if( !pointFacesPtr_ )
 274         pointFacesPtr_ = new VRWGraph();
 275     VRWGraph& pointFacesAddr = *pointFacesPtr_;
 276
 277     if( !pointInFacePtr_ )
 278         pointInFacePtr_ = new VRWGraph();
 279     VRWGraph& pointInFaceAddr = *pointInFacePtr_;
 280
 281     const labelList& bPoints = this->boundaryPoints();
 282     const faceList::subList& bFaces = this->boundaryFaces();
 283
 284     //- create boundary points
 285     const labelList& bp = this->bp();
 286
 287     labelLongList npf;
 288
 289     # ifdef USE_OMP
 290     label nThreads = 3 * omp_get_num_procs();
 291     if( bPoints.size() < 1000 )
 292         nThreads = 1;
 293     # else
 294     const label nThreads(1);
 295     # endif
 296
 297     label minRow(INT_MAX), maxRow(0);
 298     List<List<LongList<labelPair> > > dataForOtherThreads(nThreads);
 299
 300     # ifdef USE_OMP
 301     # pragma omp parallel num_threads(nThreads)
 302     # endif
 303     {
 304         # ifdef USE_OMP
 305         const label threadI = omp_get_thread_num();
 306         # else
 307         const label threadI(0);
 308         # endif
 309
 310         List<LongList<labelPair> >& dot = dataForOtherThreads[threadI];
 311         dot.setSize(nThreads);
 312
 313         //- find min and max entry in the graph
 314         //- they are used for assigning ranges of values local for each process
 315         label localMinRow(minRow), localMaxRow(0);
 316         # ifdef USE_OMP
 317         # pragma omp for schedule(static)
 318         # endif
 319         forAll(bFaces, bfI)
 320         {
 321             const face& bf = bFaces[bfI];
 322             forAll(bf, pI)
 323             {
 324                 const label bpI = bp[bf[pI]];
 325                 localMaxRow = Foam::max(localMaxRow, bpI);
 326                 localMinRow = Foam::min(localMinRow, bpI);
 327             }
 328         }
 329
 330         ++localMaxRow;
 331
 332         # ifdef USE_OMP
 333         # pragma omp critical
 334         # endif
 335         {
 336             minRow = Foam::min(minRow, localMinRow);
 337             minRow = Foam::max(minRow, 0);
 338             maxRow = Foam::max(maxRow, localMaxRow);
 339
 340             npf.setSize(maxRow);
 341         }
 342
 343         # ifdef USE_OMP
 344         # pragma omp barrier
 345         # endif
 346
 347         //- initialise appearances
 348         # ifdef USE_OMP
 349         # pragma omp for schedule(static)
 350         # endif
 351         for(label i=0;i<maxRow;++i)
 352             npf[i] = 0;
 353
 354         # ifdef USE_OMP
 355         # pragma omp barrier
 356         # endif
 357
 358         const label range = (maxRow - minRow) / nThreads + 1;
 359         const label localMin = minRow + threadI * range;
 360         const label localMax = Foam::min(localMin + range, maxRow);
 361
 362         //- find the number of appearances of each element in the original graph
 363         # ifdef USE_OMP
 364         # pragma omp for schedule(static)
 365         # endif
 366         forAll(bFaces, bfI)
 367         {
 368             const face& bf = bFaces[bfI];
 369
 370             forAll(bf, pI)
 371             {
 372                 const label bpI = bp[bf[pI]];
 373
 374                 const label threadNo = (bpI - minRow) / range;
 375
 376                 if( threadNo == threadI )
 377                 {
 378                     ++npf[bpI];
 379                 }
 380                 else
 381                 {
 382                     dot[threadNo].append(labelPair(bpI, bfI));
 383                 }
 384             }
 385         }
 386
 387         # ifdef USE_OMP
 388         # pragma omp barrier
 389         # endif
 390
 391         //- count the appearances which are not local to the processor
 392         for(label i=0;i<nThreads;++i)
 393         {
 394             const LongList<labelPair>& data =
 395                 dataForOtherThreads[i][threadI];
 396
 397             forAll(data, j)
 398                 ++npf[data[j].first()];
 399         }
 400
 401         # ifdef USE_OMP
 402         # pragma omp barrier
 403         # endif
 404
 405         //- allocate graph
 406         # ifdef USE_OMP
 407         # pragma omp master
 408         # endif
 409         {
 410             VRWGraphSMPModifier(pointFacesAddr).setSizeAndRowSize(npf);
 411             VRWGraphSMPModifier(pointInFaceAddr).setSizeAndRowSize(npf);
 412         }
 413
 414         # ifdef USE_OMP
 415         # pragma omp barrier
 416         # endif
 417
 418         for(label i=localMin;i<localMax;++i)
 419             npf[i] = 0;
 420
 421         //- start filling reverse addressing graph
 422         //- update data from processors with smaller labels
 423         for(label i=0;i<threadI;++i)
 424         {
 425             const LongList<labelPair>& data =
 426                 dataForOtherThreads[i][threadI];
 427
 428             forAll(data, j)
 429             {
 430                 const label bpI = data[j].first();
 431                 const label bfI = data[j].second();
 432
 433                 pointFacesAddr(bpI, npf[bpI]) = bfI;
 434                 pointInFaceAddr(bpI, npf[bpI]) =
 435                     bFaces[bfI].which(bPoints[bpI]);
 436
 437                 ++npf[bpI];
 438             }
 439         }
 440
 441         //- update data local to the processor
 442         # ifdef USE_OMP
 443         # pragma omp for schedule(static)
 444         # endif
 445         forAll(bFaces, bfI)
 446         {
 447             const face& bf = bFaces[bfI];
 448
 449             forAll(bf, pI)
 450             {
 451                 const label bpI = bp[bf[pI]];
 452
 453                 if( (bpI >= localMin) && (bpI < localMax) )
 454                 {
 455                     pointInFaceAddr(bpI, npf[bpI]) = pI;
 456                     pointFacesAddr(bpI, npf[bpI]++) = bfI;
 457                 }
 458             }
 459         }
 460
 461         //- update data from the processors with higher labels
 462         for(label i=threadI+1;i<nThreads;++i)
 463         {
 464             const LongList<labelPair>& data =
 465                 dataForOtherThreads[i][threadI];
 466
 467             forAll(data, j)
 468             {
 469                 const label bpI = data[j].first();
 470                 const label bfI = data[j].second();
 471
 472                 pointFacesAddr(bpI, npf[bpI]) = bfI;
 473                 pointInFaceAddr(bpI, npf[bpI]) =
 474                     bFaces[bfI].which(bPoints[bpI]);
 475
 476                 ++npf[bpI];
 477             }
 478         }
 479     }
 480
 481     pointFacesAddr.setSize(bPoints.size());
 482     pointInFaceAddr.setSize(bPoints.size());
 483 }
 484
 485 void meshSurfaceEngine::calculatePointPatches() const
 486 {
 487     if( !pointPatchesPtr_ )
 488         pointPatchesPtr_ = new VRWGraph();
 489     VRWGraph& pPatches = *pointPatchesPtr_;
 490
 491     const labelList& facePatch = boundaryFacePatches();
 492     const VRWGraph& pFaces = pointFaces();
 493
 494     # ifdef USE_OMP
 495     const label nThreads = 3 * omp_get_num_procs();
 496     # endif
 497
 498     labelList npPatches(pFaces.size());
 499
 500     # ifdef USE_OMP
 501     # pragma omp parallel num_threads(nThreads)
 502     # endif
 503     {
 504         # ifdef USE_OMP
 505         # pragma omp for schedule(static)
 506         # endif
 507         forAll(npPatches, i)
 508             npPatches[i] = 0;
 509
 510         # ifdef USE_OMP
 511         # pragma omp for schedule(static)
 512         # endif
 513         forAll(pFaces, bpI)
 514         {
 515             DynList<label> pf;
 516             forAllRow(pFaces, bpI, pfI)
 517                 pf.appendIfNotIn(facePatch[pFaces(bpI, pfI)]);
 518
 519             npPatches[bpI] = pf.size();
 520         }
 521
 522         # ifdef USE_OMP
 523         # pragma omp barrier
 524
 525         # pragma omp master
 526         # endif
 527         VRWGraphSMPModifier(pPatches).setSizeAndRowSize(npPatches);
 528
 529         # ifdef USE_OMP
 530         # pragma omp barrier
 531
 532         # pragma omp for schedule(static)
 533         # endif
 534         forAll(pFaces, bpI)
 535         {
 536             DynList<label> pf;
 537             forAllRow(pFaces, bpI, pfI)
 538                 pf.appendIfNotIn(facePatch[pFaces(bpI, pfI)]);
 539
 540             pPatches.setRow(bpI, pf);
 541         }
 542     }
 543
 544     if( Pstream::parRun() )
 545     {
 546         const labelList& globalPointLabel = globalBoundaryPointLabel();
 547         const VRWGraph& bpAtProcs = this->bpAtProcs();
 548         const Map<label>& globalToLocal = globalToLocalBndPointAddressing();
 549
 550         std::map<label, labelLongList> exchangeData;
 551
 552         forAllConstIter(Map<label>, globalToLocal, iter)
 553         {
 554             const label bpI = iter();
 555
 556             forAllRow(bpAtProcs, bpI, procI)
 557             {
 558                 const label neiProc = bpAtProcs(bpI, procI);
 559                 if( neiProc == Pstream::myProcNo() )
 560                     continue;
 561
 562                 if( exchangeData.find(neiProc) == exchangeData.end() )
 563                 {
 564                     exchangeData.insert
 565                     (
 566                         std::make_pair(neiProc, labelLongList())
 567                     );
 568                 }
 569
 570                 labelLongList& dataToSend = exchangeData[neiProc];
 571
 572                 //- prepare data which will be sent
 573                 //- data is sent as follows
 574                 //- 1. global point label
 575                 //- 2. number of local patches for point
 576                 //- 3. patch labels for a given point
 577                 dataToSend.append(globalPointLabel[bpI]);
 578                 dataToSend.append(pPatches.sizeOfRow(bpI));
 579                 forAllRow(pPatches, bpI, patchI)
 580                     dataToSend.append(pPatches(bpI, patchI));
 581             }
 582         }
 583
 584         //- exchange data with other processors
 585         labelLongList receivedData;
 586         help::exchangeMap(exchangeData, receivedData);
 587
 588         label counter(0);
 589         while( counter < receivedData.size() )
 590         {
 591             const label bpI = globalToLocal[receivedData[counter++]];
 592             const label nPatches = receivedData[counter++];
 593             for(label i=0;i<nPatches;++i)
 594                 pPatches.appendIfNotIn(bpI, receivedData[counter++]);
 595         }
 596     }
 597 }
 598
 599 void meshSurfaceEngine::calculatePointPoints() const
 600 {
 601     if( !pointPointsPtr_ )
 602         pointPointsPtr_ = new VRWGraph();
 603
 604     VRWGraph& pointPoints = *pointPointsPtr_;
 605
 606     const labelList& boundaryPoints = this->boundaryPoints();
 607     const faceList::subList& bFaces = this->boundaryFaces();
 608     const VRWGraph& pFaces = this->pointFaces();
 609     const labelList& bp = this->bp();
 610
 611     # ifdef USE_OMP
 612     const label nThreads = 3 * omp_get_num_procs();
 613     # endif
 614
 615     labelList npp(boundaryPoints.size());
 616
 617     # ifdef USE_OMP
 618     # pragma omp parallel num_threads(nThreads)
 619     # endif
 620     {
 621         # ifdef USE_OMP
 622         # pragma omp for schedule(static)
 623         # endif
 624         forAll(npp, i)
 625             npp[i] = 0;
 626
 627         # ifdef USE_OMP
 628         # pragma omp for schedule(static)
 629         # endif
 630         forAll(pFaces, bpI)
 631         {
 632             DynList<label> pPoints;
 633
 634             forAllRow(pFaces, bpI, pfI)
 635             {
 636                 const face& bf = bFaces[pFaces(bpI, pfI)];
 637
 638                 const label pos = bf.which(boundaryPoints[bpI]);
 639
 640                 pPoints.appendIfNotIn(bp[bf.nextLabel(pos)]);
 641                 pPoints.appendIfNotIn(bp[bf.prevLabel(pos)]);
 642             }
 643
 644             npp[bpI] = pPoints.size();
 645         }
 646
 647         # ifdef USE_OMP
 648         # pragma omp barrier
 649
 650         # pragma omp master
 651         # endif
 652         VRWGraphSMPModifier(pointPoints).setSizeAndRowSize(npp);
 653
 654         # ifdef USE_OMP
 655         # pragma omp barrier
 656
 657         # pragma omp for schedule(static)
 658         # endif
 659         forAll(pFaces, bpI)
 660         {
 661             DynList<label> pPoints;
 662
 663             forAllRow(pFaces, bpI, pfI)
 664             {
 665                 const face& bf = bFaces[pFaces(bpI, pfI)];
 666
 667                 const label pos = bf.which(boundaryPoints[bpI]);
 668
 669                 pPoints.appendIfNotIn(bp[bf.nextLabel(pos)]);
 670                 pPoints.appendIfNotIn(bp[bf.prevLabel(pos)]);
 671             }
 672
 673             pointPoints.setRow(bpI, pPoints);
 674         }
 675     }
 676
 677     if( Pstream::parRun() )
 678     {
 679         //- this is needed to make the connection matrix symmetric
 680         //- on all processors. In some cases the points on a given processor
 681         //- may not be connected because of a single layer of faces on some
 682         //- other processor. P0, P0 | P1 | P0 P0
 683         const labelList& globalPointLabel = this->globalBoundaryPointLabel();
 684         const Map<label>& globalToLocal =
 685             this->globalToLocalBndPointAddressing();
 686         const VRWGraph& bpAtProcs = this->bpAtProcs();
 687
 688         std::map<label, labelLongList> exchangeData;
 689         forAllConstIter(Map<label>, globalToLocal, iter)
 690         {
 691             const label bpI = iter();
 692
 693             DynList<label> neiToSend;
 694             forAllRow(pointPoints, bpI, j)
 695             {
 696                 const label bpJ = pointPoints(bpI, j);
 697                 if( bpAtProcs.sizeOfRow(bpJ) != 0 )
 698                     neiToSend.append(bpJ);
 699             }
 700
 701             forAllRow(bpAtProcs, bpI, procI)
 702             {
 703                 const label neiProc = bpAtProcs(bpI, procI);
 704                 if( neiProc == Pstream::myProcNo() )
 705                     continue;
 706
 707                 if( exchangeData.find(neiProc) == exchangeData.end() )
 708                     exchangeData.insert(std::make_pair(neiProc,labelLongList()));
 709
 710                 if( neiToSend.size() != 0 )
 711                 {
 712                     labelLongList& dts = exchangeData[neiProc];
 713                     dts.append(globalPointLabel[bpI]);
 714                     dts.append(neiToSend.size());
 715                     forAll(neiToSend, i)
 716                         dts.append(globalPointLabel[neiToSend[i]]);
 717                 }
 718             }
 719         }
 720
 721         labelLongList receivedData;
 722         help::exchangeMap(exchangeData, receivedData);
 723
 724         label counter(0);
 725         while( counter < receivedData.size() )
 726         {
 727             const label bpI = globalToLocal[receivedData[counter++]];
 728             const label size = receivedData[counter++];
 729             for(label i=0;i<size;++i)
 730             {
 731                 const label gpI = receivedData[counter++];
 732                 if( globalToLocal.found(gpI) )
 733                     pointPoints.appendIfNotIn(bpI, globalToLocal[gpI]);
 734             }
 735         }
 736     }
 737 }
 738
 739 void meshSurfaceEngine::calculatePointNormals() const
 740 {
 741     const VRWGraph& pFaces = pointFaces();
 742     const vectorField& fNormals = faceNormals();
 743
 744     pointNormalsPtr_ = new vectorField(pFaces.size());
 745
 746     # ifdef USE_OMP
 747     # pragma omp parallel for if( pFaces.size() > 1000 ) schedule(dynamic, 50)
 748     # endif
 749     forAll(pFaces, pI)
 750     {
 751         vector normal(vector::zero);
 752
 753         forAllRow(pFaces, pI, pfI)
 754             normal += fNormals[pFaces(pI, pfI)];
 755
 756         const scalar d = mag(normal);
 757         if( d > VSMALL )
 758         {
 759             normal /= d;
 760         }
 761         else
 762         {
 763             normal = vector::zero;
 764         }
 765
 766         (*pointNormalsPtr_)[pI] = normal;
 767     }
 768
 769     updatePointNormalsAtProcBoundaries();
 770 }
 771
 772 void meshSurfaceEngine::calculateFaceNormals() const
 773 {
 774     const faceList::subList& bFaces = this->boundaryFaces();
 775     const pointFieldPMG& points = mesh_.points();
 776
 777     faceNormalsPtr_ = new vectorField(bFaces.size());
 778
 779     # ifdef USE_OMP
 780     # pragma omp parallel for if( bFaces.size() > 1000 )
 781     # endif
 782     forAll(bFaces, bfI)
 783     {
 784         const face& bf = bFaces[bfI];
 785
 786         faceNormalsPtr_->operator[](bfI) = bf.normal(points);
 787     }
 788 }
 789
 790 void meshSurfaceEngine::calculateFaceCentres() const
 791 {
 792     const faceList::subList& bFaces = this->boundaryFaces();
 793     const pointFieldPMG& points = mesh_.points();
 794
 795     faceCentresPtr_ = new vectorField(bFaces.size());
 796
 797     # ifdef USE_OMP
 798     # pragma omp parallel for if( bFaces.size() > 1000 )
 799     # endif
 800     forAll(bFaces, bfI)
 801         faceCentresPtr_->operator[](bfI) = bFaces[bfI].centre(points);
 802 }
 803
 804 void meshSurfaceEngine::updatePointNormalsAtProcBoundaries() const
 805 {
 806     if( !Pstream::parRun() )
 807         return;
 808
 809     const VRWGraph& pFaces = pointFaces();
 810     const vectorField& fNormals = faceNormals();
 811     const labelList& globalPointLabel = this->globalBoundaryPointLabel();
 812     const Map<label>& globalToLocal =
 813         this->globalToLocalBndPointAddressing();
 814     const VRWGraph& bpAtProcs = this->bpAtProcs();
 815
 816     vectorField& pNormals = *pointNormalsPtr_;
 817
 818     //- create data which will be sent to other processors
 819     std::map<label, LongList<labelledPoint> > exchangeData;
 820
 821     forAllConstIter(Map<label>, globalToLocal, iter)
 822     {
 823         const label bpI = iter();
 824
 825         vector& n = pNormals[bpI];
 826         n = vector::zero;
 827
 828         forAllRow(pFaces, bpI, pfI)
 829             n += fNormals[pFaces(bpI, pfI)];
 830
 831         forAllRow(bpAtProcs, bpI, procI)
 832         {
 833             const label neiProc = bpAtProcs(bpI, procI);
 834             if( neiProc == Pstream::myProcNo() )
 835                 continue;
 836             if( exchangeData.find(neiProc) == exchangeData.end() )
 837             {
 838                 exchangeData.insert
 839                 (
 840                     std::make_pair(neiProc, LongList<labelledPoint>())
 841                 );
 842             }
 843
 844             exchangeData[neiProc].append
 845             (
 846                 labelledPoint(globalPointLabel[bpI], n)
 847             );
 848         }
 849     }
 850
 851     //- exchange data with other procs
 852     LongList<labelledPoint> receivedData;
 853     help::exchangeMap(exchangeData, receivedData);
 854
 855     forAll(receivedData, i)
 856     {
 857         const label bpI = globalToLocal[receivedData[i].pointLabel()];
 858         pNormals[bpI] += receivedData[i].coordinates();
 859     }
 860
 861     //- normalize vectors
 862     # ifdef USE_OMP
 863     # pragma omp parallel for if( bpAtProcs.size() > 1000 ) \
 864     schedule(guided)
 865     # endif
 866     forAll(bpAtProcs, bpI)
 867     {
 868         if( bpAtProcs.sizeOfRow(bpI) == 0 )
 869             continue;
 870
 871         vector normal = pNormals[bpI];
 872         const scalar d = mag(normal);
 873         if( d > VSMALL )
 874         {
 875             normal /= d;
 876         }
 877         else
 878         {
 879             normal = vector::zero;
 880         }
 881
 882         pNormals[bpI] = normal;
 883     }
 884 }
 885
 886 void meshSurfaceEngine::calculateEdgesAndAddressing() const
 887 {
 888     const VRWGraph& pFaces = pointFaces();
 889     const faceList::subList& bFaces = boundaryFaces();
 890     const labelList& bp = this->bp();
 891
 892     edgesPtr_ = new edgeList();
 893     edgeList& edges = *edgesPtr_;
 894
 895     bpEdgesPtr_ = new VRWGraph();
 896     VRWGraph& bpEdges = *bpEdgesPtr_;
 897
 898     # ifdef USE_OMP
 899     label nThreads = 3 * omp_get_num_procs();
 900     if( pFaces.size() < 1000 )
 901         nThreads = 1;
 902     # else
 903     const label nThreads(1);
 904     # endif
 905
 906     labelList nEdgesForThread(nThreads);
 907
 908     label edgeI(0);
 909
 910     # ifdef USE_OMP
 911     # pragma omp parallel num_threads(nThreads)
 912     # endif
 913     {
 914         LongList<edge> edgesHelper;
 915
 916         # ifdef USE_OMP
 917         # pragma omp for schedule(static)
 918         # endif
 919         forAll(bFaces, bfI)
 920         {
 921             const face& bf = bFaces[bfI];
 922
 923             forAll(bf, pI)
 924             {
 925                 const edge fe = bf.faceEdge(pI);
 926                 const label bpI = bp[fe.start()];
 927                 const label e = fe.end();
 928
 929                 DynList<label> edgeFaces;
 930
 931                 bool store(true);
 932
 933                 //- find all faces attached to this edge
 934                 //- store the edge in case the face faceI is the face
 935                 //- with the smallest label
 936                 forAllRow(pFaces, bpI, pfI)
 937                 {
 938                     const label ofI = pFaces(bpI, pfI);
 939                     const face& of = bFaces[ofI];
 940
 941                     if( of.which(e) < 0 )
 942                         continue;
 943                     if( ofI < bfI )
 944                     {
 945                         store = false;
 946                         break;
 947                     }
 948
 949                     edgeFaces.append(ofI);
 950                 }
 951
 952                 if( store )
 953                     edgesHelper.append(fe);
 954             }
 955         }
 956
 957         //- this enables other threads to see the number of edges
 958         //- generated by each thread
 959         # ifdef USE_OMP
 960         const label threadI = omp_get_thread_num();
 961         # else
 962         const label threadI(0);
 963         # endif
 964         nEdgesForThread[threadI] = edgesHelper.size();
 965
 966         # ifdef USE_OMP
 967         # pragma omp critical
 968         # endif
 969         edgeI += edgesHelper.size();
 970
 971         # ifdef USE_OMP
 972         # pragma omp barrier
 973
 974         # pragma omp master
 975         # endif
 976         edgesPtr_->setSize(edgeI);
 977
 978         # ifdef USE_OMP
 979         # pragma omp barrier
 980         # endif
 981
 982         //- find the starting position of the edges generated by this thread
 983         //- in the global list of edges
 984         label localStart(0);
 985         for(label i=0;i<threadI;++i)
 986             localStart += nEdgesForThread[i];
 987
 988         //- store edges into the global list
 989         forAll(edgesHelper, i)
 990             edgesPtr_->operator[](localStart+i) = edgesHelper[i];
 991     }
 992
 993     //- set the bpEdges
 994     VRWGraphSMPModifier(bpEdges).reverseAddressing(bp, edges);
 995     bpEdges.setSize(pFaces.size());
 996
 997     if( !Pstream::parRun() )
 998         return;
 999
1000     bool addEdges;
1001     do
1002     {
1003         addEdges = false;
1004
1005         //- mark boundary edges for processors which do not contain
1006         //- boundary faces. This procedure is needed to identify boundary
1007         //- edges which are not part of any boundary face on their processor
1008         const faceListPMG& faces = mesh_.faces();
1009         const PtrList<processorBoundaryPatch>& procBoundaries =
1010             mesh_.procBoundaries();
1011
1012         //- send boundary edges to neighbour processors
1013         forAll(procBoundaries, patchI)
1014         {
1015             const label start = procBoundaries[patchI].patchStart();
1016             const label end = start + procBoundaries[patchI].patchSize();
1017
1018             labelLongList dts;
1019             for(label faceI=start;faceI<end;++faceI)
1020             {
1021                 const face& f = faces[faceI];
1022                 forAll(f, eI)
1023                 {
1024                     const edge e = f.faceEdge(eI);
1025                     const label s = bp[e.start()];
1026                     if( s < 0 )
1027                         continue;
1028
1029                     forAllRow(bpEdges, s, peI)
1030                         if( edges[bpEdges(s, peI)] == e )
1031                         {
1032                             dts.append(faceI-start);
1033                             dts.append((f.size()-1-eI)%f.size());
1034                             break;
1035                         }
1036                 }
1037             }
1038
1039             OPstream toOtherProc
1040             (
1041                 Pstream::blocking,
1042                 procBoundaries[patchI].neiProcNo(),
1043                 dts.byteSize()
1044             );
1045             toOtherProc << dts;
1046         }
1047
1048         //- receive data from other processors. Mark edges which are not yet
1049         //- marked as boundary edges
1050         forAll(procBoundaries, patchI)
1051         {
1052             labelList receivedEdges;
1053             IPstream fromOtherProc
1054             (
1055                 Pstream::blocking,
1056                 procBoundaries[patchI].neiProcNo()
1057             );
1058             fromOtherProc >> receivedEdges;
1059
1060             const label start = procBoundaries[patchI].patchStart();
1061             label nReceivedEdges(0);
1062             while( nReceivedEdges < receivedEdges.size() )
1063             {
1064                 const face& f = faces[start+receivedEdges[nReceivedEdges++]];
1065                 const label eI = receivedEdges[nReceivedEdges++];
1066
1067                 const edge e = f.faceEdge(eI);
1068                 const label s = bp[e.start()];
1069
1070                 bool found(false);
1071                 forAllRow(bpEdges, s, peI)
1072                     if( edges[bpEdges(s, peI)] == e )
1073                     {
1074                         found = true;
1075                         break;
1076                     }
1077
1078                 if( !found )
1079                 {
1080                     //- create a new edge
1081                     addEdges = true;
1082                     edges.newElmt(edgeI) = e;
1083
1084                     bpEdges.append(bp[e.start()], edgeI);
1085                     bpEdges.append(bp[e.end()], edgeI);
1086                     ++edgeI;
1087                 }
1088             }
1089         }
1090
1091         reduce(addEdges, maxOp<bool>());
1092     } while( addEdges );
1093
1094     edges.setSize(edgeI);
1095 }
1096
1097 void meshSurfaceEngine::calculateFaceEdgesAddressing() const
1098 {
1099     const faceList::subList& bFaces = this->boundaryFaces();
1100     const labelList& bp = this->bp();
1101     const edgeList& edges = this->edges();
1102     const VRWGraph& pointFaces = this->pointFaces();
1103
1104     faceEdgesPtr_ = new VRWGraph(bFaces.size());
1105     VRWGraph& faceEdges = *faceEdgesPtr_;
1106
1107     labelList nfe(bFaces.size());
1108
1109     # ifdef USE_OMP
1110     const label nThreads = 3 * omp_get_num_procs();
1111
1112     # pragma omp parallel num_threads(nThreads)
1113     # endif
1114     {
1115         # ifdef USE_OMP
1116         # pragma omp for schedule(static)
1117         # endif
1118         forAll(bFaces, bfI)
1119             nfe[bfI] = bFaces[bfI].size();
1120
1121         # ifdef USE_OMP
1122         # pragma omp barrier
1123
1124         # pragma omp master
1125         # endif
1126         VRWGraphSMPModifier(faceEdges).setSizeAndRowSize(nfe);
1127
1128         # ifdef USE_OMP
1129         # pragma omp barrier
1130
1131         # pragma omp for schedule(static)
1132         # endif
1133         forAll(edges, edgeI)
1134         {
1135             const edge ee = edges[edgeI];
1136             const label bpI = bp[ee.start()];
1137
1138             forAllRow(pointFaces, bpI, pfI)
1139             {
1140                 const label bfI = pointFaces(bpI, pfI);
1141
1142                 const face& bf = bFaces[bfI];
1143                 forAll(bf, eI)
1144                 {
1145                     if( bf.faceEdge(eI) == ee )
1146                     {
1147                         faceEdges[bfI][eI] = edgeI;
1148                         break;
1149                     }
1150                 }
1151             }
1152         }
1153     }
1154 }
1155
1156 void meshSurfaceEngine::calculateEdgeFacesAddressing() const
1157 {
1158     const faceList::subList& bFaces = this->boundaryFaces();
1159     const VRWGraph& pointFaces = this->pointFaces();
1160     const edgeList& edges = this->edges();
1161     const labelList& bp = this->bp();
1162
1163     edgeFacesPtr_ = new VRWGraph();
1164     VRWGraph& edgeFaces = *edgeFacesPtr_;
1165
1166     labelList nef(edges.size());
1167
1168     # ifdef USE_OMP
1169     const label nThreads = 3 * omp_get_num_procs();
1170
1171     # pragma omp parallel num_threads(nThreads)
1172     # endif
1173     {
1174         # ifdef USE_OMP
1175         # pragma omp for schedule(static)
1176         # endif
1177         forAll(nef, edgeI)
1178             nef[edgeI] = 0;
1179
1180         # ifdef USE_OMP
1181         # pragma omp for schedule(static)
1182         # endif
1183         forAll(edges, edgeI)
1184         {
1185             const edge& ee = edges[edgeI];
1186             const label bpI = bp[ee.start()];
1187
1188             forAllRow(pointFaces, bpI, pfI)
1189             {
1190                 const label bfI = pointFaces(bpI, pfI);
1191
1192                 const face& bf = bFaces[bfI];
1193
1194                 forAll(bf, eI)
1195                 {
1196                     if( bf.faceEdge(eI) == ee )
1197                     {
1198                         ++nef[edgeI];
1199                         break;
1200                     }
1201                 }
1202             }
1203         }
1204
1205         # ifdef USE_OMP
1206         # pragma omp barrier
1207
1208         # pragma omp master
1209         # endif
1210         VRWGraphSMPModifier(edgeFaces).setSizeAndRowSize(nef);
1211
1212         # ifdef USE_OMP
1213         # pragma omp barrier
1214
1215         # pragma omp for schedule(static)
1216         # endif
1217         forAll(edges, edgeI)
1218         {
1219             const edge& ee = edges[edgeI];
1220             const label bpI = bp[ee.start()];
1221
1222             DynList<label> eFaces;
1223             forAllRow(pointFaces, bpI, pfI)
1224             {
1225                 const label bfI = pointFaces(bpI, pfI);
1226
1227                 const face& bf = bFaces[bfI];
1228
1229                 forAll(bf, eI)
1230                 {
1231                     if( bf.faceEdge(eI) == ee )
1232                     {
1233                         eFaces.append(bfI);
1234                         break;
1235                     }
1236                 }
1237             }
1238
1239             edgeFaces.setRow(edgeI, eFaces);
1240         }
1241     }
1242 }
1243
1244 void meshSurfaceEngine::calculateEdgePatchesAddressing() const
1245 {
1246     edgePatchesPtr_ = new VRWGraph();
1247     VRWGraph& edgePatches = *edgePatchesPtr_;
1248
1249     const VRWGraph& edgeFaces = this->edgeFaces();
1250     const labelList& facePatch = this->boundaryFacePatches();
1251
1252     edgePatches.setSize(edgeFaces.size());
1253
1254     forAll(edgeFaces, eI)
1255     {
1256         DynList<label> ePatches;
1257
1258         forAllRow(edgeFaces, eI, i)
1259         {
1260             const label patchI = facePatch[edgeFaces(eI, i)];
1261
1262             ePatches.appendIfNotIn(patchI);
1263         }
1264
1265         edgePatches.setRow(eI, ePatches);
1266     }
1267
1268     if( Pstream::parRun() )
1269     {
1270         const Map<label>& globalToLocal = globalToLocalBndEdgeAddressing();
1271         const Map<label>& otherPatch = this->otherEdgeFacePatch();
1272
1273         forAllConstIter(Map<label>, globalToLocal, it)
1274         {
1275             const label beI = it();
1276
1277             edgePatches.appendIfNotIn(beI, otherPatch[beI]);
1278         }
1279     }
1280 }
1281
1282 void meshSurfaceEngine::calculateFaceFacesAddressing() const
1283 {
1284     const VRWGraph& edgeFaces = this->edgeFaces();
1285
1286     const faceList::subList& bFaces = boundaryFaces();
1287     faceFacesPtr_ = new VRWGraph(bFaces.size());
1288     VRWGraph& faceFaces = *faceFacesPtr_;
1289
1290     forAll(bFaces, bfI)
1291         faceFaces.setRowSize(bfI, bFaces[bfI].size());
1292
1293     labelList nAppearances(bFaces.size(), 0);
1294
1295     forAll(edgeFaces, efI)
1296     {
1297         if( edgeFaces.sizeOfRow(efI) == 2 )
1298         {
1299             const label f0 = edgeFaces(efI, 0);
1300             const label f1 = edgeFaces(efI, 1);
1301
1302             faceFaces(f0, nAppearances[f0]++) = f1;
1303             faceFaces(f1, nAppearances[f1]++) = f0;
1304         }
1305         else if( Pstream::parRun() && (edgeFaces.sizeOfRow(efI) == 1) )
1306         {
1307             const label f0 = edgeFaces(efI, 0);
1308             faceFaces(f0, nAppearances[f0]++) = -1;
1309         }
1310         else if( Pstream::parRun() && (edgeFaces.sizeOfRow(efI) != 0 ) )
1311         {
1312             FatalErrorIn
1313             (
1314                 "void meshSurfaceEngine::calculateFaceFacesAddressing() const"
1315             ) << "The surface of the mesh is invalid!"
1316                 << " The number of faces containing edge " << efI
1317                 << " is " << edgeFaces.sizeOfRow(efI)
1318                 << " Cannot continue" << exit(FatalError);
1319         }
1320     }
1321 }
1322
1323 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1324
1325 } // End namespace Foam
1326
1327 // ************************************************************************* //