src/mesh/cfMesh/utilities/boundaryLayers/detectBoundaryLayers/detectBoundaryLayersFunctions.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 "detectBoundaryLayers.H"
  29 #include "meshSurfacePartitioner.H"
  30 #include "helperFunctions.H"
  31 #include "polyMeshGenAddressing.H"
  32 #include "polyMeshGen2DEngine.H"
  33 #include "VRWGraphList.H"
  34
  35 #include "labelledPair.H"
  36 #include "labelledScalar.H"
  37
  38 # ifdef USE_OMP
  39 #include <omp.h>
  40 # endif
  41
  42 //#define DEBUGLayer
  43
  44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  45
  46 namespace Foam
  47 {
  48
  49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  50
  51 namespace bndLayerOps
  52 {
  53
  54 class meshBndLayerNeighbourOperator
  55 {
  56     const meshSurfaceEngine& mse_;
  57     const label size_;
  58
  59 public:
  60
  61     meshBndLayerNeighbourOperator(const meshSurfaceEngine& mse)
  62     :
  63         mse_(mse),
  64         size_(mse.boundaryFaces().size())
  65     {}
  66
  67     label size() const
  68     {
  69         return size_;
  70     }
  71
  72     void operator()(const label bfI, DynList<label>& neighbourFaces) const
  73     {
  74         neighbourFaces.clear();
  75
  76         const cellListPMG& cells = mse_.cells();
  77
  78         const labelList& faceOwner = mse_.faceOwners();
  79         const label cellI = faceOwner[bfI];
  80         const cell& c = cells[cellI];
  81
  82         const VRWGraph& faceEdges = mse_.faceEdges();
  83         const VRWGraph& edgeFaces = mse_.edgeFaces();
  84
  85         forAllRow(faceEdges, bfI, feI)
  86         {
  87             const label edgeI = faceEdges(bfI, feI);
  88
  89             if( edgeFaces.sizeOfRow(edgeI) == 2 )
  90             {
  91                 label nei = edgeFaces(edgeI, 0);
  92
  93                 if( nei == bfI )
  94                     nei = edgeFaces(edgeI, 1);
  95
  96                 //- faces must not be part of the same cell
  97                 if( faceOwner[nei] == cellI )
  98                     continue;
  99
 100                 //- owner cell of the other face must
 101                 //- have cellI as its neighbour
 102                 const cell& neiC = cells[faceOwner[nei]];
 103                 bool sharedFace(false);
 104                 forAll(c, fI)
 105                 {
 106                     forAll(neiC, fJ)
 107                     {
 108                         if( c[fI] == neiC[fJ] )
 109                         {
 110                             sharedFace = true;
 111                             break;
 112                         }
 113                     }
 114
 115                     if( sharedFace )
 116                         break;
 117                 }
 118
 119                 if( sharedFace )
 120                     neighbourFaces.append(nei);
 121             }
 122         }
 123     }
 124
 125     template<class labelListType>
 126     void collectGroups
 127     (
 128         std::map<label, DynList<label> >& neiGroups,
 129         const labelListType& elementInGroup,
 130         const DynList<label>& localGroupLabel
 131     ) const
 132     {
 133         const polyMeshGen& mesh = mse_.mesh();
 134         const faceListPMG& faces = mesh.faces();
 135         const cellListPMG& cells = mesh.cells();
 136
 137         const edgeList& edges = mse_.edges();
 138         const labelList& faceOwner = mse_.faceOwners();
 139         const VRWGraph& edgeFaces = mse_.edgeFaces();
 140         const Map<label>& otherProc = mse_.otherEdgeFaceAtProc();
 141         const Map<label>& globalToLocal =
 142             mse_.globalToLocalBndEdgeAddressing();
 143
 144         std::map<label, LongList<labelPair> > exchangeData;
 145         forAll(mse_.beNeiProcs(), procI)
 146             exchangeData[mse_.beNeiProcs()[procI]].clear();
 147
 148         forAllConstIter(Map<label>, globalToLocal, it)
 149         {
 150             const label beI = it();
 151
 152             //- combine data if the cell attached to this face has a face
 153             //- attached to the inter-processor boundary
 154             //- this must hold for boundary layer cells
 155             const cell& c = cells[faceOwner[edgeFaces(beI, 0)]];
 156
 157             bool validCell(false);
 158             forAll(c, fI)
 159             {
 160                 const face& f = faces[c[fI]];
 161
 162                 forAll(f, eI)
 163                 {
 164                     const edge fe = f.faceEdge(eI);
 165
 166                     if
 167                     (
 168                         (fe == edges[beI]) &&
 169                         (mesh.faceIsInProcPatch(c[fI]) >= 0)
 170                     )
 171                     {
 172                         validCell = true;
 173                         break;
 174                     }
 175                 }
 176
 177                 if( validCell )
 178                     break;
 179             }
 180
 181             if( !validCell )
 182                 continue;
 183
 184             const label groupI = elementInGroup[edgeFaces(beI, 0)];
 185
 186             if( groupI < 0 )
 187                 continue;
 188
 189             const label lgI = localGroupLabel[groupI];
 190             exchangeData[otherProc[beI]].append(labelPair(it.key(), lgI));
 191         }
 192
 193         LongList<labelPair> receivedData;
 194         help::exchangeMap(exchangeData, receivedData);
 195
 196         forAll(receivedData, i)
 197         {
 198             const labelPair& lp = receivedData[i];
 199
 200             const label beI = globalToLocal[lp.first()];
 201
 202             const cell& c = cells[faceOwner[edgeFaces(beI, 0)]];
 203
 204             //- combine data if the cell attached to this face has a face
 205             //- attached to the inter-processor boundary
 206             //- this must hold for boundary layer cells
 207             bool validCell(false);
 208             forAll(c, fI)
 209             {
 210                 const face& f = faces[c[fI]];
 211
 212                 forAll(f, eI)
 213                 {
 214                     const edge fe = f.faceEdge(eI);
 215
 216                     if
 217                     (
 218                         (fe == edges[beI]) &&
 219                         (mesh.faceIsInProcPatch(c[fI]) >= 0)
 220                     )
 221                     {
 222                         validCell = true;
 223                         break;
 224                     }
 225                 }
 226
 227                 if( validCell )
 228                     break;
 229             }
 230
 231             if( !validCell )
 232                 continue;
 233
 234             const label groupI = elementInGroup[edgeFaces(beI, 0)];
 235
 236             if( groupI < 0 )
 237                 continue;
 238
 239             DynList<label>& ng = neiGroups[localGroupLabel[groupI]];
 240
 241             //- store the connection over the inter-processor boundary
 242             ng.appendIfNotIn(lp.second());
 243         }
 244     }
 245 };
 246
 247 class meshBndLayerSelectorOperator
 248 {
 249     const meshSurfaceEngine& mse_;
 250
 251 public:
 252
 253     meshBndLayerSelectorOperator(const meshSurfaceEngine& mse)
 254     :
 255         mse_(mse)
 256     {}
 257
 258     bool operator()(const label bfI) const
 259     {
 260         const labelList& faceOwner = mse_.faceOwners();
 261         const polyMeshGen& mesh = mse_.mesh();
 262         const faceListPMG& faces = mesh.faces();
 263
 264         const cell& c = mesh.cells()[faceOwner[bfI]];
 265         const PtrList<boundaryPatch>& boundaries = mesh.boundaries();
 266         const label start = boundaries[0].patchStart();
 267
 268         label nBndFaces(0), baseFace(-1), otherBase(-1), nQuads(0);
 269         forAll(c, fI)
 270         {
 271             if( faces[c[fI]].size() == 4 )
 272                 ++nQuads;
 273
 274             if( (c[fI] - start) == bfI )
 275             {
 276                 baseFace = fI;
 277                 ++nBndFaces;
 278             }
 279         }
 280
 281         if( nQuads == 6 )
 282         {
 283             //- cell is a hex
 284             return true;
 285         }
 286
 287         if( (nQuads + 2) != c.size() )
 288             return false;
 289
 290         if( nBndFaces != 1 )
 291             return false;
 292
 293         label nQuadsAttachedToBaseFace(0);
 294         forAll(c, fI)
 295         {
 296             if( fI == baseFace )
 297                 continue;
 298
 299             const bool sEdge =
 300                 help::shareAnEdge(faces[c[baseFace]], faces[c[fI]]);
 301
 302             if( (faces[c[fI]].size() == 4) && sEdge )
 303             {
 304                 ++nQuadsAttachedToBaseFace;
 305             }
 306             else if( !sEdge )
 307             {
 308                 if( otherBase != -1 )
 309                     return false;
 310
 311                 otherBase = fI;
 312             }
 313         }
 314
 315         if(
 316             (nQuads == 6) ||
 317             (
 318                 ((nQuadsAttachedToBaseFace + 2) == c.size()) &&
 319                 otherBase != -1 &&
 320                 !help::shareAnEdge(faces[c[baseFace]], faces[c[otherBase]])
 321             )
 322         )
 323             return true;
 324
 325         return false;
 326     }
 327 };
 328
 329 } // End namespace bndLayerOps
 330
 331 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 332
 333 void detectBoundaryLayers::analyseLayers()
 334 {
 335     Info << "Analysing mesh for bnd layer existence" << endl;
 336
 337     const meshSurfaceEngine& mse = meshSurface_.surfaceEngine();
 338     const polyMeshGen& mesh = mse.mesh();
 339     const PtrList<boundaryPatch>& boundaries = mesh.boundaries();
 340
 341     //- allocate data needed in parallel loops
 342     mse.faceOwners();
 343     mse.faceEdges();
 344     mse.edgeFaces();
 345     mse.edges();
 346     mse.boundaryPointEdges();
 347
 348     //- find layers in patch
 349     nFirstLayers_ =
 350         help::groupMarking
 351         (
 352             layerAtBndFace_,
 353             bndLayerOps::meshBndLayerNeighbourOperator(mse),
 354             bndLayerOps::meshBndLayerSelectorOperator(mse)
 355         );
 356
 357     # ifdef DEBUGLayer
 358     labelList layerSubsetId(nFirstLayers_);
 359     polyMeshGen& pmg = const_cast<polyMeshGen&>(mesh);
 360     forAll(layerSubsetId, i)
 361         layerSubsetId[i] = pmg.addCellSubset("bndLayer"+help::scalarToText(i));
 362
 363
 364     forAll(layerAtBndFace_, bfI)
 365     {
 366         if( layerAtBndFace_[bfI] < 0 )
 367             continue;
 368
 369         pmg.addCellToSubset
 370         (
 371             layerSubsetId[layerAtBndFace_[bfI]],
 372             mse.faceOwners()[bfI]
 373         );
 374     }
 375     # endif
 376
 377     if( is2DMesh_ )
 378     {
 379         polyMeshGen2DEngine mesh2DEngine(mse.mesh());
 380         const boolList& zMinPoint = mesh2DEngine.zMinPoints();
 381         const boolList& zMaxPoint = mesh2DEngine.zMaxPoints();
 382
 383         const faceList::subList& bFaces = mse.boundaryFaces();
 384
 385         # ifdef USE_OMP
 386         # pragma omp parallel for schedule(dynamic, 50)
 387         # endif
 388         forAll(bFaces, bfI)
 389         {
 390             const face& bf = bFaces[bfI];
 391
 392             bool allZMin(true), allZMax(true);
 393             forAll(bf, pI)
 394             {
 395                 if( !zMinPoint[bf[pI]] )
 396                     allZMin = false;
 397                 if( !zMaxPoint[bf[pI]] )
 398                     allZMax = false;
 399             }
 400
 401             if( allZMax ^ allZMin )
 402                 layerAtBndFace_[bfI] = -1;
 403         }
 404     }
 405
 406     //- check the number of layer found at a patch
 407     typedef std::map<label, DynList<label> > patchToLayerType;
 408     patchToLayerType patchToLayer;
 409
 410     const labelList& facePatch = meshSurface_.boundaryFacePatches();
 411
 412     forAll(facePatch, bfI)
 413     {
 414         patchToLayer[facePatch[bfI]].appendIfNotIn(layerAtBndFace_[bfI]);
 415     }
 416
 417     //- all faces of a patch must be in the same layer
 418     layerAtPatch_.setSize(boundaries.size());
 419     forAll(layerAtPatch_, i)
 420         layerAtPatch_[i].clear();
 421
 422     for
 423     (
 424         patchToLayerType::const_iterator it=patchToLayer.begin();
 425         it!=patchToLayer.end();
 426         ++it
 427     )
 428     {
 429         const DynList<label>& layersAtPatch = it->second;
 430
 431         forAll(layersAtPatch, i)
 432         {
 433             if( layersAtPatch[i] < 0 )
 434             {
 435                 layerAtPatch_[it->first].clear();
 436                 break;
 437             }
 438             else
 439             {
 440                 layerAtPatch_[it->first].append(layersAtPatch[i]);
 441             }
 442         }
 443     }
 444
 445     //- set the layer ID to -1 for all faces where the patch is set to -1
 446     forAll(facePatch, bfI)
 447     {
 448         if( layerAtPatch_[facePatch[bfI]].size() == 0 )
 449             layerAtBndFace_[bfI] = -1;
 450     }
 451
 452     # ifdef DEBUGLayer
 453     Info << "Layer at patch " << layerAtPatch_ << endl;
 454     forAll(layerAtBndFace_, bfI)
 455         if( layerAtBndFace_[bfI] < 0 )
 456             Info << "0.2 No layer at boundary face " << bfI << endl;
 457     # endif
 458 }
 459
 460 bool detectBoundaryLayers::findHairsForFace
 461 (
 462     const label bfI,
 463     DynList<edge>& hairEdges
 464 ) const
 465 {
 466     const meshSurfaceEngine& mse = meshSurface_.surfaceEngine();
 467     const polyMeshGen& mesh = mse.mesh();
 468
 469     const label nInternalFaces = mesh.boundaries()[0].patchStart();
 470
 471     const labelList& faceOwner = mse.faceOwners();
 472
 473     const faceListPMG& faces = mesh.faces();
 474     const cell& c = mesh.cells()[faceOwner[bfI]];
 475
 476     //- check cell topology
 477     DynList<edge, 48> edges;
 478     DynList<DynList<label, 2>, 48> edgeFaces;
 479     DynList<DynList<label, 10>, 24> faceEdges;
 480     faceEdges.setSize(c.size());
 481     label baseFace(-1);
 482     forAll(c, fI)
 483     {
 484         if( c[fI] - nInternalFaces == bfI )
 485         {
 486             baseFace = fI;
 487         }
 488
 489         const face& f = faces[c[fI]];
 490         faceEdges[fI].setSize(f.size());
 491
 492         forAll(f, eI)
 493         {
 494             const edge e = f.faceEdge(eI);
 495
 496             label pos = edges.containsAtPosition(e);
 497
 498             if( pos < 0 )
 499             {
 500                 pos = edges.size();
 501                 edges.append(e);
 502                 edgeFaces.setSize(pos+1);
 503             }
 504
 505             edgeFaces[pos].append(fI);
 506             faceEdges[fI][eI] = pos;
 507         }
 508     }
 509
 510     //- check does the base face exist and is the number of faces
 511     //- in the cell corresponding to a prism cell
 512     if( (baseFace < 0) || ((c.size() - faces[c[baseFace]].size()) != 2) )
 513         return false;
 514
 515     //- check if all faces attached to the base face are quads
 516     bool isPrism(true);
 517
 518     const face& bf = faces[c[baseFace]];
 519     hairEdges.setSize(bf.size());
 520
 521     forAll(bf, pI)
 522     {
 523         const label nextEdge = faceEdges[baseFace][pI];
 524         const label prevEdge = faceEdges[baseFace][bf.rcIndex(pI)];
 525
 526         if( edgeFaces[nextEdge].size() != 2 || edgeFaces[prevEdge].size() != 2 )
 527         {
 528             isPrism = false;
 529             break;
 530         }
 531
 532         //- find the face attached to the edge after the current point
 533         label otherNextFace = edgeFaces[nextEdge][0];
 534         if( otherNextFace == baseFace )
 535             otherNextFace = edgeFaces[nextEdge][1];
 536
 537         //- find the face attached to the edge before the current point
 538         label otherPrevFace = edgeFaces[prevEdge][0];
 539         if( otherPrevFace == baseFace )
 540             otherPrevFace = edgeFaces[prevEdge][1];
 541
 542         label commonEdge;
 543         for(commonEdge=0;commonEdge<edges.size();++commonEdge)
 544             if(
 545                 edgeFaces[commonEdge].contains(otherNextFace) &&
 546                 edgeFaces[commonEdge].contains(otherPrevFace)
 547             )
 548                 break;
 549
 550         if( commonEdge == edges.size() )
 551         {
 552             isPrism = false;
 553             break;
 554         }
 555
 556         //- there exists a common edge which shall be used as a hair
 557         if( edges[commonEdge].start() == bf[pI] )
 558         {
 559             hairEdges[pI] = edges[commonEdge];
 560         }
 561         else
 562         {
 563             hairEdges[pI] = edges[commonEdge].reverseEdge();
 564         }
 565     }
 566
 567     return isPrism;
 568 }
 569
 570 void detectBoundaryLayers::generateHairEdges()
 571 {
 572     hairEdges_.clear();
 573     hairEdgesAtBoundaryPoint_.clear();
 574
 575     const meshSurfaceEngine& mse = meshSurface_.surfaceEngine();
 576     const faceList::subList& bFaces = mse.boundaryFaces();
 577     mse.faceOwners();
 578     const VRWGraph& pFaces = mse.pointFaces();
 579     const labelList& bp = mse.bp();
 580
 581     # ifdef USE_OMP
 582     # pragma omp parallel if( bFaces.size() > 1000 )
 583     # endif
 584     {
 585         edgeLongList localEdges;
 586
 587         # ifdef USE_OMP
 588         # pragma omp for schedule(dynamic, 100)
 589         # endif
 590         forAll(bFaces, bfI)
 591         {
 592             if( layerAtBndFace_[bfI] < 0 )
 593                 continue;
 594
 595             //- find hair edges for this face
 596             DynList<edge> hairEdges;
 597             if( !findHairsForFace(bfI, hairEdges) )
 598                 continue;
 599
 600             const face& bf = bFaces[bfI];
 601
 602             forAll(bf, pI)
 603             {
 604                 //- store hair edges in a list
 605                 const edge& he = hairEdges[pI];
 606
 607                 if( he.start() != bf[pI] )
 608                     FatalErrorIn
 609                     (
 610                         "void detectBoundaryLayers::generateHairEdges()"
 611                     ) << "Wrong starting point" << abort(FatalError);
 612
 613                 localEdges.append(he);
 614             }
 615         }
 616
 617         # ifdef USE_OMP
 618         //- find the starting element for this thread
 619         label startEl;
 620         # pragma omp critical
 621         {
 622             startEl = hairEdges_.size();
 623
 624             hairEdges_.setSize(startEl+localEdges.size());
 625         }
 626
 627         # pragma omp barrier
 628
 629         //- copy the local data to splitEdges_
 630         forAll(localEdges, i)
 631             hairEdges_[startEl++] = localEdges[i];
 632         # else
 633         //- just transfer the data to splitEdges_
 634         hairEdges_.transfer(localEdges);
 635         # endif
 636     }
 637
 638     //- filter out duplicate edges
 639     VRWGraph pHairEdges;
 640     pHairEdges.reverseAddressing(hairEdges_);
 641
 642     boolList duplicateEdge(hairEdges_.size(), false);
 643     # ifdef USE_OMP
 644     # pragma omp parallel for schedule(dynamic, 50)
 645     # endif
 646     forAll(pHairEdges, pointI)
 647     {
 648         forAllRow(pHairEdges, pointI, pheI)
 649         {
 650             const label heI = pHairEdges(pointI, pheI);
 651             const edge& he = hairEdges_[heI];
 652
 653             for(label pheJ=pheI+1;pheJ<pHairEdges.sizeOfRow(pointI);++pheJ)
 654             {
 655                 const label heJ = pHairEdges(pointI, pheJ);
 656                 const edge& nhe = hairEdges_[heJ];
 657
 658                 if( he == nhe )
 659                     duplicateEdge[heJ] = true;
 660             }
 661         }
 662     }
 663
 664     label counter(0);
 665     forAll(hairEdges_, heI)
 666     {
 667         if( !duplicateEdge[heI] )
 668         {
 669             if( heI > counter )
 670             {
 671                 hairEdges_[counter++] = hairEdges_[heI];
 672             }
 673             else
 674             {
 675                 ++counter;
 676             }
 677         }
 678     }
 679
 680     hairEdges_.setSize(counter);
 681
 682     //- create point to split edges addressing
 683     hairEdgesAtBoundaryPoint_.setSize(pFaces.size());
 684
 685     forAll(hairEdges_, heI)
 686     {
 687         const edge& he = hairEdges_[heI];
 688         hairEdgesAtBoundaryPoint_.append(bp[he.start()], heI);
 689     }
 690 }
 691
 692 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 693
 694 } // End namespace Foam
 695
 696 // ************************************************************************* //