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