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