src/surfMesh/MeshedSurface/MeshedSurface.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "MeshedSurface.H"
  27 #include "UnsortedMeshedSurface.H"
  28 #include "MeshedSurfaceProxy.H"
  29 #include "mergePoints.H"
  30 #include "foamTime.H"
  31 #include "ListOps.H"
  32 #include "polyBoundaryMesh.H"
  33 #include "polyMesh.H"
  34 #include "surfMesh.H"
  35 #include "Xfer.H"
  36 #include "surfZone.H"
  37 #include "primitivePatch.H"
  38 #include "addToRunTimeSelectionTable.H"
  39
  40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  41
  42 template<class Face>
  43 inline bool Foam::MeshedSurface<Face>::isTri()
  44 {
  45     return false;
  46 }
  47
  48
  49 template<class Face>
  50 Foam::wordHashSet Foam::MeshedSurface<Face>::readTypes()
  51 {
  52     return wordHashSet(*fileExtensionConstructorTablePtr_);
  53 }
  54
  55
  56 template<class Face>
  57 Foam::wordHashSet Foam::MeshedSurface<Face>::writeTypes()
  58 {
  59     return wordHashSet(*writefileExtensionMemberFunctionTablePtr_);
  60 }
  61
  62
  63 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
  64
  65 template<class Face>
  66 bool Foam::MeshedSurface<Face>::canReadType
  67 (
  68     const word& ext,
  69     const bool verbose
  70 )
  71 {
  72     return checkSupport
  73     (
  74         readTypes() | FriendType::readTypes(),
  75         ext,
  76         verbose,
  77         "reading"
  78    );
  79 }
  80
  81
  82 template<class Face>
  83 bool Foam::MeshedSurface<Face>::canWriteType
  84 (
  85     const word& ext,
  86     const bool verbose
  87 )
  88 {
  89     return checkSupport
  90     (
  91         writeTypes() | ProxyType::writeTypes(),
  92         ext,
  93         verbose,
  94         "writing"
  95     );
  96 }
  97
  98
  99 template<class Face>
 100 bool Foam::MeshedSurface<Face>::canRead
 101 (
 102     const fileName& name,
 103     const bool verbose
 104 )
 105 {
 106     word ext = name.ext();
 107     if (ext == "gz")
 108     {
 109         ext = name.lessExt().ext();
 110     }
 111     return canReadType(ext, verbose);
 112 }
 113
 114
 115 template<class Face>
 116 void Foam::MeshedSurface<Face>::write
 117 (
 118     const fileName& name,
 119     const MeshedSurface<Face>& surf
 120 )
 121 {
 122     if (debug)
 123     {
 124         Info<< "MeshedSurface::write"
 125             "(const fileName&, const MeshedSurface&) : "
 126             "writing to " << name
 127             << endl;
 128     }
 129
 130     const word ext = name.ext();
 131
 132     typename writefileExtensionMemberFunctionTable::iterator mfIter =
 133         writefileExtensionMemberFunctionTablePtr_->find(ext);
 134
 135     if (mfIter == writefileExtensionMemberFunctionTablePtr_->end())
 136     {
 137         // no direct writer, delegate to proxy if possible
 138         wordHashSet supported = ProxyType::writeTypes();
 139
 140         if (supported.found(ext))
 141         {
 142             MeshedSurfaceProxy<Face>(surf).write(name);
 143         }
 144         else
 145         {
 146             FatalErrorIn
 147             (
 148                 "MeshedSurface::write"
 149                 "(const fileName&, const MeshedSurface&)"
 150             )   << "Unknown file extension " << ext << nl << nl
 151                 << "Valid types are :" << endl
 152                 << (supported | writeTypes())
 153                 << exit(FatalError);
 154         }
 155     }
 156     else
 157     {
 158         mfIter()(name, surf);
 159     }
 160 }
 161
 162
 163 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 164
 165 template<class Face>
 166 Foam::MeshedSurface<Face>::MeshedSurface()
 167 :
 168     ParentType(List<Face>(), pointField())
 169 {}
 170
 171
 172 template<class Face>
 173 Foam::MeshedSurface<Face>::MeshedSurface
 174 (
 175     const Xfer<pointField >& pointLst,
 176     const Xfer<List<Face> >& faceLst,
 177     const Xfer<surfZoneList>& zoneLst
 178 )
 179 :
 180     ParentType(List<Face>(), pointField()),
 181     zones_()
 182 {
 183     reset(pointLst, faceLst, zoneLst);
 184 }
 185
 186
 187 template<class Face>
 188 Foam::MeshedSurface<Face>::MeshedSurface
 189 (
 190     const Xfer<pointField>& pointLst,
 191     const Xfer<List<Face> >& faceLst,
 192     const UList<label>& zoneSizes,
 193     const UList<word>& zoneNames
 194 )
 195 :
 196     ParentType(List<Face>(), pointField())
 197 {
 198     reset(pointLst, faceLst, Xfer<surfZoneList>());
 199
 200     if (zoneSizes.size())
 201     {
 202         if (zoneNames.size())
 203         {
 204             addZones(zoneSizes, zoneNames);
 205         }
 206         else
 207         {
 208             addZones(zoneSizes);
 209         }
 210     }
 211 }
 212
 213
 214 template<class Face>
 215 Foam::MeshedSurface<Face>::MeshedSurface
 216 (
 217     const MeshedSurface<Face>& surf
 218 )
 219 :
 220     ParentType(surf.faces(), surf.points()),
 221     zones_(surf.surfZones())
 222 {}
 223
 224
 225 template<class Face>
 226 Foam::MeshedSurface<Face>::MeshedSurface
 227 (
 228     const UnsortedMeshedSurface<Face>& surf
 229 )
 230 :
 231     ParentType(List<Face>(), surf.points())
 232 {
 233     labelList faceMap;
 234     this->storedZones().transfer(surf.sortedZones(faceMap));
 235
 236     const List<Face>& origFaces = surf.faces();
 237     List<Face> newFaces(origFaces.size());
 238
 239     // this is somewhat like ListOps reorder and/or IndirectList
 240     forAll(newFaces, faceI)
 241     {
 242         newFaces[faceI] = origFaces[faceMap[faceI]];
 243     }
 244
 245     this->storedFaces().transfer(newFaces);
 246 }
 247
 248
 249 template<class Face>
 250 Foam::MeshedSurface<Face>::MeshedSurface(const surfMesh& mesh)
 251 :
 252     ParentType(List<Face>(), pointField())
 253 {
 254     // same face type as surfMesh
 255     MeshedSurface<face> surf
 256     (
 257         xferCopy(mesh.points()),
 258         xferCopy(mesh.faces()),
 259         xferCopy(mesh.surfZones())
 260     );
 261
 262     this->transcribe(surf);
 263 }
 264
 265
 266 template<class Face>
 267 Foam::MeshedSurface<Face>::MeshedSurface
 268 (
 269     const polyBoundaryMesh& bMesh,
 270     const bool useGlobalPoints
 271 )
 272 :
 273     ParentType(List<Face>(), pointField())
 274 {
 275     const polyMesh& mesh = bMesh.mesh();
 276     const polyPatchList& bPatches = bMesh;
 277
 278     // Get a single patch for all boundaries
 279     primitivePatch allBoundary
 280     (
 281         SubList<face>
 282         (
 283             mesh.faces(),
 284             mesh.nFaces() - mesh.nInternalFaces(),
 285             mesh.nInternalFaces()
 286         ),
 287         mesh.points()
 288     );
 289
 290     // use global/local points:
 291     const pointField& bPoints =
 292     (
 293         useGlobalPoints ? mesh.points() : allBoundary.localPoints()
 294     );
 295
 296     // global/local face addressing:
 297     const List<Face>& bFaces =
 298     (
 299         useGlobalPoints ? allBoundary : allBoundary.localFaces()
 300     );
 301
 302
 303     // create zone list
 304     surfZoneList newZones(bPatches.size());
 305
 306     label startFaceI = 0;
 307     label nZone = 0;
 308     forAll(bPatches, patchI)
 309     {
 310         const polyPatch& p = bPatches[patchI];
 311
 312         if (p.size())
 313         {
 314             newZones[nZone] = surfZone
 315             (
 316                 p.name(),
 317                 p.size(),
 318                 startFaceI,
 319                 nZone
 320             );
 321
 322             nZone++;
 323             startFaceI += p.size();
 324         }
 325     }
 326
 327     newZones.setSize(nZone);
 328
 329     // same face type as the polyBoundaryMesh
 330     MeshedSurface<face> surf
 331     (
 332         xferCopy(bPoints),
 333         xferCopy(bFaces),
 334         xferMove(newZones)
 335     );
 336
 337     this->transcribe(surf);
 338 }
 339
 340
 341 template<class Face>
 342 Foam::MeshedSurface<Face>::MeshedSurface
 343 (
 344     const fileName& name,
 345     const word& ext
 346 )
 347 :
 348     ParentType(List<Face>(), pointField())
 349 {
 350     read(name, ext);
 351 }
 352
 353
 354 template<class Face>
 355 Foam::MeshedSurface<Face>::MeshedSurface(const fileName& name)
 356 :
 357     ParentType(List<Face>(), pointField())
 358 {
 359     read(name);
 360 }
 361
 362
 363 template<class Face>
 364 Foam::MeshedSurface<Face>::MeshedSurface
 365 (
 366     const Time& t,
 367     const word& surfName
 368 )
 369 :
 370     ParentType(List<Face>(), pointField())
 371 {
 372     surfMesh mesh
 373     (
 374         IOobject
 375         (
 376             "dummyName",
 377             t.timeName(),
 378             t,
 379             IOobject::MUST_READ,
 380             IOobject::NO_WRITE,
 381             false
 382         ),
 383         surfName
 384     );
 385
 386     // same face type as surfMesh
 387     MeshedSurface<face> surf
 388     (
 389         xferMove(mesh.storedPoints()),
 390         xferMove(mesh.storedFaces()),
 391         xferMove(mesh.storedZones())
 392     );
 393
 394     this->transcribe(surf);
 395 }
 396
 397
 398 template<class Face>
 399 Foam::MeshedSurface<Face>::MeshedSurface
 400 (
 401     const Xfer< UnsortedMeshedSurface<Face> >& surf
 402 )
 403 :
 404     ParentType(List<Face>(), pointField())
 405 {
 406     transfer(surf());
 407 }
 408
 409
 410 template<class Face>
 411 Foam::MeshedSurface<Face>::MeshedSurface
 412 (
 413     const Xfer< MeshedSurface<Face> >& surf
 414 )
 415 :
 416     ParentType(List<Face>(), pointField())
 417 {
 418     transfer(surf());
 419 }
 420
 421
 422
 423 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 424
 425 template<class Face>
 426 Foam::MeshedSurface<Face>::~MeshedSurface()
 427 {}
 428
 429
 430 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 431
 432
 433 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 434
 435 template<class Face>
 436 void Foam::MeshedSurface<Face>::remapFaces
 437 (
 438     const UList<label>& faceMap
 439 )
 440 {
 441     // recalculate the zone start/size
 442     if (!faceMap.empty())
 443     {
 444         surfZoneList& zones = storedZones();
 445
 446         if (zones.size() == 1)
 447         {
 448             // optimized for single zone case
 449             zones[0].size() = faceMap.size();
 450         }
 451         else if (zones.size())
 452         {
 453             label newFaceI = 0;
 454             label origEndI = 0;
 455             forAll(zones, zoneI)
 456             {
 457                 surfZone& zone = zones[zoneI];
 458
 459                 // adjust zone start
 460                 zone.start() = newFaceI;
 461                 origEndI += zone.size();
 462
 463                 for (label faceI = newFaceI; faceI < faceMap.size(); ++faceI)
 464                 {
 465                     if (faceMap[faceI] < origEndI)
 466                     {
 467                         ++newFaceI;
 468                     }
 469                     else
 470                     {
 471                         break;
 472                     }
 473                 }
 474
 475                 // adjust zone size
 476                 zone.size() = newFaceI - zone.start();
 477             }
 478         }
 479     }
 480 }
 481
 482
 483 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 484
 485 template<class Face>
 486 void Foam::MeshedSurface<Face>::clear()
 487 {
 488     ParentType::clearOut();
 489
 490     storedPoints().clear();
 491     storedFaces().clear();
 492     storedZones().clear();
 493 }
 494
 495
 496 template<class Face>
 497 void Foam::MeshedSurface<Face>::movePoints(const pointField& newPoints)
 498 {
 499     // Remove all geometry dependent data
 500     ParentType::clearTopology();
 501
 502     // Adapt for new point position
 503     ParentType::movePoints(newPoints);
 504
 505     // Copy new points
 506     storedPoints() = newPoints;
 507 }
 508
 509
 510 template<class Face>
 511 void Foam::MeshedSurface<Face>::scalePoints(const scalar& scaleFactor)
 512 {
 513     // avoid bad scaling
 514     if (scaleFactor > 0 && scaleFactor != 1.0)
 515     {
 516         // Remove all geometry dependent data
 517         ParentType::clearTopology();
 518
 519         // Adapt for new point position
 520         ParentType::movePoints(pointField());
 521
 522         storedPoints() *= scaleFactor;
 523     }
 524 }
 525
 526
 527 template<class Face>
 528 void Foam::MeshedSurface<Face>::reset
 529 (
 530     const Xfer<pointField>& pointLst,
 531     const Xfer<List<Face> >& faceLst,
 532     const Xfer<surfZoneList>& zoneLst
 533 )
 534 {
 535     ParentType::clearOut();
 536
 537     // Take over new primitive data.
 538     // Optimized to avoid overwriting data at all
 539     if (!pointLst->empty())
 540     {
 541         storedPoints().transfer(pointLst());
 542     }
 543
 544     if (!faceLst->empty())
 545     {
 546         storedFaces().transfer(faceLst());
 547     }
 548
 549     if (!zoneLst->empty())
 550     {
 551         storedZones().transfer(zoneLst());
 552     }
 553 }
 554
 555
 556 template<class Face>
 557 void Foam::MeshedSurface<Face>::reset
 558 (
 559     const Xfer<List<point> >& pointLst,
 560     const Xfer<List<Face> >& faceLst,
 561     const Xfer<surfZoneList>& zoneLst
 562 )
 563 {
 564     ParentType::clearOut();
 565
 566     // Take over new primitive data.
 567     // Optimized to avoid overwriting data at all
 568     if (!pointLst->empty())
 569     {
 570         storedPoints().transfer(pointLst());
 571     }
 572
 573     if (!faceLst->empty())
 574     {
 575         storedFaces().transfer(faceLst());
 576     }
 577
 578     if (!zoneLst->empty())
 579     {
 580         storedZones().transfer(zoneLst());
 581     }
 582 }
 583
 584
 585 // Remove badly degenerate faces, double faces.
 586 template<class Face>
 587 void Foam::MeshedSurface<Face>::cleanup(const bool verbose)
 588 {
 589     // merge points (already done for STL, TRI)
 590     stitchFaces(SMALL, verbose);
 591
 592     checkFaces(verbose);
 593     this->checkTopology(verbose);
 594 }
 595
 596
 597 template<class Face>
 598 bool Foam::MeshedSurface<Face>::stitchFaces
 599 (
 600     const scalar tol,
 601     const bool verbose
 602 )
 603 {
 604     pointField& pointLst = this->storedPoints();
 605
 606     // Merge points
 607     labelList  pointMap(pointLst.size());
 608     pointField newPoints(pointLst.size());
 609
 610     bool hasMerged = mergePoints(pointLst, tol, verbose, pointMap, newPoints);
 611
 612     if (!hasMerged)
 613     {
 614         return false;
 615     }
 616
 617     if (verbose)
 618     {
 619         Info<< "MeshedSurface::stitchFaces : Renumbering all faces"
 620             << endl;
 621     }
 622
 623     // Set the coordinates to the merged ones
 624     pointLst.transfer(newPoints);
 625
 626     List<Face>& faceLst = this->storedFaces();
 627
 628     List<label> faceMap(faceLst.size());
 629
 630     // Reset the point labels to the unique points array
 631     label newFaceI = 0;
 632     forAll(faceLst, faceI)
 633     {
 634         Face& f = faceLst[faceI];
 635         forAll(f, fp)
 636         {
 637             f[fp] = pointMap[f[fp]];
 638         }
 639
 640         // for extra safety: collapse face as well
 641         if (f.collapse() >= 3)
 642         {
 643             if (newFaceI != faceI)
 644             {
 645                 faceLst[newFaceI] = f;
 646             }
 647             faceMap[newFaceI] = faceI;
 648             newFaceI++;
 649         }
 650         else if (verbose)
 651         {
 652             Pout<< "MeshedSurface::stitchFaces : "
 653                 << "Removing collapsed face " << faceI << endl
 654                 << "    vertices   :" << f << endl;
 655         }
 656     }
 657     pointMap.clear();
 658
 659     if (newFaceI != faceLst.size())
 660     {
 661         if (verbose)
 662         {
 663             Pout<< "MeshedSurface::stitchFaces : "
 664                 << "Removed " << faceLst.size() - newFaceI
 665                 << " faces" << endl;
 666         }
 667         faceLst.setSize(newFaceI);
 668         remapFaces(faceMap);
 669     }
 670     faceMap.clear();
 671
 672     // Merging points might have changed geometric factors
 673     ParentType::clearOut();
 674     return true;
 675 }
 676
 677
 678 // Remove badly degenerate faces and double faces.
 679 template<class Face>
 680 bool Foam::MeshedSurface<Face>::checkFaces
 681 (
 682     const bool verbose
 683 )
 684 {
 685     bool changed = false;
 686     List<Face>& faceLst = this->storedFaces();
 687
 688     List<label> faceMap(faceLst.size());
 689
 690     label newFaceI = 0;
 691     // Detect badly labelled faces and mark degenerate faces
 692     const label maxPointI = this->points().size() - 1;
 693     forAll(faceLst, faceI)
 694     {
 695         Face& f = faceLst[faceI];
 696
 697         // avoid degenerate faces
 698         if (f.collapse() >= 3)
 699         {
 700             forAll(f, fp)
 701             {
 702                 if (f[fp] < 0 || f[fp] > maxPointI)
 703                 {
 704                     FatalErrorIn("MeshedSurface::checkFaces(bool)")
 705                         << "face " << f
 706                         << " uses point indices outside point range 0.."
 707                     << maxPointI
 708                         << exit(FatalError);
 709                 }
 710             }
 711
 712             faceMap[faceI] = faceI;
 713             newFaceI++;
 714         }
 715         else
 716         {
 717             // mark as bad face
 718             faceMap[faceI] = -1;
 719
 720             changed = true;
 721             if (verbose)
 722             {
 723                 WarningIn
 724                 (
 725                     "MeshedSurface::checkFaces(bool verbose)"
 726                 )   << "face[" << faceI << "] = " << f
 727                     << " does not have three unique vertices" << endl;
 728             }
 729         }
 730     }
 731
 732     // Detect doubled faces
 733     // do not touch the faces
 734     const labelListList& fFaces = this->faceFaces();
 735     newFaceI = 0;
 736     forAll(faceLst, faceI)
 737     {
 738         // skip already collapsed faces:
 739         if (faceMap[faceI] < 0)
 740         {
 741             continue;
 742         }
 743
 744         const Face& f = faceLst[faceI];
 745
 746         // duplicate face check
 747         bool okay = true;
 748         const labelList& neighbours = fFaces[faceI];
 749
 750         // Check if faceNeighbours use same points as this face.
 751         // Note: discards normal information - sides of baffle are merged.
 752         forAll(neighbours, neighI)
 753         {
 754             const label neiFaceI = neighbours[neighI];
 755
 756             if (neiFaceI <= faceI || faceMap[neiFaceI] < 0)
 757             {
 758                 // lower numbered faces already checked
 759                 // skip neighbours that are themselves collapsed
 760                 continue;
 761             }
 762
 763             const Face& nei = faceLst[neiFaceI];
 764
 765             if (f == nei)
 766             {
 767                 okay = false;
 768
 769                 if (verbose)
 770                 {
 771                     WarningIn
 772                     (
 773                         "MeshedSurface::checkFaces(bool verbose)"
 774                     )   << "faces share the same vertices:" << nl
 775                         << "    face[" << faceI << "] : " << f << nl
 776                         << "    face[" << neiFaceI << "] : " << nei << endl;
 777                     // printFace(Warning, "    ", f, points());
 778                     // printFace(Warning, "    ", nei, points());
 779                 }
 780
 781                 break;
 782             }
 783         }
 784
 785         if (okay)
 786         {
 787             faceMap[faceI] = faceI;
 788             newFaceI++;
 789         }
 790         else
 791         {
 792             faceMap[faceI] = -1;
 793         }
 794     }
 795
 796     // Phase 1: pack
 797     // Done to keep numbering constant in phase 1
 798
 799     if (changed || newFaceI < faceLst.size())
 800     {
 801         changed = true;
 802
 803         if (verbose)
 804         {
 805             WarningIn
 806             (
 807                 "MeshedSurface::checkFaces(bool verbose)"
 808             )   << "Removed " << faceLst.size() - newFaceI
 809                 << " illegal faces." << endl;
 810         }
 811
 812         // compress the face list
 813         newFaceI = 0;
 814         forAll(faceLst, faceI)
 815         {
 816             if (faceMap[faceI] >= 0)
 817             {
 818                 if (newFaceI != faceI)
 819                 {
 820                     faceLst[newFaceI] = faceLst[faceI];
 821                 }
 822                 faceMap[newFaceI] = faceI;
 823                 newFaceI++;
 824             }
 825         }
 826
 827         faceLst.setSize(newFaceI);
 828         remapFaces(faceMap);
 829     }
 830     faceMap.clear();
 831
 832     // Topology can change because of renumbering
 833     ParentType::clearOut();
 834     return changed;
 835 }
 836
 837
 838 template<class Face>
 839 Foam::label Foam::MeshedSurface<Face>::triangulate()
 840 {
 841     return triangulate
 842     (
 843         const_cast<List<label>&>(List<label>::null())
 844     );
 845 }
 846
 847
 848 template<class Face>
 849 Foam::label Foam::MeshedSurface<Face>::triangulate
 850 (
 851     List<label>& faceMapOut
 852 )
 853 {
 854     label nTri = 0;
 855     label maxTri = 0;  // the maximum number of triangles for any single face
 856     List<Face>& faceLst = this->storedFaces();
 857
 858     // determine how many triangles will be needed
 859     forAll(faceLst, faceI)
 860     {
 861         const label n = faceLst[faceI].nTriangles();
 862         if (maxTri < n)
 863         {
 864             maxTri = n;
 865         }
 866         nTri += n;
 867     }
 868
 869     // nothing to do
 870     if (nTri <= faceLst.size())
 871     {
 872         if (!faceMapOut.empty())
 873         {
 874             faceMapOut.clear();
 875         }
 876         return 0;
 877     }
 878
 879     List<Face>  newFaces(nTri);
 880     List<label> faceMap;
 881
 882     // reuse storage from optional faceMap
 883     if (!faceMapOut.empty())
 884     {
 885         faceMap.transfer(faceMapOut);
 886     }
 887     faceMap.setSize(nTri);
 888
 889     // remember the number of *additional* faces
 890     nTri -= faceLst.size();
 891
 892     if (this->points().empty())
 893     {
 894         // triangulate without points
 895         // simple face triangulation around f[0]
 896         label newFaceI = 0;
 897         forAll(faceLst, faceI)
 898         {
 899             const Face& f = faceLst[faceI];
 900
 901             for (label fp = 1; fp < f.size() - 1; ++fp)
 902             {
 903                 label fp1 = f.fcIndex(fp);
 904
 905                 newFaces[newFaceI] = triFace(f[0], f[fp], f[fp1]);
 906                 faceMap[newFaceI] = faceI;
 907                 newFaceI++;
 908             }
 909         }
 910     }
 911     else
 912     {
 913         // triangulate with points
 914         List<face> tmpTri(maxTri);
 915
 916         label newFaceI = 0;
 917         forAll(faceLst, faceI)
 918         {
 919             // 'face' not '<Face>'
 920             const face& f = faceLst[faceI];
 921
 922             label nTmp = 0;
 923             f.triangles(this->points(), nTmp, tmpTri);
 924             for (label triI = 0; triI < nTmp; triI++)
 925             {
 926                 newFaces[newFaceI] = Face
 927                 (
 928                     static_cast<UList<label>&>(tmpTri[triI])
 929                 );
 930                 faceMap[newFaceI] = faceI;
 931                 newFaceI++;
 932             }
 933         }
 934     }
 935
 936     faceLst.transfer(newFaces);
 937     remapFaces(faceMap);
 938
 939     // optionally return the faceMap
 940     faceMapOut.transfer(faceMap);
 941     faceMap.clear();
 942
 943     // Topology can change because of renumbering
 944     ParentType::clearOut();
 945     return nTri;
 946 }
 947
 948
 949
 950
 951 template<class Face>
 952 Foam::MeshedSurface<Face> Foam::MeshedSurface<Face>::subsetMesh
 953 (
 954     const labelHashSet& include,
 955     labelList& pointMap,
 956     labelList& faceMap
 957 ) const
 958 {
 959     const pointField& locPoints = this->localPoints();
 960     const List<Face>& locFaces  = this->localFaces();
 961
 962
 963     // Fill pointMap, faceMap
 964     PatchTools::subsetMap(*this, include, pointMap, faceMap);
 965
 966     // Create compact coordinate list and forward mapping array
 967     pointField newPoints(pointMap.size());
 968     labelList oldToNew(locPoints.size());
 969     forAll(pointMap, pointI)
 970     {
 971         newPoints[pointI] = locPoints[pointMap[pointI]];
 972         oldToNew[pointMap[pointI]] = pointI;
 973     }
 974
 975     // create/copy a new zones list, each zone with zero size
 976     surfZoneList newZones(this->surfZones());
 977     forAll(newZones, zoneI)
 978     {
 979         newZones[zoneI].size() = 0;
 980     }
 981
 982     // Renumber face node labels
 983     List<Face> newFaces(faceMap.size());
 984     forAll(faceMap, faceI)
 985     {
 986         const label origFaceI = faceMap[faceI];
 987         newFaces[faceI] = Face(locFaces[origFaceI]);
 988
 989         // Renumber labels for face
 990         Face& f = newFaces[faceI];
 991         forAll(f, fp)
 992         {
 993             f[fp] = oldToNew[f[fp]];
 994         }
 995     }
 996     oldToNew.clear();
 997
 998     // recalculate the zones start/size
 999     label newFaceI = 0;
1000     label origEndI = 0;
1001
1002     // adjust zone sizes
1003     forAll(newZones, zoneI)
1004     {
1005         surfZone& zone = newZones[zoneI];
1006
1007         // adjust zone start
1008         zone.start() = newFaceI;
1009         origEndI += zone.size();
1010
1011         for (label faceI = newFaceI; faceI < faceMap.size(); ++faceI)
1012         {
1013             if (faceMap[faceI] < origEndI)
1014             {
1015                 ++newFaceI;
1016             }
1017             else
1018             {
1019                 break;
1020             }
1021         }
1022
1023         // adjust zone size
1024         zone.size() = newFaceI - zone.start();
1025     }
1026
1027
1028     // construct a sub-surface
1029     return MeshedSurface
1030     (
1031         xferMove(newPoints),
1032         xferMove(newFaces),
1033         xferMove(newZones)
1034     );
1035 }
1036
1037
1038 template<class Face>
1039 Foam::MeshedSurface<Face>
1040 Foam::MeshedSurface<Face>::subsetMesh
1041 (
1042     const labelHashSet& include
1043 ) const
1044 {
1045     labelList pointMap, faceMap;
1046     return subsetMesh(include, pointMap, faceMap);
1047 }
1048
1049
1050
1051 template<class Face>
1052 void Foam::MeshedSurface<Face>::transfer
1053 (
1054     MeshedSurface<Face>& surf
1055 )
1056 {
1057     reset
1058     (
1059         xferMove(surf.storedPoints()),
1060         xferMove(surf.storedFaces()),
1061         xferMove(surf.storedZones())
1062     );
1063 }
1064
1065
1066 template<class Face>
1067 void Foam::MeshedSurface<Face>::transfer
1068 (
1069     UnsortedMeshedSurface<Face>& surf
1070 )
1071 {
1072     clear();
1073
1074     labelList faceMap;
1075     surfZoneList zoneLst = surf.sortedZones(faceMap);
1076
1077     if (zoneLst.size() <= 1)
1078     {
1079         reset
1080         (
1081             xferMove(surf.storedPoints()),
1082             xferMove(surf.storedFaces()),
1083             Xfer<surfZoneList>()
1084         );
1085     }
1086     else
1087     {
1088         List<Face>& oldFaces = surf.storedFaces();
1089         List<Face> newFaces(faceMap.size());
1090
1091         forAll(faceMap, faceI)
1092         {
1093             newFaces[faceI].transfer(oldFaces[faceMap[faceI]]);
1094         }
1095
1096         reset
1097         (
1098             xferMove(surf.storedPoints()),
1099             xferMove(newFaces),
1100             xferMove(zoneLst)
1101         );
1102     }
1103
1104     faceMap.clear();
1105     surf.clear();
1106 }
1107
1108
1109 template<class Face>
1110 Foam::Xfer< Foam::MeshedSurface<Face> >
1111 Foam::MeshedSurface<Face>::xfer()
1112 {
1113     return xferMove(*this);
1114 }
1115
1116
1117 // Read from file, determine format from extension
1118 template<class Face>
1119 bool Foam::MeshedSurface<Face>::read(const fileName& name)
1120 {
1121     word ext = name.ext();
1122     if (ext == "gz")
1123     {
1124         fileName unzipName = name.lessExt();
1125         return read(unzipName, unzipName.ext());
1126     }
1127     else
1128     {
1129         return read(name, ext);
1130     }
1131 }
1132
1133
1134 // Read from file in given format
1135 template<class Face>
1136 bool Foam::MeshedSurface<Face>::read
1137 (
1138     const fileName& name,
1139     const word& ext
1140 )
1141 {
1142     clear();
1143
1144     // read via selector mechanism
1145     transfer(New(name, ext)());
1146     return true;
1147 }
1148
1149
1150 template<class Face>
1151 void Foam::MeshedSurface<Face>::write
1152 (
1153     const Time& t,
1154     const word& surfName
1155 ) const
1156 {
1157     MeshedSurfaceProxy<Face>(*this).write(t, surfName);
1158 }
1159
1160 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
1161
1162 template<class Face>
1163 void Foam::MeshedSurface<Face>::operator=(const MeshedSurface& surf)
1164 {
1165     clear();
1166
1167     this->storedPoints() = surf.points();
1168     this->storedFaces()  = surf.faces();
1169     this->storedZones()  = surf.surfZones();
1170 }
1171
1172
1173 template<class Face>
1174 Foam::MeshedSurface<Face>::operator
1175 Foam::MeshedSurfaceProxy<Face>() const
1176 {
1177     return MeshedSurfaceProxy<Face>
1178     (
1179         this->points(),
1180         this->faces(),
1181         this->surfZones()
1182     );
1183 }
1184
1185 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
1186
1187 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
1188
1189 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1190
1191 #include "MeshedSurfaceZones.C"
1192 #include "MeshedSurfaceIO.C"
1193 #include "MeshedSurfaceNew.C"
1194
1195 // ************************************************************************* //