src/surfMesh/MeshedSurface/MeshedSurface.C

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