src/dynamicMesh/boundaryMesh/boundaryMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "boundaryMesh.H"
  28 #include "Time.H"
  29 #include "polyMesh.H"
  30 #include "repatchPolyTopoChanger.H"
  31 #include "faceList.H"
  32 #include "octree.H"
  33 #include "octreeDataFaceList.H"
  34 #include "triSurface.H"
  35 #include "SortableList.H"
  36 #include "OFstream.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 defineTypeNameAndDebug(Foam::boundaryMesh, 0);
  41
  42 // Normal along which to divide faces into categories (used in getNearest)
  43 const Foam::vector Foam::boundaryMesh::splitNormal_(3, 2, 1);
  44
  45 // Distance to face tolerance for getNearest
  46 const Foam::scalar Foam::boundaryMesh::distanceTol_ = 1E-2;
  47
  48
  49 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  50
  51 // Returns number of feature edges connected to pointI
  52 Foam::label Foam::boundaryMesh::nFeatureEdges(label pointI) const
  53 {
  54     label nFeats = 0;
  55
  56     const labelList& pEdges = mesh().pointEdges()[pointI];
  57
  58     forAll(pEdges, pEdgeI)
  59     {
  60         label edgeI = pEdges[pEdgeI];
  61
  62         if (edgeToFeature_[edgeI] != -1)
  63         {
  64             nFeats++;
  65         }
  66     }
  67     return nFeats;
  68 }
  69
  70
  71 // Returns next feature edge connected to pointI
  72 Foam::label Foam::boundaryMesh::nextFeatureEdge
  73 (
  74     const label edgeI,
  75     const label vertI
  76 ) const
  77 {
  78     const labelList& pEdges = mesh().pointEdges()[vertI];
  79
  80     forAll(pEdges, pEdgeI)
  81     {
  82         label nbrEdgeI = pEdges[pEdgeI];
  83
  84         if (nbrEdgeI != edgeI)
  85         {
  86             label featI = edgeToFeature_[nbrEdgeI];
  87
  88             if (featI != -1)
  89             {
  90                 return nbrEdgeI;
  91             }
  92         }
  93     }
  94
  95     return -1;
  96 }
  97
  98
  99 // Finds connected feature edges, starting from startPointI and returns
 100 // feature labels (not edge labels). Marks feature edges handled in
 101 // featVisited.
 102 Foam::labelList Foam::boundaryMesh::collectSegment
 103 (
 104     const boolList& isFeaturePoint,
 105     const label startEdgeI,
 106     boolList& featVisited
 107 ) const
 108 {
 109     // Find starting feature point on edge.
 110
 111     label edgeI = startEdgeI;
 112
 113     const edge& e = mesh().edges()[edgeI];
 114
 115     label vertI = e.start();
 116
 117     while (!isFeaturePoint[vertI])
 118     {
 119         // Step to next feature edge
 120
 121         edgeI = nextFeatureEdge(edgeI, vertI);
 122
 123         if ((edgeI == -1) || (edgeI == startEdgeI))
 124         {
 125             break;
 126         }
 127
 128         // Step to next vertex on edge
 129
 130         const edge& e = mesh().edges()[edgeI];
 131
 132         vertI = e.otherVertex(vertI);
 133     }
 134
 135     //
 136     // Now we have:
 137     //    edgeI : first edge on this segment
 138     //    vertI : one of the endpoints of this segment
 139     //
 140     // Start walking other way and storing edges as we go along.
 141     //
 142
 143     // Untrimmed storage for current segment
 144     labelList featLabels(featureEdges_.size());
 145
 146     label featLabelI = 0;
 147
 148     label initEdgeI = edgeI;
 149
 150     do
 151     {
 152         // Mark edge as visited
 153         label featI = edgeToFeature_[edgeI];
 154
 155         if (featI == -1)
 156         {
 157             FatalErrorIn("boundaryMesh::collectSegment")
 158                 << "Problem" << abort(FatalError);
 159         }
 160         featLabels[featLabelI++] = featI;
 161
 162         featVisited[featI] = true;
 163
 164         // Step to next vertex on edge
 165
 166         const edge& e = mesh().edges()[edgeI];
 167
 168         vertI = e.otherVertex(vertI);
 169
 170         // Step to next feature edge
 171
 172         edgeI = nextFeatureEdge(edgeI, vertI);
 173
 174         if ((edgeI == -1) || (edgeI == initEdgeI))
 175         {
 176             break;
 177         }
 178     }
 179     while (!isFeaturePoint[vertI]);
 180
 181
 182     // Trim to size
 183     featLabels.setSize(featLabelI);
 184
 185     return featLabels;
 186 }
 187
 188
 189 void Foam::boundaryMesh::markEdges
 190 (
 191     const label maxDistance,
 192     const label edgeI,
 193     const label distance,
 194     labelList& minDistance,
 195     DynamicList<label>& visited
 196 ) const
 197 {
 198     if (distance < maxDistance)
 199     {
 200         // Don't do anything if reached beyond maxDistance.
 201
 202         if (minDistance[edgeI] == -1)
 203         {
 204             // First visit of edge. Store edge label.
 205             visited.append(edgeI);
 206         }
 207         else if (minDistance[edgeI] <= distance)
 208         {
 209             // Already done this edge
 210             return;
 211         }
 212
 213         minDistance[edgeI] = distance;
 214
 215         const edge& e = mesh().edges()[edgeI];
 216
 217         // Do edges connected to e.start
 218         const labelList& startEdges = mesh().pointEdges()[e.start()];
 219
 220         forAll(startEdges, pEdgeI)
 221         {
 222             markEdges
 223             (
 224                 maxDistance,
 225                 startEdges[pEdgeI],
 226                 distance+1,
 227                 minDistance,
 228                 visited
 229             );
 230         }
 231
 232         // Do edges connected to e.end
 233         const labelList& endEdges = mesh().pointEdges()[e.end()];
 234
 235         forAll(endEdges, pEdgeI)
 236         {
 237             markEdges
 238             (
 239                 maxDistance,
 240                 endEdges[pEdgeI],
 241                 distance+1,
 242                 minDistance,
 243                 visited
 244             );
 245         }
 246     }
 247 }
 248
 249
 250 Foam::label Foam::boundaryMesh::findPatchID
 251 (
 252     const polyPatchList& patches,
 253     const word& patchName
 254 ) const
 255 {
 256     forAll(patches, patchI)
 257     {
 258         if (patches[patchI].name() == patchName)
 259         {
 260             return patchI;
 261         }
 262     }
 263
 264     return -1;
 265 }
 266
 267
 268 Foam::wordList Foam::boundaryMesh::patchNames() const
 269 {
 270     wordList names(patches_.size());
 271
 272     forAll(patches_, patchI)
 273     {
 274         names[patchI] = patches_[patchI].name();
 275     }
 276     return names;
 277 }
 278
 279
 280 Foam::label Foam::boundaryMesh::whichPatch
 281 (
 282     const polyPatchList& patches,
 283     const label faceI
 284 ) const
 285 {
 286     forAll(patches, patchI)
 287     {
 288         const polyPatch& pp = patches[patchI];
 289
 290         if ((faceI >= pp.start()) && (faceI < (pp.start() + pp.size())))
 291         {
 292             return patchI;
 293         }
 294     }
 295     return -1;
 296 }
 297
 298
 299 // Gets labels of changed faces and propagates them to the edges. Returns
 300 // labels of edges changed.
 301 Foam::labelList Foam::boundaryMesh::faceToEdge
 302 (
 303     const boolList& regionEdge,
 304     const label region,
 305     const labelList& changedFaces,
 306     labelList& edgeRegion
 307 ) const
 308 {
 309     labelList changedEdges(mesh().nEdges(), -1);
 310     label changedI = 0;
 311
 312     forAll(changedFaces, i)
 313     {
 314         label faceI = changedFaces[i];
 315
 316         const labelList& fEdges = mesh().faceEdges()[faceI];
 317
 318         forAll(fEdges, fEdgeI)
 319         {
 320             label edgeI = fEdges[fEdgeI];
 321
 322             if (!regionEdge[edgeI] && (edgeRegion[edgeI] == -1))
 323             {
 324                 edgeRegion[edgeI] = region;
 325
 326                 changedEdges[changedI++] = edgeI;
 327             }
 328         }
 329     }
 330
 331     changedEdges.setSize(changedI);
 332
 333     return changedEdges;
 334 }
 335
 336
 337 // Reverse of faceToEdge: gets edges and returns faces
 338 Foam::labelList Foam::boundaryMesh::edgeToFace
 339 (
 340     const label region,
 341     const labelList& changedEdges,
 342     labelList& faceRegion
 343 ) const
 344 {
 345     labelList changedFaces(mesh().size(), -1);
 346     label changedI = 0;
 347
 348     forAll(changedEdges, i)
 349     {
 350         label edgeI = changedEdges[i];
 351
 352         const labelList& eFaces = mesh().edgeFaces()[edgeI];
 353
 354         forAll(eFaces, eFaceI)
 355         {
 356             label faceI = eFaces[eFaceI];
 357
 358             if (faceRegion[faceI] == -1)
 359             {
 360                 faceRegion[faceI] = region;
 361
 362                 changedFaces[changedI++] = faceI;
 363             }
 364         }
 365     }
 366
 367     changedFaces.setSize(changedI);
 368
 369     return changedFaces;
 370 }
 371
 372
 373 // Finds area, starting at faceI, delimited by borderEdge
 374 void Foam::boundaryMesh::markZone
 375 (
 376     const boolList& borderEdge,
 377     label faceI,
 378     label currentZone,
 379     labelList& faceZone
 380 ) const
 381 {
 382     faceZone[faceI] = currentZone;
 383
 384     // List of faces whose faceZone has been set.
 385     labelList changedFaces(1, faceI);
 386     // List of edges whose faceZone has been set.
 387     labelList changedEdges;
 388
 389     // Zones on all edges.
 390     labelList edgeZone(mesh().nEdges(), -1);
 391
 392     while(1)
 393     {
 394         changedEdges = faceToEdge
 395         (
 396             borderEdge,
 397             currentZone,
 398             changedFaces,
 399             edgeZone
 400         );
 401
 402         if (debug)
 403         {
 404             Pout<< "From changedFaces:" << changedFaces.size()
 405                 << " to changedEdges:" << changedEdges.size()
 406                 << endl;
 407         }
 408
 409         if (changedEdges.empty())
 410         {
 411             break;
 412         }
 413
 414         changedFaces = edgeToFace(currentZone, changedEdges, faceZone);
 415
 416         if (debug)
 417         {
 418             Pout<< "From changedEdges:" << changedEdges.size()
 419                 << " to changedFaces:" << changedFaces.size()
 420                 << endl;
 421         }
 422
 423         if (changedFaces.empty())
 424         {
 425             break;
 426         }
 427     }
 428 }
 429
 430
 431 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 432
 433 // Null constructor
 434 Foam::boundaryMesh::boundaryMesh()
 435 :
 436     meshPtr_(NULL),
 437     patches_(),
 438     meshFace_(),
 439     featurePoints_(),
 440     featureEdges_(),
 441     featureToEdge_(),
 442     edgeToFeature_(),
 443     featureSegments_(),
 444     extraEdges_()
 445 {}
 446
 447
 448 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 449
 450 Foam::boundaryMesh::~boundaryMesh()
 451 {
 452     clearOut();
 453 }
 454
 455
 456 void Foam::boundaryMesh::clearOut()
 457 {
 458     if (meshPtr_)
 459     {
 460         delete meshPtr_;
 461
 462         meshPtr_ = NULL;
 463     }
 464 }
 465
 466
 467 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 468
 469 void Foam::boundaryMesh::read(const polyMesh& mesh)
 470 {
 471     patches_.clear();
 472
 473     patches_.setSize(mesh.boundaryMesh().size());
 474
 475     // Number of boundary faces
 476     label nBFaces = mesh.nFaces() - mesh.nInternalFaces();
 477
 478     faceList bFaces(nBFaces);
 479
 480     meshFace_.setSize(nBFaces);
 481
 482     label bFaceI = 0;
 483
 484     // Collect all boundary faces.
 485     forAll(mesh.boundaryMesh(), patchI)
 486     {
 487         const polyPatch& pp = mesh.boundaryMesh()[patchI];
 488
 489         patches_.set
 490         (
 491             patchI,
 492             new boundaryPatch
 493             (
 494                 pp.name(),
 495                 patchI,
 496                 pp.size(),
 497                 bFaceI,
 498                 pp.type()
 499             )
 500         );
 501
 502         // Collect all faces in global numbering.
 503         forAll(pp, patchFaceI)
 504         {
 505             meshFace_[bFaceI] = pp.start() + patchFaceI;
 506
 507             bFaces[bFaceI] = pp[patchFaceI];
 508
 509             bFaceI++;
 510         }
 511     }
 512
 513
 514     if (debug)
 515     {
 516         Pout<< "read : patches now:" << endl;
 517
 518         forAll(patches_, patchI)
 519         {
 520             const boundaryPatch& bp = patches_[patchI];
 521
 522             Pout<< "    name  : " << bp.name() << endl
 523                 << "    size  : " << bp.size() << endl
 524                 << "    start : " << bp.start() << endl
 525                 << "    type  : " << bp.physicalType() << endl
 526                 << endl;
 527         }
 528     }
 529
 530     //
 531     // Construct single patch for all of boundary
 532     //
 533
 534     // Temporary primitivePatch to calculate compact points & faces.
 535     PrimitivePatch<face, List, const pointField&> globalPatch
 536     (
 537         bFaces,
 538         mesh.points()
 539     );
 540
 541     // Store in local(compact) addressing
 542     clearOut();
 543
 544     meshPtr_ = new bMesh(globalPatch.localFaces(), globalPatch.localPoints());
 545
 546
 547     if (debug & 2)
 548     {
 549         const bMesh& msh = *meshPtr_;
 550
 551         Pout<< "** Start of Faces **" << endl;
 552
 553         forAll(msh, faceI)
 554         {
 555             const face& f = msh[faceI];
 556
 557             point ctr(vector::zero);
 558
 559             forAll(f, fp)
 560             {
 561                 ctr += msh.points()[f[fp]];
 562             }
 563             ctr /= f.size();
 564
 565             Pout<< "    " << faceI
 566                 << " ctr:" << ctr
 567                 << " verts:" << f
 568                 << endl;
 569         }
 570
 571         Pout<< "** End of Faces **" << endl;
 572
 573         Pout<< "** Start of Points **" << endl;
 574
 575         forAll(msh.points(), pointI)
 576         {
 577             Pout<< "    " << pointI
 578                 << " coord:" << msh.points()[pointI]
 579                 << endl;
 580         }
 581
 582         Pout<< "** End of Points **" << endl;
 583     }
 584
 585     // Clear edge storage
 586     featurePoints_.setSize(0);
 587     featureEdges_.setSize(0);
 588
 589     featureToEdge_.setSize(0);
 590     edgeToFeature_.setSize(meshPtr_->nEdges());
 591     edgeToFeature_ = -1;
 592
 593     featureSegments_.setSize(0);
 594
 595     extraEdges_.setSize(0);
 596 }
 597
 598
 599 void Foam::boundaryMesh::readTriSurface(const fileName& fName)
 600 {
 601     triSurface surf(fName);
 602
 603     if (surf.empty())
 604     {
 605         return;
 606     }
 607
 608     // Sort according to region
 609     SortableList<label> regions(surf.size());
 610
 611     forAll(surf, triI)
 612     {
 613         regions[triI] = surf[triI].region();
 614     }
 615     regions.sort();
 616
 617     // Determine region mapping.
 618     Map<label> regionToBoundaryPatch;
 619
 620     label oldRegion = -1111;
 621     label boundPatch = 0;
 622
 623     forAll(regions, i)
 624     {
 625         if (regions[i] != oldRegion)
 626         {
 627             regionToBoundaryPatch.insert(regions[i], boundPatch);
 628
 629             oldRegion = regions[i];
 630             boundPatch++;
 631         }
 632     }
 633
 634     const geometricSurfacePatchList& surfPatches = surf.patches();
 635
 636     patches_.clear();
 637
 638     if (surfPatches.size() == regionToBoundaryPatch.size())
 639     {
 640         // There are as many surface patches as region numbers in triangles
 641         // so use the surface patches
 642
 643         patches_.setSize(surfPatches.size());
 644
 645         // Take over patches, setting size to 0 for now.
 646         forAll(surfPatches, patchI)
 647         {
 648             const geometricSurfacePatch& surfPatch = surfPatches[patchI];
 649
 650             patches_.set
 651             (
 652                 patchI,
 653                 new boundaryPatch
 654                 (
 655                     surfPatch.name(),
 656                     patchI,
 657                     0,
 658                     0,
 659                     surfPatch.geometricType()
 660                 )
 661             );
 662         }
 663     }
 664     else
 665     {
 666         // There are not enough surface patches. Make up my own.
 667
 668         patches_.setSize(regionToBoundaryPatch.size());
 669
 670         forAll(patches_, patchI)
 671         {
 672             patches_.set
 673             (
 674                 patchI,
 675                 new boundaryPatch
 676                 (
 677                     "patch" + name(patchI),
 678                     patchI,
 679                     0,
 680                     0,
 681                     "empty"
 682                 )
 683             );
 684         }
 685     }
 686
 687     //
 688     // Copy according into bFaces according to regions
 689     //
 690
 691     const labelList& indices = regions.indices();
 692
 693     faceList bFaces(surf.size());
 694
 695     meshFace_.setSize(surf.size());
 696
 697     label bFaceI = 0;
 698
 699     // Current region number
 700     label surfRegion = regions[0];
 701     label foamRegion = regionToBoundaryPatch[surfRegion];
 702
 703     Pout<< "Surface region " << surfRegion << " becomes boundary patch "
 704         << foamRegion << " with name " << patches_[foamRegion].name() << endl;
 705
 706
 707     // Index in bFaces of start of current patch
 708     label startFaceI = 0;
 709
 710     forAll(indices, indexI)
 711     {
 712         label triI = indices[indexI];
 713
 714         const labelledTri& tri = surf.localFaces()[triI];
 715
 716         if (tri.region() != surfRegion)
 717         {
 718             // Change of region. We now know the size of the previous one.
 719             boundaryPatch& bp = patches_[foamRegion];
 720
 721             bp.size() = bFaceI - startFaceI;
 722             bp.start() = startFaceI;
 723
 724             surfRegion = tri.region();
 725             foamRegion = regionToBoundaryPatch[surfRegion];
 726
 727             Pout<< "Surface region " << surfRegion << " becomes boundary patch "
 728                 << foamRegion << " with name " << patches_[foamRegion].name()
 729                 << endl;
 730
 731             startFaceI = bFaceI;
 732         }
 733
 734         meshFace_[bFaceI] = triI;
 735
 736         bFaces[bFaceI++] = face(tri);
 737     }
 738
 739     // Final region
 740     boundaryPatch& bp = patches_[foamRegion];
 741
 742     bp.size() = bFaceI - startFaceI;
 743     bp.start() = startFaceI;
 744
 745     //
 746     // Construct single primitivePatch for all of boundary
 747     //
 748
 749     clearOut();
 750
 751     // Store compact.
 752     meshPtr_ = new bMesh(bFaces, surf.localPoints());
 753
 754     // Clear edge storage
 755     featurePoints_.setSize(0);
 756     featureEdges_.setSize(0);
 757
 758     featureToEdge_.setSize(0);
 759     edgeToFeature_.setSize(meshPtr_->nEdges());
 760     edgeToFeature_ = -1;
 761
 762     featureSegments_.setSize(0);
 763
 764     extraEdges_.setSize(0);
 765 }
 766
 767
 768 void Foam::boundaryMesh::writeTriSurface(const fileName& fName) const
 769 {
 770     geometricSurfacePatchList surfPatches(patches_.size());
 771
 772     forAll(patches_, patchI)
 773     {
 774         const boundaryPatch& bp = patches_[patchI];
 775
 776         surfPatches[patchI] =
 777             geometricSurfacePatch
 778             (
 779                 bp.physicalType(),
 780                 bp.name(),
 781                 patchI
 782             );
 783     }
 784
 785     //
 786     // Simple triangulation.
 787     //
 788
 789     // Get number of triangles per face
 790     labelList nTris(mesh().size());
 791
 792     label totalNTris = getNTris(0, mesh().size(), nTris);
 793
 794     // Determine per face the starting triangle.
 795     labelList startTri(mesh().size());
 796
 797     label triI = 0;
 798
 799     forAll(mesh(), faceI)
 800     {
 801         startTri[faceI] = triI;
 802
 803         triI += nTris[faceI];
 804     }
 805
 806     // Triangulate
 807     labelList triVerts(3*totalNTris);
 808
 809     triangulate(0, mesh().size(), totalNTris, triVerts);
 810
 811     // Convert to labelledTri
 812
 813     List<labelledTri> tris(totalNTris);
 814
 815     triI = 0;
 816
 817     forAll(patches_, patchI)
 818     {
 819         const boundaryPatch& bp = patches_[patchI];
 820
 821         forAll(bp, patchFaceI)
 822         {
 823             label faceI = bp.start() + patchFaceI;
 824
 825             label triVertI = 3*startTri[faceI];
 826
 827             for (label faceTriI = 0; faceTriI < nTris[faceI]; faceTriI++)
 828             {
 829                 label v0 = triVerts[triVertI++];
 830                 label v1 = triVerts[triVertI++];
 831                 label v2 = triVerts[triVertI++];
 832
 833                 tris[triI++] = labelledTri(v0, v1, v2, patchI);
 834             }
 835         }
 836     }
 837
 838     triSurface surf(tris, surfPatches, mesh().points());
 839
 840     OFstream surfStream(fName);
 841
 842     surf.write(surfStream);
 843 }
 844
 845
 846 // Get index in this (boundaryMesh) of face nearest to each boundary face in
 847 // pMesh.
 848 // Origininally all triangles/faces of boundaryMesh would be bunged into
 849 // one big octree. Problem was that faces on top of each other, differing
 850 // only in sign of normal, could not be found separately. It would always
 851 // find only one. We could detect that it was probably finding the wrong one
 852 // (based on normal) but could not 'tell' the octree to retrieve the other
 853 // one (since they occupy exactly the same space)
 854 // So now faces get put into different octrees depending on normal.
 855 // !It still will not be possible to differentiate between two faces on top
 856 // of each other having the same normal
 857 Foam::labelList Foam::boundaryMesh::getNearest
 858 (
 859     const primitiveMesh& pMesh,
 860     const vector& searchSpan
 861 ) const
 862 {
 863
 864     // Divide faces into two bins acc. to normal
 865     // - left of splitNormal
 866     // - right ,,
 867     DynamicList<label> leftFaces(mesh().size()/2);
 868     DynamicList<label> rightFaces(mesh().size()/2);
 869
 870     forAll(mesh(), bFaceI)
 871     {
 872         scalar sign = mesh().faceNormals()[bFaceI] & splitNormal_;
 873
 874         if (sign > -1E-5)
 875         {
 876             rightFaces.append(bFaceI);
 877         }
 878         if (sign < 1E-5)
 879         {
 880             leftFaces.append(bFaceI);
 881         }
 882     }
 883
 884     leftFaces.shrink();
 885     rightFaces.shrink();
 886
 887     if (debug)
 888     {
 889         Pout<< "getNearest :"
 890             << " rightBin:" << rightFaces.size()
 891             << " leftBin:" << leftFaces.size()
 892             << endl;
 893     }
 894
 895
 896     // Overall bb
 897     treeBoundBox overallBb(mesh().localPoints());
 898
 899     // Extend domain slightly (also makes it 3D if was 2D)
 900     // Note asymmetry to avoid having faces align with octree cubes.
 901     scalar tol = 1E-6 * overallBb.avgDim();
 902
 903     point& bbMin = overallBb.min();
 904     bbMin.x() -= tol;
 905     bbMin.y() -= tol;
 906     bbMin.z() -= tol;
 907
 908     point& bbMax = overallBb.max();
 909     bbMax.x() += 2*tol;
 910     bbMax.y() += 2*tol;
 911     bbMax.z() += 2*tol;
 912
 913     // Create the octrees
 914     octree<octreeDataFaceList> leftTree
 915     (
 916         overallBb,
 917         octreeDataFaceList
 918         (
 919             mesh(),
 920             leftFaces
 921         ),
 922         1,
 923         20,
 924         10
 925     );
 926     octree<octreeDataFaceList> rightTree
 927     (
 928         overallBb,
 929         octreeDataFaceList
 930         (
 931             mesh(),
 932             rightFaces
 933         ),
 934         1,
 935         20,
 936         10
 937     );
 938
 939     if (debug)
 940     {
 941         Pout<< "getNearest : built trees" << endl;
 942     }
 943
 944
 945     const vectorField& ns = mesh().faceNormals();
 946
 947
 948     //
 949     // Search nearest triangle centre for every polyMesh boundary face
 950     //
 951
 952     labelList nearestBFaceI(pMesh.nFaces() - pMesh.nInternalFaces());
 953
 954     treeBoundBox tightest;
 955
 956     const scalar searchDim = mag(searchSpan);
 957
 958     forAll(nearestBFaceI, patchFaceI)
 959     {
 960         label meshFaceI = pMesh.nInternalFaces() + patchFaceI;
 961
 962         const point& ctr = pMesh.faceCentres()[meshFaceI];
 963
 964         if (debug && (patchFaceI % 1000) == 0)
 965         {
 966             Pout<< "getNearest : patchFace:" << patchFaceI
 967                 << " meshFaceI:" << meshFaceI << " ctr:" << ctr << endl;
 968         }
 969
 970
 971         // Get normal from area vector
 972         vector n = pMesh.faceAreas()[meshFaceI];
 973         scalar area = mag(n);
 974         n /= area;
 975
 976         scalar typDim = -GREAT;
 977         const face& f = pMesh.faces()[meshFaceI];
 978
 979         forAll(f, fp)
 980         {
 981             typDim = max(typDim, mag(pMesh.points()[f[fp]] - ctr));
 982         }
 983
 984         // Search right tree
 985         tightest.min() = ctr - searchSpan;
 986         tightest.max() = ctr + searchSpan;
 987         scalar rightDist = searchDim;
 988         label rightI = rightTree.findNearest(ctr, tightest, rightDist);
 989
 990
 991         // Search left tree. Note: could start from rightDist bounding box
 992         // instead of starting from top.
 993         tightest.min() = ctr - searchSpan;
 994         tightest.max() = ctr + searchSpan;
 995         scalar leftDist = searchDim;
 996         label leftI = leftTree.findNearest(ctr, tightest, leftDist);
 997
 998
 999         if (rightI == -1)
1000         {
1001             // No face found in right tree.
1002
1003             if (leftI == -1)
1004             {
1005                 // No face found in left tree.
1006                 nearestBFaceI[patchFaceI] = -1;
1007             }
1008             else
1009             {
1010                 // Found in left but not in right. Choose left regardless if
1011                 // correct sign. Note: do we want this?
1012                 nearestBFaceI[patchFaceI] = leftFaces[leftI];
1013             }
1014         }
1015         else
1016         {
1017             if (leftI == -1)
1018             {
1019                 // Found in right but not in left. Choose right regardless if
1020                 // correct sign. Note: do we want this?
1021                 nearestBFaceI[patchFaceI] = rightFaces[rightI];
1022             }
1023             else
1024             {
1025                 // Found in both trees. Compare normals.
1026
1027                 scalar rightSign = n & ns[rightFaces[rightI]];
1028                 scalar leftSign = n & ns[leftFaces[leftI]];
1029
1030                 if
1031                 (
1032                     (rightSign > 0 && leftSign > 0)
1033                  || (rightSign < 0 && leftSign < 0)
1034                 )
1035                 {
1036                     // Both same sign. Choose nearest.
1037                     if (rightDist < leftDist)
1038                     {
1039                         nearestBFaceI[patchFaceI] = rightFaces[rightI];
1040                     }
1041                     else
1042                     {
1043                         nearestBFaceI[patchFaceI] = leftFaces[leftI];
1044                     }
1045                 }
1046                 else
1047                 {
1048                     // Differing sign.
1049                     // - if both near enough choose one with correct sign
1050                     // - otherwise choose nearest.
1051
1052                     // Get local dimension as max of distance between ctr and
1053                     // any face vertex.
1054
1055                     typDim *= distanceTol_;
1056
1057                     if (rightDist < typDim && leftDist < typDim)
1058                     {
1059                         // Different sign and nearby. Choosing matching normal
1060                         if (rightSign > 0)
1061                         {
1062                             nearestBFaceI[patchFaceI] = rightFaces[rightI];
1063                         }
1064                         else
1065                         {
1066                             nearestBFaceI[patchFaceI] = leftFaces[leftI];
1067                         }
1068                     }
1069                     else
1070                     {
1071                         // Different sign but faraway. Choosing nearest.
1072                         if (rightDist < leftDist)
1073                         {
1074                             nearestBFaceI[patchFaceI] = rightFaces[rightI];
1075                         }
1076                         else
1077                         {
1078                             nearestBFaceI[patchFaceI] = leftFaces[leftI];
1079                         }
1080                     }
1081                 }
1082             }
1083         }
1084     }
1085
1086     return nearestBFaceI;
1087 }
1088
1089
1090 void Foam::boundaryMesh::patchify
1091 (
1092     const labelList& nearest,
1093     const polyBoundaryMesh& oldPatches,
1094     polyMesh& newMesh
1095 ) const
1096 {
1097
1098     // 2 cases to be handled:
1099     // A- patches in boundaryMesh patches_
1100     // B- patches not in boundaryMesh patches_ but in polyMesh
1101
1102     // Create maps from patch name to new patch index.
1103     HashTable<label> nameToIndex(2*patches_.size());
1104
1105     Map<word> indexToName(2*patches_.size());
1106
1107
1108     label nNewPatches = patches_.size();
1109
1110     forAll(oldPatches, oldPatchI)
1111     {
1112         const polyPatch& patch = oldPatches[oldPatchI];
1113
1114         label newPatchI = findPatchID(patch.name());
1115
1116         if (newPatchI != -1)
1117         {
1118             nameToIndex.insert(patch.name(), newPatchI);
1119
1120             indexToName.insert(newPatchI, patch.name());
1121         }
1122     }
1123
1124     // Include all boundaryPatches not yet in nameToIndex (i.e. not in old
1125     // patches)
1126     forAll(patches_, bPatchI)
1127     {
1128         const boundaryPatch& bp = patches_[bPatchI];
1129
1130         if (!nameToIndex.found(bp.name()))
1131         {
1132             nameToIndex.insert(bp.name(), bPatchI);
1133
1134             indexToName.insert(bPatchI, bp.name());
1135         }
1136     }
1137
1138     // Pass1:
1139     // Copy names&type of patches (with zero size) from old mesh as far as
1140     // possible. First patch created gets all boundary faces; all others get
1141     // zero faces (repatched later on). Exception is coupled patches which
1142     // keep their size.
1143
1144     List<polyPatch*> newPatchPtrList(nNewPatches);
1145
1146     label meshFaceI = newMesh.nInternalFaces();
1147
1148     // First patch gets all non-coupled faces
1149     label facesToBeDone = newMesh.nFaces() - newMesh.nInternalFaces();
1150
1151     forAll(patches_, bPatchI)
1152     {
1153         const boundaryPatch& bp = patches_[bPatchI];
1154
1155         label newPatchI = nameToIndex[bp.name()];
1156
1157         // Find corresponding patch in polyMesh
1158         label oldPatchI = findPatchID(oldPatches, bp.name());
1159
1160         if (oldPatchI == -1)
1161         {
1162             // Newly created patch. Gets all or zero faces.
1163             if (debug)
1164             {
1165                 Pout<< "patchify : Creating new polyPatch:" << bp.name()
1166                     << " type:" << bp.physicalType() << endl;
1167             }
1168
1169             newPatchPtrList[newPatchI] = polyPatch::New
1170             (
1171                 bp.physicalType(),
1172                 bp.name(),
1173                 facesToBeDone,
1174                 meshFaceI,
1175                 newPatchI,
1176                 newMesh.boundaryMesh()
1177             ).ptr();
1178
1179             meshFaceI += facesToBeDone;
1180
1181             // first patch gets all boundary faces; all others get 0.
1182             facesToBeDone = 0;
1183         }
1184         else
1185         {
1186             // Existing patch. Gets all or zero faces.
1187             const polyPatch& oldPatch = oldPatches[oldPatchI];
1188
1189             if (debug)
1190             {
1191                 Pout<< "patchify : Cloning existing polyPatch:"
1192                     << oldPatch.name() << endl;
1193             }
1194
1195             newPatchPtrList[newPatchI] = oldPatch.clone
1196             (
1197                 newMesh.boundaryMesh(),
1198                 newPatchI,
1199                 facesToBeDone,
1200                 meshFaceI
1201             ).ptr();
1202
1203             meshFaceI += facesToBeDone;
1204
1205             // first patch gets all boundary faces; all others get 0.
1206             facesToBeDone = 0;
1207         }
1208     }
1209
1210
1211     if (debug)
1212     {
1213         Pout<< "Patchify : new polyPatch list:" << endl;
1214
1215         forAll(newPatchPtrList, patchI)
1216         {
1217             const polyPatch& newPatch = *newPatchPtrList[patchI];
1218
1219             if (debug)
1220             {
1221                 Pout<< "polyPatch:" << newPatch.name() << endl
1222                     << "    type :" << newPatch.typeName << endl
1223                     << "    size :" << newPatch.size() << endl
1224                     << "    start:" << newPatch.start() << endl
1225                     << "    index:" << patchI << endl;
1226             }
1227         }
1228     }
1229
1230     // Actually add new list of patches
1231     repatchPolyTopoChanger polyMeshRepatcher(newMesh);
1232     polyMeshRepatcher.changePatches(newPatchPtrList);
1233
1234
1235     // Pass2:
1236     // Change patch type for face
1237
1238     if (newPatchPtrList.size())
1239     {
1240         List<DynamicList<label> > patchFaces(nNewPatches);
1241
1242         // Give reasonable estimate for size of patches
1243         label nAvgFaces =
1244             (newMesh.nFaces() - newMesh.nInternalFaces())
1245           / nNewPatches;
1246
1247         forAll(patchFaces, newPatchI)
1248         {
1249             patchFaces[newPatchI].setCapacity(nAvgFaces);
1250         }
1251
1252         //
1253         // Sort faces acc. to new patch index. Can loop over all old patches
1254         // since will contain all faces.
1255         //
1256
1257         forAll(oldPatches, oldPatchI)
1258         {
1259             const polyPatch& patch = oldPatches[oldPatchI];
1260
1261             forAll(patch, patchFaceI)
1262             {
1263                 // Put face into region given by nearest boundary face
1264
1265                 label meshFaceI = patch.start() + patchFaceI;
1266
1267                 label bFaceI = meshFaceI - newMesh.nInternalFaces();
1268
1269                 patchFaces[whichPatch(nearest[bFaceI])].append(meshFaceI);
1270             }
1271         }
1272
1273         forAll(patchFaces, newPatchI)
1274         {
1275             patchFaces[newPatchI].shrink();
1276         }
1277
1278
1279         // Change patch > 0. (since above we put all faces into the zeroth
1280         // patch)
1281
1282         for (label newPatchI = 1; newPatchI < patchFaces.size(); newPatchI++)
1283         {
1284             const labelList& pFaces = patchFaces[newPatchI];
1285
1286             forAll(pFaces, pFaceI)
1287             {
1288                 polyMeshRepatcher.changePatchID(pFaces[pFaceI], newPatchI);
1289             }
1290         }
1291
1292         polyMeshRepatcher.repatch();
1293     }
1294 }
1295
1296
1297 void Foam::boundaryMesh::setFeatureEdges(const scalar minCos)
1298 {
1299     edgeToFeature_.setSize(mesh().nEdges());
1300
1301     edgeToFeature_ = -1;
1302
1303     // 1. Mark feature edges
1304
1305     // Storage for edge labels that are features. Trim later.
1306     featureToEdge_.setSize(mesh().nEdges());
1307
1308     label featureI = 0;
1309
1310     if (minCos >= 0.9999)
1311     {
1312         // Select everything
1313         forAll(mesh().edges(), edgeI)
1314         {
1315             edgeToFeature_[edgeI] = featureI;
1316             featureToEdge_[featureI++] = edgeI;
1317         }
1318     }
1319     else
1320     {
1321         forAll(mesh().edges(), edgeI)
1322         {
1323             const labelList& eFaces = mesh().edgeFaces()[edgeI];
1324
1325             if (eFaces.size() == 2)
1326             {
1327                 label face0I = eFaces[0];
1328
1329                 label face1I = eFaces[1];
1330
1331                 ////- Uncomment below code if you want to include patch
1332                 ////  boundaries in feature edges.
1333                 //if (whichPatch(face0I) != whichPatch(face1I))
1334                 //{
1335                 //    edgeToFeature_[edgeI] = featureI;
1336                 //    featureToEdge_[featureI++] = edgeI;
1337                 //}
1338                 //else
1339                 {
1340                     const vector& n0 = mesh().faceNormals()[face0I];
1341
1342                     const vector& n1 = mesh().faceNormals()[face1I];
1343
1344                     float cosAng = n0 & n1;
1345
1346                     if (cosAng < minCos)
1347                     {
1348                         edgeToFeature_[edgeI] = featureI;
1349                         featureToEdge_[featureI++] = edgeI;
1350                     }
1351                 }
1352             }
1353             else
1354             {
1355                 //Should not occur: 0 or more than two faces
1356                 edgeToFeature_[edgeI] = featureI;
1357                 featureToEdge_[featureI++] = edgeI;
1358             }
1359         }
1360     }
1361
1362     // Trim featureToEdge_ to actual number of edges.
1363     featureToEdge_.setSize(featureI);
1364
1365     //
1366     // Compact edges i.e. relabel vertices.
1367     //
1368
1369     featureEdges_.setSize(featureI);
1370     featurePoints_.setSize(mesh().nPoints());
1371
1372     labelList featToMeshPoint(mesh().nPoints(), -1);
1373
1374     label featPtI = 0;
1375
1376     forAll(featureToEdge_, fEdgeI)
1377     {
1378         label edgeI = featureToEdge_[fEdgeI];
1379
1380         const edge& e = mesh().edges()[edgeI];
1381
1382         label start = featToMeshPoint[e.start()];
1383
1384         if (start == -1)
1385         {
1386             featToMeshPoint[e.start()] = featPtI;
1387
1388             featurePoints_[featPtI] = mesh().points()[e.start()];
1389
1390             start = featPtI;
1391
1392             featPtI++;
1393         }
1394
1395         label end = featToMeshPoint[e.end()];
1396
1397         if (end == -1)
1398         {
1399             featToMeshPoint[e.end()] = featPtI;
1400
1401             featurePoints_[featPtI] = mesh().points()[e.end()];
1402
1403             end = featPtI;
1404
1405             featPtI++;
1406         }
1407
1408         // Store with renumbered vertices.
1409         featureEdges_[fEdgeI] = edge(start, end);
1410     }
1411
1412     // Compact points
1413     featurePoints_.setSize(featPtI);
1414
1415
1416     //
1417     // 2. Mark endpoints of feature segments. These are points with
1418     // != 2 feature edges connected.
1419     // Note: can add geometric constraint here as well that if 2 feature
1420     // edges the angle between them should be less than xxx.
1421     //
1422
1423     boolList isFeaturePoint(mesh().nPoints(), false);
1424
1425     forAll(featureToEdge_, featI)
1426     {
1427         label edgeI = featureToEdge_[featI];
1428
1429         const edge& e = mesh().edges()[edgeI];
1430
1431         if (nFeatureEdges(e.start()) != 2)
1432         {
1433             isFeaturePoint[e.start()] = true;
1434         }
1435
1436         if (nFeatureEdges(e.end()) != 2)
1437         {
1438             isFeaturePoint[e.end()] = true;
1439         }
1440     }
1441
1442
1443     //
1444     // 3: Split feature edges into segments:
1445     // find point with not 2 feature edges -> start of feature segment
1446     //
1447
1448     DynamicList<labelList> segments;
1449
1450
1451     boolList featVisited(featureToEdge_.size(), false);
1452
1453     do
1454     {
1455         label startFeatI = -1;
1456
1457         forAll(featVisited, featI)
1458         {
1459             if (!featVisited[featI])
1460             {
1461                 startFeatI = featI;
1462
1463                 break;
1464             }
1465         }
1466
1467         if (startFeatI == -1)
1468         {
1469             // No feature lines left.
1470             break;
1471         }
1472
1473         segments.append
1474         (
1475             collectSegment
1476             (
1477                 isFeaturePoint,
1478                 featureToEdge_[startFeatI],
1479                 featVisited
1480             )
1481         );
1482     }
1483     while (true);
1484
1485
1486     //
1487     // Store in *this
1488     //
1489     featureSegments_.setSize(segments.size());
1490
1491     forAll(featureSegments_, segmentI)
1492     {
1493         featureSegments_[segmentI] = segments[segmentI];
1494     }
1495 }
1496
1497
1498 void Foam::boundaryMesh::setExtraEdges(const label edgeI)
1499 {
1500     labelList minDistance(mesh().nEdges(), -1);
1501
1502     // All edge labels encountered
1503     DynamicList<label> visitedEdges;
1504
1505     // Floodfill from edgeI starting from distance 0. Stop at distance.
1506     markEdges(8, edgeI, 0, minDistance, visitedEdges);
1507
1508     // Set edge labels to display
1509     extraEdges_.transfer(visitedEdges);
1510 }
1511
1512
1513 Foam::label Foam::boundaryMesh::whichPatch(const label faceI) const
1514 {
1515     forAll(patches_, patchI)
1516     {
1517         const boundaryPatch& bp = patches_[patchI];
1518
1519         if ((faceI >= bp.start()) && (faceI < (bp.start() + bp.size())))
1520         {
1521             return patchI;
1522         }
1523     }
1524
1525     FatalErrorIn("boundaryMesh::whichPatch(const label)")
1526         << "Cannot find face " << faceI << " in list of boundaryPatches "
1527         << patches_
1528         << abort(FatalError);
1529
1530     return -1;
1531 }
1532
1533
1534 Foam::label Foam::boundaryMesh::findPatchID(const word& patchName) const
1535 {
1536     forAll(patches_, patchI)
1537     {
1538         if (patches_[patchI].name() == patchName)
1539         {
1540             return patchI;
1541         }
1542     }
1543
1544     return -1;
1545 }
1546
1547
1548 void Foam::boundaryMesh::addPatch(const word& patchName)
1549 {
1550     patches_.setSize(patches_.size() + 1);
1551
1552     // Add empty patch at end of patch list.
1553
1554     label patchI = patches_.size()-1;
1555
1556     boundaryPatch* bpPtr = new boundaryPatch
1557     (
1558         patchName,
1559         patchI,
1560         0,
1561         mesh().size(),
1562         "empty"
1563     );
1564
1565     patches_.set(patchI, bpPtr);
1566
1567     if (debug)
1568     {
1569         Pout<< "addPatch : patches now:" << endl;
1570
1571         forAll(patches_, patchI)
1572         {
1573             const boundaryPatch& bp = patches_[patchI];
1574
1575             Pout<< "    name  : " << bp.name() << endl
1576                 << "    size  : " << bp.size() << endl
1577                 << "    start : " << bp.start() << endl
1578                 << "    type  : " << bp.physicalType() << endl
1579                 << endl;
1580         }
1581     }
1582 }
1583
1584
1585 void Foam::boundaryMesh::deletePatch(const word& patchName)
1586 {
1587     label delPatchI = findPatchID(patchName);
1588
1589     if (delPatchI == -1)
1590     {
1591         FatalErrorIn("boundaryMesh::deletePatch(const word&")
1592             << "Can't find patch named " << patchName
1593             << abort(FatalError);
1594     }
1595
1596     if (patches_[delPatchI].size())
1597     {
1598         FatalErrorIn("boundaryMesh::deletePatch(const word&")
1599             << "Trying to delete non-empty patch " << patchName
1600             << endl << "Current size:" << patches_[delPatchI].size()
1601             << abort(FatalError);
1602     }
1603
1604     PtrList<boundaryPatch> newPatches(patches_.size() - 1);
1605
1606     for (label patchI = 0; patchI < delPatchI; patchI++)
1607     {
1608         newPatches.set(patchI, patches_[patchI].clone());
1609     }
1610
1611     // Move patches down, starting from delPatchI.
1612
1613     for (label patchI = delPatchI + 1; patchI < patches_.size(); patchI++)
1614     {
1615         newPatches.set(patchI - 1, patches_[patchI].clone());
1616     }
1617
1618     patches_.clear();
1619
1620     patches_ = newPatches;
1621
1622     if (debug)
1623     {
1624         Pout<< "deletePatch : patches now:" << endl;
1625
1626         forAll(patches_, patchI)
1627         {
1628             const boundaryPatch& bp = patches_[patchI];
1629
1630             Pout<< "    name  : " << bp.name() << endl
1631                 << "    size  : " << bp.size() << endl
1632                 << "    start : " << bp.start() << endl
1633                 << "    type  : " << bp.physicalType() << endl
1634                 << endl;
1635         }
1636     }
1637 }
1638
1639
1640 void Foam::boundaryMesh::changePatchType
1641 (
1642     const word& patchName,
1643     const word& patchType
1644 )
1645 {
1646     label changeI = findPatchID(patchName);
1647
1648     if (changeI == -1)
1649     {
1650         FatalErrorIn("boundaryMesh::changePatchType(const word&, const word&)")
1651             << "Can't find patch named " << patchName
1652             << abort(FatalError);
1653     }
1654
1655
1656     // Cause we can't reassign to individual PtrList elems ;-(
1657     // work on copy
1658
1659
1660     PtrList<boundaryPatch> newPatches(patches_.size());
1661
1662     forAll(patches_, patchI)
1663     {
1664         if (patchI == changeI)
1665         {
1666             // Create copy but for type
1667             const boundaryPatch& oldBp = patches_[patchI];
1668
1669             boundaryPatch* bpPtr = new boundaryPatch
1670             (
1671                 oldBp.name(),
1672                 oldBp.index(),
1673                 oldBp.size(),
1674                 oldBp.start(),
1675                 patchType
1676             );
1677
1678             newPatches.set(patchI, bpPtr);
1679         }
1680         else
1681         {
1682             // Create copy
1683             newPatches.set(patchI, patches_[patchI].clone());
1684         }
1685     }
1686
1687     patches_ = newPatches;
1688 }
1689
1690
1691 void Foam::boundaryMesh::changeFaces
1692 (
1693     const labelList& patchIDs,
1694     labelList& oldToNew
1695 )
1696 {
1697     if (patchIDs.size() != mesh().size())
1698     {
1699         FatalErrorIn("boundaryMesh::changeFaces(const labelList& patchIDs)")
1700             << "List of patchIDs not equal to number of faces." << endl
1701             << "PatchIDs size:" << patchIDs.size()
1702             << " nFaces::" << mesh().size()
1703             << abort(FatalError);
1704     }
1705
1706     // Count number of faces for each patch
1707
1708     labelList nFaces(patches_.size(), 0);
1709
1710     forAll(patchIDs, faceI)
1711     {
1712         label patchID = patchIDs[faceI];
1713
1714         if (patchID < 0 || patchID >= patches_.size())
1715         {
1716             FatalErrorIn("boundaryMesh::changeFaces(const labelList&)")
1717                 << "PatchID " << patchID << " out of range"
1718                 << abort(FatalError);
1719         }
1720         nFaces[patchID]++;
1721     }
1722
1723
1724     // Determine position in faces_ for each patch
1725
1726     labelList startFace(patches_.size());
1727
1728     startFace[0] = 0;
1729
1730     for (label patchI = 1; patchI < patches_.size(); patchI++)
1731     {
1732         startFace[patchI] = startFace[patchI-1] + nFaces[patchI-1];
1733     }
1734
1735     // Update patch info
1736     PtrList<boundaryPatch> newPatches(patches_.size());
1737
1738     forAll(patches_, patchI)
1739     {
1740         const boundaryPatch& bp = patches_[patchI];
1741
1742         newPatches.set
1743         (
1744             patchI,
1745             new boundaryPatch
1746             (
1747                 bp.name(),
1748                 patchI,
1749                 nFaces[patchI],
1750                 startFace[patchI],
1751                 bp.physicalType()
1752             )
1753         );
1754     }
1755     patches_ = newPatches;
1756
1757     if (debug)
1758     {
1759         Pout<< "changeFaces : patches now:" << endl;
1760
1761         forAll(patches_, patchI)
1762         {
1763             const boundaryPatch& bp = patches_[patchI];
1764
1765             Pout<< "    name  : " << bp.name() << endl
1766                 << "    size  : " << bp.size() << endl
1767                 << "    start : " << bp.start() << endl
1768                 << "    type  : " << bp.physicalType() << endl
1769                 << endl;
1770         }
1771     }
1772
1773
1774     // Construct face mapping array
1775     oldToNew.setSize(patchIDs.size());
1776
1777     forAll(patchIDs, faceI)
1778     {
1779         int patchID = patchIDs[faceI];
1780
1781         oldToNew[faceI] = startFace[patchID]++;
1782     }
1783
1784     // Copy faces into correct position and maintain label of original face
1785     faceList newFaces(mesh().size());
1786
1787     labelList newMeshFace(mesh().size());
1788
1789     forAll(oldToNew, faceI)
1790     {
1791         newFaces[oldToNew[faceI]] = mesh()[faceI];
1792         newMeshFace[oldToNew[faceI]] = meshFace_[faceI];
1793     }
1794
1795     // Reconstruct 'mesh' from new faces and (copy of) existing points.
1796     bMesh* newMeshPtr_ = new bMesh(newFaces, mesh().points());
1797
1798     // Reset meshFace_ to new ordering.
1799     meshFace_.transfer(newMeshFace);
1800
1801
1802     // Remove old PrimitivePatch on meshPtr_.
1803     clearOut();
1804
1805     // And insert new 'mesh'.
1806     meshPtr_ = newMeshPtr_;
1807 }
1808
1809
1810 Foam::label Foam::boundaryMesh::getNTris(const label faceI) const
1811 {
1812     const face& f = mesh()[faceI];
1813
1814     return f.nTriangles(mesh().points());
1815 }
1816
1817
1818 Foam::label Foam::boundaryMesh::getNTris
1819 (
1820     const label startFaceI,
1821     const label nFaces,
1822     labelList& nTris
1823 ) const
1824 {
1825     label totalNTris = 0;
1826
1827     nTris.setSize(nFaces);
1828
1829     for (label i = 0; i < nFaces; i++)
1830     {
1831         label faceNTris = getNTris(startFaceI + i);
1832
1833         nTris[i] = faceNTris;
1834
1835         totalNTris += faceNTris;
1836     }
1837     return totalNTris;
1838 }
1839
1840
1841 // Simple triangulation of face subset. Stores vertices in tris[] as three
1842 // consecutive vertices per triangle.
1843 void Foam::boundaryMesh::triangulate
1844 (
1845     const label startFaceI,
1846     const label nFaces,
1847     const label totalNTris,
1848     labelList& triVerts
1849 ) const
1850 {
1851     // Triangulate faces.
1852     triVerts.setSize(3*totalNTris);
1853
1854     label vertI = 0;
1855
1856     for (label i = 0; i < nFaces; i++)
1857     {
1858         label faceI = startFaceI + i;
1859
1860         const face& f = mesh()[faceI];
1861
1862         // Have face triangulate itself (results in faceList)
1863         faceList triFaces(f.nTriangles(mesh().points()));
1864
1865         label nTri = 0;
1866
1867         f.triangles(mesh().points(), nTri, triFaces);
1868
1869         // Copy into triVerts
1870
1871         forAll(triFaces, triFaceI)
1872         {
1873             const face& triF = triFaces[triFaceI];
1874
1875             triVerts[vertI++] = triF[0];
1876             triVerts[vertI++] = triF[1];
1877             triVerts[vertI++] = triF[2];
1878         }
1879     }
1880 }
1881
1882
1883 // Number of local points in subset.
1884 Foam::label Foam::boundaryMesh::getNPoints
1885 (
1886     const label startFaceI,
1887     const label nFaces
1888 ) const
1889 {
1890     SubList<face> patchFaces(mesh(), nFaces, startFaceI);
1891
1892     primitivePatch patch(patchFaces, mesh().points());
1893
1894     return patch.nPoints();
1895 }
1896
1897
1898 // Triangulation of face subset in local coords.
1899 void Foam::boundaryMesh::triangulateLocal
1900 (
1901     const label startFaceI,
1902     const label nFaces,
1903     const label totalNTris,
1904     labelList& triVerts,
1905     labelList& localToGlobal
1906 ) const
1907 {
1908     SubList<face> patchFaces(mesh(), nFaces, startFaceI);
1909
1910     primitivePatch patch(patchFaces, mesh().points());
1911
1912     // Map from local to mesh().points() addressing
1913     localToGlobal = patch.meshPoints();
1914
1915     // Triangulate patch faces.
1916     triVerts.setSize(3*totalNTris);
1917
1918     label vertI = 0;
1919
1920     for (label i = 0; i < nFaces; i++)
1921     {
1922         // Local face
1923         const face& f = patch.localFaces()[i];
1924
1925         // Have face triangulate itself (results in faceList)
1926         faceList triFaces(f.nTriangles(patch.localPoints()));
1927
1928         label nTri = 0;
1929
1930         f.triangles(patch.localPoints(), nTri, triFaces);
1931
1932         // Copy into triVerts
1933
1934         forAll(triFaces, triFaceI)
1935         {
1936             const face& triF = triFaces[triFaceI];
1937
1938             triVerts[vertI++] = triF[0];
1939             triVerts[vertI++] = triF[1];
1940             triVerts[vertI++] = triF[2];
1941         }
1942     }
1943 }
1944
1945
1946 void Foam::boundaryMesh::markFaces
1947 (
1948     const labelList& protectedEdges,
1949     const label seedFaceI,
1950     boolList& visited
1951 ) const
1952 {
1953     boolList protectedEdge(mesh().nEdges(), false);
1954
1955     forAll(protectedEdges, i)
1956     {
1957         protectedEdge[protectedEdges[i]] = true;
1958     }
1959
1960
1961     // Initialize zone for all faces to -1
1962     labelList currentZone(mesh().size(), -1);
1963
1964     // Mark with 0 all faces reachable from seedFaceI
1965     markZone(protectedEdge, seedFaceI, 0, currentZone);
1966
1967     // Set in visited all reached ones.
1968     visited.setSize(mesh().size());
1969
1970     forAll(currentZone, faceI)
1971     {
1972         if (currentZone[faceI] == 0)
1973         {
1974             visited[faceI] = true;
1975         }
1976         else
1977         {
1978             visited[faceI] = false;
1979         }
1980     }
1981 }
1982
1983
1984 // ************************************************************************* //