src/mesh/cfMesh/utilities/surfaceTools/correctEdgesBetweenPatches/correctEdgesBetweenPatchesDistributeFaces.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 "demandDrivenData.H"
  29 #include "correctEdgesBetweenPatches.H"
  30 #include "decomposeFaces.H"
  31 #include "triSurface.H"
  32 #include "meshSurfaceEngine.H"
  33 #include "meshSurfaceCheckEdgeTypes.H"
  34 #include "meshSurfacePartitioner.H"
  35 #include "helperFunctions.H"
  36 #include "helperFunctionsPar.H"
  37
  38 #include <map>
  39
  40 # ifdef USE_OMP
  41 #include <omp.h>
  42 # endif
  43
  44 //#define DEBUGMapping
  45
  46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  47
  48 namespace Foam
  49 {
  50
  51 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  52
  53 void correctEdgesBetweenPatches::decomposeProblematicFaces()
  54 {
  55     Info << "Decomposing problematic faces" << endl;
  56     const meshSurfaceEngine& mse = meshSurface();
  57     const labelList& bp = mse.bp();
  58     const edgeList& edges = mse.edges();
  59     const VRWGraph& pointEdges = mse.boundaryPointEdges();
  60     const VRWGraph& faceEdges = mse.faceEdges();
  61     const VRWGraph& edgeFaces = mse.edgeFaces();
  62     const labelList& facePatches = mse.boundaryFacePatches();
  63
  64     //- mark feature edges
  65     boolList featureBndEdge(edgeFaces.size(), false);
  66
  67     forAll(edgeFaces, beI)
  68     {
  69         if( edgeFaces.sizeOfRow(beI) != 2 )
  70             continue;
  71
  72         if( facePatches[edgeFaces(beI, 0)] != facePatches[edgeFaces(beI, 1)] )
  73             featureBndEdge[beI] = true;
  74     }
  75
  76     if( Pstream::parRun() )
  77     {
  78         //- find feature edges at parallel boundaries and propagate the
  79         //- information to all processors
  80         const Map<label>& globalToLocalEdge =
  81             mse.globalToLocalBndEdgeAddressing();
  82         const VRWGraph& beAtProcs = mse.beAtProcs();
  83         const Map<label>& otherProcPatch = mse.otherEdgeFacePatch();
  84         forAllConstIter(Map<label>, globalToLocalEdge, it)
  85         {
  86             const label beI = it();
  87
  88             if( edgeFaces.sizeOfRow(beI) != 1 )
  89                 continue;
  90             if( facePatches[edgeFaces(beI, 0)] != otherProcPatch[beI] )
  91                 featureBndEdge[beI] = true;
  92         }
  93
  94         //- propagate information to all processors that need this information
  95         std::map<label, labelLongList> exchangeData;
  96         forAll(mse.beNeiProcs(), i)
  97             exchangeData.insert
  98             (
  99                 std::make_pair(mse.beNeiProcs()[i], labelLongList())
 100             );
 101
 102         //- append labels of feature edges that need to be sent to other
 103         //- processors sharing that edge
 104         forAllConstIter(Map<label>, globalToLocalEdge, it)
 105         {
 106             const label beI = it();
 107
 108             if( featureBndEdge[beI] )
 109             {
 110                 forAllRow(beAtProcs, beI, i)
 111                 {
 112                     const label procI = beAtProcs(beI, i);
 113                     if( procI == Pstream::myProcNo() )
 114                         continue;
 115
 116                     exchangeData[procI].append(it.key());
 117                 }
 118             }
 119         }
 120
 121         labelLongList receivedData;
 122         help::exchangeMap(exchangeData, receivedData);
 123
 124         label counter(0);
 125         while( counter < receivedData.size() )
 126         {
 127             featureBndEdge[globalToLocalEdge[receivedData[counter++]]] = true;
 128         }
 129     }
 130
 131     const polyMeshGen& mesh = mse.mesh();
 132     const faceListPMG& faces = mesh.faces();
 133     const labelList& owner = mesh.owner();
 134     const labelList& neighbour = mesh.neighbour();
 135
 136     boolList decomposeFace(faces.size(), false);
 137     label nDecomposedFaces(0);
 138
 139     //- decompose internal faces with more than one feature edge
 140     const label nIntFaces = mesh.nInternalFaces();
 141     # ifdef USE_OMP
 142     # pragma omp parallel for schedule(guided) reduction(+ : nDecomposedFaces)
 143     # endif
 144     for(label faceI=0;faceI<nIntFaces;++faceI)
 145     {
 146         const face& f = faces[faceI];
 147
 148         label nFeatureEdges(0);
 149
 150         forAll(f, eI)
 151         {
 152             const edge e = f.faceEdge(eI);
 153
 154             const label bs = bp[e[0]];
 155             const label be = bp[e[1]];
 156             if( (bs != -1) && (be != -1) )
 157             {
 158                 //- check if this edge is a boundary edges and a feature edge
 159                 forAllRow(pointEdges, bs, i)
 160                 {
 161                     const label beI = pointEdges(bs, i);
 162
 163                     if( (edges[beI] == e) && featureBndEdge[beI] )
 164                         ++nFeatureEdges;
 165                 }
 166             }
 167         }
 168
 169         if( nFeatureEdges > 1 )
 170         {
 171             ++nDecomposedFaces;
 172             decomposeFace[faceI] = true;
 173             decomposeCell_[owner[faceI]] = true;
 174             decomposeCell_[neighbour[faceI]] = true;
 175         }
 176     }
 177
 178     //- decompose boundary faces in case the feature edges are not connected
 179     //- into a single open chain of edges
 180     # ifdef USE_OMP
 181     # pragma omp parallel for schedule(guided) reduction(+ : nDecomposedFaces)
 182     # endif
 183     forAll(faceEdges, bfI)
 184     {
 185         boolList featureEdge(faceEdges.sizeOfRow(bfI), false);
 186
 187         forAllRow(faceEdges, bfI, feI)
 188             if( featureBndEdge[faceEdges(bfI, feI)] )
 189                 featureEdge[feI] = true;
 190
 191         if( !help::areElementsInChain(featureEdge) )
 192         {
 193             ++nDecomposedFaces;
 194             decomposeFace[nIntFaces+bfI] = true;
 195             decomposeCell_[owner[nIntFaces+bfI]] = true;
 196         }
 197     }
 198
 199     if( Pstream::parRun() )
 200     {
 201         //- decompose processor faces having more than one feature edge
 202         const PtrList<processorBoundaryPatch>& procBoundaries =
 203             mesh.procBoundaries();
 204
 205         forAll(procBoundaries, patchI)
 206         {
 207             const label start = procBoundaries[patchI].patchStart();
 208             const label end = start + procBoundaries[patchI].patchSize();
 209
 210             # ifdef USE_OMP
 211             # pragma omp parallel for schedule(guided) \
 212             reduction(+ : nDecomposedFaces)
 213             # endif
 214             for(label faceI=start;faceI<end;++faceI)
 215             {
 216                 const face& f = faces[faceI];
 217
 218                 label nFeatureEdges(0);
 219
 220                 forAll(f, eI)
 221                 {
 222                     const edge e = f.faceEdge(eI);
 223
 224                     const label bs = bp[e[0]];
 225                     const label be = bp[e[1]];
 226                     if( (bs != -1) && (be != -1) )
 227                     {
 228                         //- check if this edge is a boundary edge
 229                         //- and a feature edge
 230                         forAllRow(pointEdges, bs, i)
 231                         {
 232                             const label beI = pointEdges(bs, i);
 233
 234                             if( (edges[beI] == e) && featureBndEdge[beI] )
 235                                 ++nFeatureEdges;
 236                         }
 237                     }
 238                 }
 239
 240                 if( nFeatureEdges > 1 )
 241                 {
 242                     ++nDecomposedFaces;
 243                     decomposeFace[faceI] = true;
 244                     decomposeCell_[owner[faceI]] = true;
 245                 }
 246             }
 247         }
 248     }
 249
 250     reduce(nDecomposedFaces, sumOp<label>());
 251
 252     if( nDecomposedFaces != 0 )
 253     {
 254         Info << nDecomposedFaces << " faces decomposed into triangles" << endl;
 255
 256         decompose_ = true;
 257         decomposeFaces df(mesh_);
 258         df.decomposeMeshFaces(decomposeFace);
 259
 260         clearMeshSurface();
 261         mesh_.clearAddressingData();
 262     }
 263
 264     Info << "Finished decomposing problematic faces" << endl;
 265 }
 266
 267 void correctEdgesBetweenPatches::decomposeConcaveFaces()
 268 {
 269     const meshSurfaceEngine& mse = meshSurface();
 270     const labelList& bPoints = mse.boundaryPoints();
 271     const labelList& bp = mse.bp();
 272     const edgeList& edges = mse.edges();
 273     const VRWGraph& edgeFaces = mse.edgeFaces();
 274     const VRWGraph& bpEdges = mse.boundaryPointEdges();
 275     const labelList& facePatch = mse.boundaryFacePatches();
 276
 277     //- classify edges at the surface
 278     meshSurfaceCheckEdgeTypes edgeChecker(mse);
 279     const List<direction>& edgeType = edgeChecker.edgeTypes();
 280
 281     //- find concave points
 282     boolList concavePoint(bPoints.size(), false);
 283
 284     labelList edgeInPatch(edges.size(), -1);
 285
 286     direction problematicTypes = 0;
 287     problematicTypes |= meshSurfaceCheckEdgeTypes::CONCAVEEDGE;
 288     problematicTypes |= meshSurfaceCheckEdgeTypes::UNDETERMINED;
 289
 290     # ifdef USE_OMP
 291     # pragma omp parallel for schedule(dynamic, 100)
 292     # endif
 293     forAll(edgeType, eI)
 294     {
 295         if( edgeType[eI] & meshSurfaceCheckEdgeTypes::PATCHEDGE )
 296         {
 297             if( edgeFaces.sizeOfRow(eI) )
 298                 edgeInPatch[eI] = facePatch[edgeFaces(eI, 0)];
 299
 300             continue;
 301         }
 302
 303         if( edgeType[eI] & problematicTypes)
 304         {
 305             const edge& e = edges[eI];
 306
 307             concavePoint[bp[e.start()]] = true;
 308             concavePoint[bp[e.end()]] = true;
 309         }
 310     }
 311
 312     if( Pstream::parRun() )
 313     {
 314         const Map<label>& globalToLocal =
 315             mse.globalToLocalBndEdgeAddressing();
 316         const VRWGraph& beAtProcs = mse.beAtProcs();
 317
 318         std::map<label, labelLongList> exchangeData;
 319         forAll(mse.beNeiProcs(), i)
 320             exchangeData[mse.beNeiProcs()[i]].clear();
 321
 322         forAllConstIter(Map<label>, globalToLocal, it)
 323         {
 324             const label beI = it();
 325
 326             if( edgeInPatch[beI] < 0 )
 327                 continue;
 328
 329             forAllRow(beAtProcs, beI, i)
 330             {
 331                 const label neiProc = beAtProcs(beI, i);
 332
 333                 if( neiProc == Pstream::myProcNo() )
 334                     continue;
 335
 336                 labelLongList& dts = exchangeData[neiProc];
 337
 338                 dts.append(it.key());
 339                 dts.append(edgeInPatch[beI]);
 340             }
 341         }
 342
 343         labelLongList receivedData;
 344         help::exchangeMap(exchangeData, receivedData);
 345
 346         for(label i=0;i<receivedData.size();)
 347         {
 348             const label beI = globalToLocal[receivedData[i++]];
 349             const label patchI = receivedData[i++];
 350             if( edgeInPatch[beI] == -1 )
 351             {
 352                 edgeInPatch[beI] = patchI;
 353             }
 354             else if( edgeInPatch[beI] != patchI )
 355             {
 356                 FatalErrorIn
 357                 (
 358                     "void correctEdgesBetweenPatches::decomposeConcaveFaces()"
 359                 ) << "Invalid patch!" << abort(FatalError);
 360             }
 361         }
 362     }
 363
 364     //- decompose internal faces attached to concave vertices which have two
 365     //- or more edges at the boundary
 366     const faceListPMG& faces = mesh_.faces();
 367     const labelList& owner = mesh_.owner();
 368     const labelList& neighbour = mesh_.neighbour();
 369
 370     boolList decomposeFace(faces.size(), false);
 371     label nDecomposed(0);
 372
 373     # ifdef USE_OMP
 374     # pragma omp parallel for schedule(dynamic, 100) \
 375     reduction(+ : nDecomposed)
 376     # endif
 377     for(label faceI=0;faceI<mesh_.nInternalFaces();++faceI)
 378     {
 379         const face& f = faces[faceI];
 380
 381         bool hasConcave(false);
 382         label nBndEdges(0);
 383         DynList<label> bndEdgePatches;
 384
 385         forAll(f, pI)
 386         {
 387             const label bpI = bp[f[pI]];
 388
 389             if( bpI < 0 )
 390                 continue;
 391
 392             if( concavePoint[bpI] )
 393                 hasConcave = true;
 394
 395             //- points is at a concave edge
 396             //- count the number of boundary edge
 397             const edge e = f.faceEdge(pI);
 398
 399             forAllRow(bpEdges, bpI, bpeI)
 400             {
 401                 const label beI = bpEdges(bpI, bpeI);
 402                 const edge& ee = edges[beI];
 403
 404                 if( e == ee )
 405                 {
 406                     ++nBndEdges;
 407                     bndEdgePatches.appendIfNotIn(edgeInPatch[beI]);
 408                     break;
 409                 }
 410             }
 411         }
 412
 413         if( hasConcave && (nBndEdges > 1) && (bndEdgePatches.size() > 1) )
 414         {
 415             //- the face has two or more edges at the boundary
 416             //- Hence, it is marked for decomposition
 417             decomposeFace[faceI] = true;
 418
 419             decomposeCell_[owner[faceI]] = true;
 420             decomposeCell_[neighbour[faceI]] = true;
 421
 422             ++nDecomposed;
 423         }
 424     }
 425
 426     //- finally, perform decomposition of marked faces
 427     if( returnReduce(nDecomposed, sumOp<label>()) != 0 )
 428     {
 429         Info << "Decomposing " << nDecomposed << " internal faces" << endl;
 430         decomposeFaces(mesh_).decomposeMeshFaces(decomposeFace);
 431
 432         decompose_ = true;
 433
 434         clearMeshSurface();
 435         mesh_.clearAddressingData();
 436     }
 437 }
 438
 439 void correctEdgesBetweenPatches::patchCorrection()
 440 {
 441     Info << "Performing patch correction" << endl;
 442
 443     const meshSurfaceEngine& mse = meshSurface();
 444
 445     meshSurfacePartitioner surfacePartitioner(mse);
 446
 447     const labelList& bPoints = mse.boundaryPoints();
 448     const VRWGraph& faceEdges = mse.faceEdges();
 449     const VRWGraph& edgeFaces = mse.edgeFaces();
 450     const faceList::subList& bFaces = mse.boundaryFaces();
 451     const labelList& bp = mse.bp();
 452     const labelList& boundaryFaceOwners = mse.faceOwners();
 453     const labelList& facePatches = mse.boundaryFacePatches();
 454
 455     //- set flag 1 to corner vertices, flag 2 to edge vertices
 456     List<direction> nodeType(bPoints.size(), direction(0));
 457
 458     //- set corner flags
 459     const labelHashSet& corners = surfacePartitioner.corners();
 460     forAllConstIter(labelHashSet, corners, it)
 461         nodeType[it.key()] |= 1;
 462
 463     //- set flgs to edge vertices
 464     const labelHashSet& edgePoints = surfacePartitioner.edgePoints();
 465     forAllConstIter(labelHashSet, edgePoints, it)
 466         nodeType[it.key()] |= 2;
 467
 468     //- set flags for feature edges
 469     boolList featureEdge(edgeFaces.size(), false);
 470
 471     # ifdef USE_OMP
 472     # pragma omp parallel for schedule(guided)
 473     # endif
 474     forAll(edgeFaces, eI)
 475     {
 476         if( edgeFaces.sizeOfRow(eI) != 2 )
 477             continue;
 478
 479         if( facePatches[edgeFaces(eI, 0)] != facePatches[edgeFaces(eI, 1)] )
 480             featureEdge[eI] = true;
 481     }
 482
 483     if( Pstream::parRun() )
 484     {
 485         //- set flags for edges at parallel boundaries
 486         const Map<label>& otherProcPatch = mse.otherEdgeFacePatch();
 487
 488         forAllConstIter(Map<label>, otherProcPatch, it)
 489             if( facePatches[edgeFaces(it.key(), 0)] != it() )
 490                 featureEdge[it.key()] = true;
 491     }
 492
 493     //- decompose bad faces into triangles
 494     newBoundaryFaces_.clear();
 495     newBoundaryOwners_.clear();
 496     newBoundaryPatches_.clear();
 497     face triF(3);
 498
 499     label nDecomposedFaces(0);
 500     forAll(bFaces, bfI)
 501     {
 502         const face& bf = bFaces[bfI];
 503
 504         bool store(true);
 505
 506         forAll(bf, i)
 507         {
 508             if(
 509                 (nodeType[bp[bf[i]]] == direction(2)) &&
 510                 featureEdge[faceEdges(bfI, i)] &&
 511                 featureEdge[faceEdges(bfI, bf.rcIndex(i))]
 512             )
 513             {
 514                 if( bf.size() > 4 )
 515                 {
 516                     store = false;
 517                     ++nDecomposedFaces;
 518                     decomposeCell_[boundaryFaceOwners[bfI]] = true;
 519                     decompose_ = true;
 520
 521                     //- decompose into triangles
 522                     const point p = bf.centre(mesh_.points());
 523                     triF[2] = mesh_.points().size();
 524                     mesh_.points().append(p);
 525
 526                     forAll(bf, j)
 527                     {
 528                         triF[0] = bf[j];
 529                         triF[1] = bf.nextLabel(j);
 530
 531                         newBoundaryFaces_.appendList(triF);
 532                         newBoundaryOwners_.append(boundaryFaceOwners[bfI]);
 533                         newBoundaryPatches_.append(facePatches[bfI]);
 534                     }
 535
 536                     break;
 537                 }
 538                 else if( bf.size() == 4 )
 539                 {
 540                     store = false;
 541                     ++nDecomposedFaces;
 542                     decomposeCell_[boundaryFaceOwners[bfI]] = true;
 543                     decompose_ = true;
 544
 545                     //- decompose the quad into 2 triangles
 546                     triF[0] = bf[i];
 547
 548                     triF[1] = bf.nextLabel(i);
 549                     triF[2] = bf[(i+2)%4];
 550                     newBoundaryFaces_.appendList(triF);
 551                     newBoundaryOwners_.append(boundaryFaceOwners[bfI]);
 552                     newBoundaryPatches_.append(facePatches[bfI]);
 553
 554                     triF[1] = bf[(i+2)%4];
 555                     triF[2] = bf.prevLabel(i);
 556                     newBoundaryFaces_.appendList(triF);
 557                     newBoundaryOwners_.append(boundaryFaceOwners[bfI]);
 558                     newBoundaryPatches_.append(facePatches[bfI]);
 559
 560                     break;
 561                 }
 562             }
 563         }
 564
 565         if( store )
 566         {
 567             //- face has not been altered
 568             newBoundaryFaces_.appendList(bf);
 569             newBoundaryOwners_.append(boundaryFaceOwners[bfI]);
 570             newBoundaryPatches_.append(facePatches[bfI]);
 571         }
 572     }
 573
 574     reduce(decompose_, maxOp<bool>());
 575
 576     if( returnReduce(nDecomposedFaces, sumOp<label>()) != 0 )
 577     {
 578         replaceBoundary();
 579         clearMeshSurface();
 580         mesh_.clearAddressingData();
 581     }
 582
 583     Info << "Finished with patch correction" << endl;
 584 }
 585
 586 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 587
 588 } // End namespace Foam
 589
 590 // ************************************************************************* //