src/dynamicMesh/boundaryMesh/boundaryMesh.C

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