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