src/surfMesh/MeshedSurface/MeshedSurface.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "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 fileFormats::surfaceFormatsCore::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 fileFormats::surfaceFormatsCore::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 labelUList& 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_IF_MODIFIED,
 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 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 429
 430 template<class Face>
 431 void Foam::MeshedSurface<Face>::remapFaces
 432 (
 433     const labelUList& faceMap
 434 )
 435 {
 436     // recalculate the zone start/size
 437     if (&faceMap && faceMap.size())
 438     {
 439         surfZoneList& zones = storedZones();
 440
 441         if (zones.size() == 1)
 442         {
 443             // optimized for single zone case
 444             zones[0].size() = faceMap.size();
 445         }
 446         else if (zones.size())
 447         {
 448             label newFaceI = 0;
 449             label origEndI = 0;
 450             forAll(zones, zoneI)
 451             {
 452                 surfZone& zone = zones[zoneI];
 453
 454                 // adjust zone start
 455                 zone.start() = newFaceI;
 456                 origEndI += zone.size();
 457
 458                 for (label faceI = newFaceI; faceI < faceMap.size(); ++faceI)
 459                 {
 460                     if (faceMap[faceI] < origEndI)
 461                     {
 462                         ++newFaceI;
 463                     }
 464                     else
 465                     {
 466                         break;
 467                     }
 468                 }
 469
 470                 // adjust zone size
 471                 zone.size() = newFaceI - zone.start();
 472             }
 473         }
 474     }
 475 }
 476
 477
 478 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 479
 480 template<class Face>
 481 void Foam::MeshedSurface<Face>::clear()
 482 {
 483     ParentType::clearOut();
 484
 485     storedPoints().clear();
 486     storedFaces().clear();
 487     storedZones().clear();
 488 }
 489
 490
 491 template<class Face>
 492 void Foam::MeshedSurface<Face>::movePoints(const pointField& newPoints)
 493 {
 494     // Remove all geometry dependent data
 495     ParentType::clearTopology();
 496
 497     // Adapt for new point position
 498     ParentType::movePoints(newPoints);
 499
 500     // Copy new points
 501     storedPoints() = newPoints;
 502 }
 503
 504
 505 template<class Face>
 506 void Foam::MeshedSurface<Face>::scalePoints(const scalar scaleFactor)
 507 {
 508     // avoid bad scaling
 509     if (scaleFactor > 0 && scaleFactor != 1.0)
 510     {
 511         // Remove all geometry dependent data
 512         ParentType::clearTopology();
 513
 514         // Adapt for new point position
 515         ParentType::movePoints(pointField());
 516
 517         storedPoints() *= scaleFactor;
 518     }
 519 }
 520
 521
 522 template<class Face>
 523 void Foam::MeshedSurface<Face>::reset
 524 (
 525     const Xfer<pointField>& pointLst,
 526     const Xfer<List<Face> >& faceLst,
 527     const Xfer<surfZoneList>& zoneLst
 528 )
 529 {
 530     ParentType::clearOut();
 531
 532     // Take over new primitive data.
 533     // Optimized to avoid overwriting data at all
 534     if (&pointLst)
 535     {
 536         storedPoints().transfer(pointLst());
 537     }
 538
 539     if (&faceLst)
 540     {
 541         storedFaces().transfer(faceLst());
 542     }
 543
 544     if (&zoneLst)
 545     {
 546         storedZones().transfer(zoneLst());
 547     }
 548 }
 549
 550
 551 template<class Face>
 552 void Foam::MeshedSurface<Face>::reset
 553 (
 554     const Xfer<List<point> >& pointLst,
 555     const Xfer<List<Face> >& faceLst,
 556     const Xfer<surfZoneList>& zoneLst
 557 )
 558 {
 559     ParentType::clearOut();
 560
 561     // Take over new primitive data.
 562     // Optimized to avoid overwriting data at all
 563     if (&pointLst)
 564     {
 565         storedPoints().transfer(pointLst());
 566     }
 567
 568     if (&faceLst)
 569     {
 570         storedFaces().transfer(faceLst());
 571     }
 572
 573     if (&zoneLst)
 574     {
 575         storedZones().transfer(zoneLst());
 576     }
 577 }
 578
 579
 580 // Remove badly degenerate faces, double faces.
 581 template<class Face>
 582 void Foam::MeshedSurface<Face>::cleanup(const bool verbose)
 583 {
 584     // merge points (already done for STL, TRI)
 585     stitchFaces(SMALL, verbose);
 586
 587     checkFaces(verbose);
 588     this->checkTopology(verbose);
 589 }
 590
 591
 592 template<class Face>
 593 bool Foam::MeshedSurface<Face>::stitchFaces
 594 (
 595     const scalar tol,
 596     const bool verbose
 597 )
 598 {
 599     pointField& pointLst = this->storedPoints();
 600
 601     // Merge points
 602     labelList  pointMap(pointLst.size());
 603     pointField newPoints(pointLst.size());
 604
 605     bool hasMerged = mergePoints(pointLst, tol, verbose, pointMap, newPoints);
 606
 607     if (!hasMerged)
 608     {
 609         return false;
 610     }
 611
 612     if (verbose)
 613     {
 614         Info<< "MeshedSurface::stitchFaces : Renumbering all faces"
 615             << endl;
 616     }
 617
 618     // Set the coordinates to the merged ones
 619     pointLst.transfer(newPoints);
 620
 621     List<Face>& faceLst = this->storedFaces();
 622
 623     List<label> faceMap(faceLst.size());
 624
 625     // Reset the point labels to the unique points array
 626     label newFaceI = 0;
 627     forAll(faceLst, faceI)
 628     {
 629         Face& f = faceLst[faceI];
 630         forAll(f, fp)
 631         {
 632             f[fp] = pointMap[f[fp]];
 633         }
 634
 635         // for extra safety: collapse face as well
 636         if (f.collapse() >= 3)
 637         {
 638             if (newFaceI != faceI)
 639             {
 640                 faceLst[newFaceI] = f;
 641             }
 642             faceMap[newFaceI] = faceI;
 643             newFaceI++;
 644         }
 645         else if (verbose)
 646         {
 647             Pout<< "MeshedSurface::stitchFaces : "
 648                 << "Removing collapsed face " << faceI << endl
 649                 << "    vertices   :" << f << endl;
 650         }
 651     }
 652     pointMap.clear();
 653
 654     if (newFaceI != faceLst.size())
 655     {
 656         if (verbose)
 657         {
 658             Pout<< "MeshedSurface::stitchFaces : "
 659                 << "Removed " << faceLst.size() - newFaceI
 660                 << " faces" << endl;
 661         }
 662         faceLst.setSize(newFaceI);
 663         remapFaces(faceMap);
 664     }
 665     faceMap.clear();
 666
 667     // Merging points might have changed geometric factors
 668     ParentType::clearOut();
 669     return true;
 670 }
 671
 672
 673 // Remove badly degenerate faces and double faces.
 674 template<class Face>
 675 bool Foam::MeshedSurface<Face>::checkFaces
 676 (
 677     const bool verbose
 678 )
 679 {
 680     bool changed = false;
 681     List<Face>& faceLst = this->storedFaces();
 682
 683     List<label> faceMap(faceLst.size());
 684
 685     label newFaceI = 0;
 686     // Detect badly labelled faces and mark degenerate faces
 687     const label maxPointI = this->points().size() - 1;
 688     forAll(faceLst, faceI)
 689     {
 690         Face& f = faceLst[faceI];
 691
 692         // avoid degenerate faces
 693         if (f.collapse() >= 3)
 694         {
 695             forAll(f, fp)
 696             {
 697                 if (f[fp] < 0 || f[fp] > maxPointI)
 698                 {
 699                     FatalErrorIn("MeshedSurface::checkFaces(bool)")
 700                         << "face " << f
 701                         << " uses point indices outside point range 0.."
 702                     << maxPointI
 703                         << exit(FatalError);
 704                 }
 705             }
 706
 707             faceMap[faceI] = faceI;
 708             newFaceI++;
 709         }
 710         else
 711         {
 712             // mark as bad face
 713             faceMap[faceI] = -1;
 714
 715             changed = true;
 716             if (verbose)
 717             {
 718                 WarningIn
 719                 (
 720                     "MeshedSurface::checkFaces(bool verbose)"
 721                 )   << "face[" << faceI << "] = " << f
 722                     << " does not have three unique vertices" << endl;
 723             }
 724         }
 725     }
 726
 727     // Detect doubled faces
 728     // do not touch the faces
 729     const labelListList& fFaces = this->faceFaces();
 730     newFaceI = 0;
 731     forAll(faceLst, faceI)
 732     {
 733         // skip already collapsed faces:
 734         if (faceMap[faceI] < 0)
 735         {
 736             continue;
 737         }
 738
 739         const Face& f = faceLst[faceI];
 740
 741         // duplicate face check
 742         bool okay = true;
 743         const labelList& neighbours = fFaces[faceI];
 744
 745         // Check if faceNeighbours use same points as this face.
 746         // Note: discards normal information - sides of baffle are merged.
 747         forAll(neighbours, neighI)
 748         {
 749             const label neiFaceI = neighbours[neighI];
 750
 751             if (neiFaceI <= faceI || faceMap[neiFaceI] < 0)
 752             {
 753                 // lower numbered faces already checked
 754                 // skip neighbours that are themselves collapsed
 755                 continue;
 756             }
 757
 758             const Face& nei = faceLst[neiFaceI];
 759
 760             if (f == nei)
 761             {
 762                 okay = false;
 763
 764                 if (verbose)
 765                 {
 766                     WarningIn
 767                     (
 768                         "MeshedSurface::checkFaces(bool verbose)"
 769                     )   << "faces share the same vertices:" << nl
 770                         << "    face[" << faceI << "] : " << f << nl
 771                         << "    face[" << neiFaceI << "] : " << nei << endl;
 772                     // printFace(Warning, "    ", f, points());
 773                     // printFace(Warning, "    ", nei, points());
 774                 }
 775
 776                 break;
 777             }
 778         }
 779
 780         if (okay)
 781         {
 782             faceMap[faceI] = faceI;
 783             newFaceI++;
 784         }
 785         else
 786         {
 787             faceMap[faceI] = -1;
 788         }
 789     }
 790
 791     // Phase 1: pack
 792     // Done to keep numbering constant in phase 1
 793
 794     if (changed || newFaceI < faceLst.size())
 795     {
 796         changed = true;
 797
 798         if (verbose)
 799         {
 800             WarningIn
 801             (
 802                 "MeshedSurface::checkFaces(bool verbose)"
 803             )   << "Removed " << faceLst.size() - newFaceI
 804                 << " illegal faces." << endl;
 805         }
 806
 807         // compress the face list
 808         newFaceI = 0;
 809         forAll(faceLst, faceI)
 810         {
 811             if (faceMap[faceI] >= 0)
 812             {
 813                 if (newFaceI != faceI)
 814                 {
 815                     faceLst[newFaceI] = faceLst[faceI];
 816                 }
 817                 faceMap[newFaceI] = faceI;
 818                 newFaceI++;
 819             }
 820         }
 821
 822         faceLst.setSize(newFaceI);
 823         remapFaces(faceMap);
 824     }
 825     faceMap.clear();
 826
 827     // Topology can change because of renumbering
 828     ParentType::clearOut();
 829     return changed;
 830 }
 831
 832
 833 template<class Face>
 834 Foam::label Foam::MeshedSurface<Face>::triangulate()
 835 {
 836     return triangulate
 837     (
 838         const_cast<List<label>&>(List<label>::null())
 839     );
 840 }
 841
 842
 843 template<class Face>
 844 Foam::label Foam::MeshedSurface<Face>::triangulate
 845 (
 846     List<label>& faceMapOut
 847 )
 848 {
 849     label nTri = 0;
 850     label maxTri = 0;  // the maximum number of triangles for any single face
 851     List<Face>& faceLst = this->storedFaces();
 852
 853     // determine how many triangles will be needed
 854     forAll(faceLst, faceI)
 855     {
 856         const label n = faceLst[faceI].nTriangles();
 857         if (maxTri < n)
 858         {
 859             maxTri = n;
 860         }
 861         nTri += n;
 862     }
 863
 864     // nothing to do
 865     if (nTri <= faceLst.size())
 866     {
 867         if (&faceMapOut)
 868         {
 869             faceMapOut.clear();
 870         }
 871         return 0;
 872     }
 873
 874     List<Face>  newFaces(nTri);
 875     List<label> faceMap;
 876
 877     // reuse storage from optional faceMap
 878     if (&faceMapOut)
 879     {
 880         faceMap.transfer(faceMapOut);
 881     }
 882     faceMap.setSize(nTri);
 883
 884     // remember the number of *additional* faces
 885     nTri -= faceLst.size();
 886
 887     if (this->points().empty())
 888     {
 889         // triangulate without points
 890         // simple face triangulation around f[0]
 891         label newFaceI = 0;
 892         forAll(faceLst, faceI)
 893         {
 894             const Face& f = faceLst[faceI];
 895
 896             for (label fp = 1; fp < f.size() - 1; ++fp)
 897             {
 898                 label fp1 = f.fcIndex(fp);
 899
 900                 newFaces[newFaceI] = triFace(f[0], f[fp], f[fp1]);
 901                 faceMap[newFaceI] = faceI;
 902                 newFaceI++;
 903             }
 904         }
 905     }
 906     else
 907     {
 908         // triangulate with points
 909         List<face> tmpTri(maxTri);
 910
 911         label newFaceI = 0;
 912         forAll(faceLst, faceI)
 913         {
 914             // 'face' not '<Face>'
 915             const face& f = faceLst[faceI];
 916
 917             label nTmp = 0;
 918             f.triangles(this->points(), nTmp, tmpTri);
 919             for (label triI = 0; triI < nTmp; triI++)
 920             {
 921                 newFaces[newFaceI] = Face
 922                 (
 923                     static_cast<labelUList&>(tmpTri[triI])
 924                 );
 925                 faceMap[newFaceI] = faceI;
 926                 newFaceI++;
 927             }
 928         }
 929     }
 930
 931     faceLst.transfer(newFaces);
 932     remapFaces(faceMap);
 933
 934     // optionally return the faceMap
 935     if (&faceMapOut)
 936     {
 937         faceMapOut.transfer(faceMap);
 938     }
 939     faceMap.clear();
 940
 941     // Topology can change because of renumbering
 942     ParentType::clearOut();
 943     return nTri;
 944 }
 945
 946
 947
 948
 949 template<class Face>
 950 Foam::MeshedSurface<Face> Foam::MeshedSurface<Face>::subsetMesh
 951 (
 952     const labelHashSet& include,
 953     labelList& pointMap,
 954     labelList& faceMap
 955 ) const
 956 {
 957     const pointField& locPoints = this->localPoints();
 958     const List<Face>& locFaces  = this->localFaces();
 959
 960
 961     // Fill pointMap, faceMap
 962     PatchTools::subsetMap(*this, include, pointMap, faceMap);
 963
 964     // Create compact coordinate list and forward mapping array
 965     pointField newPoints(pointMap.size());
 966     labelList oldToNew(locPoints.size());
 967     forAll(pointMap, pointI)
 968     {
 969         newPoints[pointI] = locPoints[pointMap[pointI]];
 970         oldToNew[pointMap[pointI]] = pointI;
 971     }
 972
 973     // create/copy a new zones list, each zone with zero size
 974     surfZoneList newZones(this->surfZones());
 975     forAll(newZones, zoneI)
 976     {
 977         newZones[zoneI].size() = 0;
 978     }
 979
 980     // Renumber face node labels
 981     List<Face> newFaces(faceMap.size());
 982     forAll(faceMap, faceI)
 983     {
 984         const label origFaceI = faceMap[faceI];
 985         newFaces[faceI] = Face(locFaces[origFaceI]);
 986
 987         // Renumber labels for face
 988         Face& f = newFaces[faceI];
 989         forAll(f, fp)
 990         {
 991             f[fp] = oldToNew[f[fp]];
 992         }
 993     }
 994     oldToNew.clear();
 995
 996     // recalculate the zones start/size
 997     label newFaceI = 0;
 998     label origEndI = 0;
 999
1000     // adjust zone sizes
1001     forAll(newZones, zoneI)
1002     {
1003         surfZone& zone = newZones[zoneI];
1004
1005         // adjust zone start
1006         zone.start() = newFaceI;
1007         origEndI += zone.size();
1008
1009         for (label faceI = newFaceI; faceI < faceMap.size(); ++faceI)
1010         {
1011             if (faceMap[faceI] < origEndI)
1012             {
1013                 ++newFaceI;
1014             }
1015             else
1016             {
1017                 break;
1018             }
1019         }
1020
1021         // adjust zone size
1022         zone.size() = newFaceI - zone.start();
1023     }
1024
1025
1026     // construct a sub-surface
1027     return MeshedSurface
1028     (
1029         xferMove(newPoints),
1030         xferMove(newFaces),
1031         xferMove(newZones)
1032     );
1033 }
1034
1035
1036 template<class Face>
1037 Foam::MeshedSurface<Face> Foam::MeshedSurface<Face>::subsetMesh
1038 (
1039     const labelHashSet& include
1040 ) const
1041 {
1042     labelList pointMap, faceMap;
1043     return subsetMesh(include, pointMap, faceMap);
1044 }
1045
1046
1047
1048 template<class Face>
1049 void Foam::MeshedSurface<Face>::transfer
1050 (
1051     MeshedSurface<Face>& surf
1052 )
1053 {
1054     reset
1055     (
1056         xferMove(surf.storedPoints()),
1057         xferMove(surf.storedFaces()),
1058         xferMove(surf.storedZones())
1059     );
1060 }
1061
1062
1063 template<class Face>
1064 void Foam::MeshedSurface<Face>::transfer
1065 (
1066     UnsortedMeshedSurface<Face>& surf
1067 )
1068 {
1069     clear();
1070
1071     labelList faceMap;
1072     surfZoneList zoneLst = surf.sortedZones(faceMap);
1073
1074     if (zoneLst.size() <= 1)
1075     {
1076         reset
1077         (
1078             xferMove(surf.storedPoints()),
1079             xferMove(surf.storedFaces()),
1080             Xfer<surfZoneList>()
1081         );
1082     }
1083     else
1084     {
1085         List<Face>& oldFaces = surf.storedFaces();
1086         List<Face> newFaces(faceMap.size());
1087
1088         forAll(faceMap, faceI)
1089         {
1090             newFaces[faceI].transfer(oldFaces[faceMap[faceI]]);
1091         }
1092
1093         reset
1094         (
1095             xferMove(surf.storedPoints()),
1096             xferMove(newFaces),
1097             xferMove(zoneLst)
1098         );
1099     }
1100
1101     faceMap.clear();
1102     surf.clear();
1103 }
1104
1105
1106 template<class Face>
1107 Foam::Xfer<Foam::MeshedSurface<Face> > Foam::MeshedSurface<Face>::xfer()
1108 {
1109     return xferMove(*this);
1110 }
1111
1112
1113 // Read from file, determine format from extension
1114 template<class Face>
1115 bool Foam::MeshedSurface<Face>::read(const fileName& name)
1116 {
1117     word ext = name.ext();
1118     if (ext == "gz")
1119     {
1120         fileName unzipName = name.lessExt();
1121         return read(unzipName, unzipName.ext());
1122     }
1123     else
1124     {
1125         return read(name, ext);
1126     }
1127 }
1128
1129
1130 // Read from file in given format
1131 template<class Face>
1132 bool Foam::MeshedSurface<Face>::read
1133 (
1134     const fileName& name,
1135     const word& ext
1136 )
1137 {
1138     clear();
1139
1140     // read via selector mechanism
1141     transfer(New(name, ext)());
1142     return true;
1143 }
1144
1145
1146 template<class Face>
1147 void Foam::MeshedSurface<Face>::write
1148 (
1149     const Time& t,
1150     const word& surfName
1151 ) const
1152 {
1153     MeshedSurfaceProxy<Face>(*this).write(t, surfName);
1154 }
1155
1156
1157 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
1158
1159 template<class Face>
1160 void Foam::MeshedSurface<Face>::operator=(const MeshedSurface& surf)
1161 {
1162     clear();
1163
1164     this->storedPoints() = surf.points();
1165     this->storedFaces()  = surf.faces();
1166     this->storedZones()  = surf.surfZones();
1167 }
1168
1169
1170 template<class Face>
1171 Foam::MeshedSurface<Face>::operator Foam::MeshedSurfaceProxy<Face>() const
1172 {
1173     return MeshedSurfaceProxy<Face>
1174     (
1175         this->points(),
1176         this->faces(),
1177         this->surfZones()
1178     );
1179 }
1180
1181
1182 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1183
1184 #include "MeshedSurfaceZones.C"
1185 #include "MeshedSurfaceIO.C"
1186 #include "MeshedSurfaceNew.C"
1187
1188 // ************************************************************************* //