src/triSurface/triSurface/triSurface.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 "triSurface.H"
  28 #include "demandDrivenData.H"
  29 #include "IFstream.H"
  30 #include "OFstream.H"
  31 #include "Time.H"
  32 #include "boundBox.H"
  33 #include "SortableList.H"
  34 #include "PackedBoolList.H"
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39 namespace Foam
  40 {
  41     defineTypeNameAndDebug(Foam::triSurface, 0);
  42 }
  43
  44
  45 Foam::fileName Foam::triSurface::triSurfInstance(const Time& d)
  46 {
  47     fileName foamName(d.caseName() + ".ftr");
  48
  49     // Search back through the time directories list to find the time
  50     // closest to and lower than current time
  51
  52     instantList ts = d.times();
  53     label i;
  54
  55     for (i=ts.size()-1; i>=0; i--)
  56     {
  57         if (ts[i].value() <= d.timeOutputValue())
  58         {
  59             break;
  60         }
  61     }
  62
  63     // Noting that the current directory has already been searched
  64     // for mesh data, start searching from the previously stored time directory
  65
  66     if (i>=0)
  67     {
  68         for (label j=i; j>=0; j--)
  69         {
  70             if (isFile(d.path()/ts[j].name()/typeName/foamName))
  71             {
  72                 if (debug)
  73                 {
  74                     Pout<< " triSurface::triSurfInstance(const Time& d)"
  75                         << "reading " << foamName
  76                         << " from " << ts[j].name()/typeName
  77                         << endl;
  78                 }
  79
  80                 return ts[j].name();
  81             }
  82         }
  83     }
  84
  85     if (debug)
  86     {
  87         Pout<< " triSurface::triSurfInstance(const Time& d)"
  88             << "reading " << foamName
  89             << " from constant/" << endl;
  90     }
  91     return "constant";
  92 }
  93
  94
  95 Foam::List<Foam::labelledTri> Foam::triSurface::convertToTri
  96 (
  97     const faceList& faces,
  98     const label defaultRegion
  99 )
 100 {
 101     List<labelledTri> triFaces(faces.size());
 102
 103     forAll(triFaces, faceI)
 104     {
 105         const face& f = faces[faceI];
 106
 107         if (f.size() != 3)
 108         {
 109             FatalErrorIn
 110             (
 111                 "triSurface::convertToTri"
 112                 "(const faceList&, const label)"
 113             )   << "Face at position " << faceI
 114                 << " does not have three vertices:" << f
 115                 << abort(FatalError);
 116         }
 117
 118         labelledTri& tri = triFaces[faceI];
 119
 120         tri[0] = f[0];
 121         tri[1] = f[1];
 122         tri[2] = f[2];
 123         tri.region() = defaultRegion;
 124     }
 125
 126     return triFaces;
 127 }
 128
 129
 130 Foam::List<Foam::labelledTri> Foam::triSurface::convertToTri
 131 (
 132     const triFaceList& faces,
 133     const label defaultRegion
 134 )
 135 {
 136     List<labelledTri> triFaces(faces.size());
 137
 138     forAll(triFaces, faceI)
 139     {
 140         const triFace& f = faces[faceI];
 141
 142         labelledTri& tri = triFaces[faceI];
 143
 144         tri[0] = f[0];
 145         tri[1] = f[1];
 146         tri[2] = f[2];
 147         tri.region() = defaultRegion;
 148     }
 149
 150     return triFaces;
 151 }
 152
 153
 154 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 155
 156 void Foam::triSurface::printTriangle
 157 (
 158     Ostream& os,
 159     const string& pre,
 160     const labelledTri& f,
 161     const pointField& points
 162 )
 163 {
 164     os
 165         << pre.c_str() << "vertex numbers:"
 166         << f[0] << ' ' << f[1] << ' ' << f[2] << endl
 167         << pre.c_str() << "vertex coords :"
 168         << points[f[0]] << ' ' << points[f[1]] << ' ' << points[f[2]]
 169         << pre.c_str() << "region        :" << f.region() << endl
 170         << endl;
 171 }
 172
 173
 174 Foam::string Foam::triSurface::getLineNoComment(IFstream& is)
 175 {
 176     string line;
 177     do
 178     {
 179         is.getLine(line);
 180     }
 181     while ((line.empty() || line[0] == '#') && is.good());
 182
 183     return line;
 184 }
 185
 186
 187 // Remove non-triangles, double triangles.
 188 void Foam::triSurface::checkTriangles(const bool verbose)
 189 {
 190     // Simple check on indices ok.
 191     const label maxPointI = points().size() - 1;
 192
 193     forAll(*this, faceI)
 194     {
 195         const labelledTri& f = (*this)[faceI];
 196
 197         if
 198         (
 199             (f[0] < 0) || (f[0] > maxPointI)
 200          || (f[1] < 0) || (f[1] > maxPointI)
 201          || (f[2] < 0) || (f[2] > maxPointI)
 202         )
 203         {
 204             FatalErrorIn("triSurface::checkTriangles(bool)")
 205                 << "triangle " << f
 206                 << " uses point indices outside point range 0.."
 207                 << maxPointI
 208                 << exit(FatalError);
 209         }
 210     }
 211
 212     // Two phase process
 213     //   1. mark invalid faces
 214     //   2. pack
 215     // Done to keep numbering constant in phase 1
 216
 217     // List of valid triangles
 218     boolList valid(size(), true);
 219     bool hasInvalid = false;
 220
 221     const labelListList& fFaces = faceFaces();
 222
 223     forAll(*this, faceI)
 224     {
 225         const labelledTri& f = (*this)[faceI];
 226
 227         if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
 228         {
 229             // 'degenerate' triangle check
 230             valid[faceI] = false;
 231             hasInvalid = true;
 232
 233             if (verbose)
 234             {
 235                 WarningIn
 236                 (
 237                     "triSurface::checkTriangles(bool verbose)"
 238                 )   << "triangle " << faceI
 239                     << " does not have three unique vertices:\n";
 240                 printTriangle(Warning, "    ", f, points());
 241             }
 242         }
 243         else
 244         {
 245             // duplicate triangle check
 246             const labelList& neighbours = fFaces[faceI];
 247
 248             // Check if faceNeighbours use same points as this face.
 249             // Note: discards normal information - sides of baffle are merged.
 250             forAll(neighbours, neighbourI)
 251             {
 252                 if (neighbours[neighbourI] <= faceI)
 253                 {
 254                     // lower numbered faces already checked
 255                     continue;
 256                 }
 257
 258                 const labelledTri& n = (*this)[neighbours[neighbourI]];
 259
 260                 if
 261                 (
 262                     ((f[0] == n[0]) || (f[0] == n[1]) || (f[0] == n[2]))
 263                  && ((f[1] == n[0]) || (f[1] == n[1]) || (f[1] == n[2]))
 264                  && ((f[2] == n[0]) || (f[2] == n[1]) || (f[2] == n[2]))
 265                 )
 266                 {
 267                     valid[faceI] = false;
 268                     hasInvalid = true;
 269
 270                     if (verbose)
 271                     {
 272                         WarningIn
 273                         (
 274                             "triSurface::checkTriangles(bool verbose)"
 275                         )   << "triangles share the same vertices:\n"
 276                             << "    face 1 :" << faceI << endl;
 277                         printTriangle(Warning, "    ", f, points());
 278
 279                         Warning
 280                             << endl
 281                             << "    face 2 :"
 282                             << neighbours[neighbourI] << endl;
 283                         printTriangle(Warning, "    ", n, points());
 284                     }
 285
 286                     break;
 287                 }
 288             }
 289         }
 290     }
 291
 292     if (hasInvalid)
 293     {
 294         // Pack
 295         label newFaceI = 0;
 296         forAll(*this, faceI)
 297         {
 298             if (valid[faceI])
 299             {
 300                 const labelledTri& f = (*this)[faceI];
 301                 (*this)[newFaceI++] = f;
 302             }
 303         }
 304
 305         if (verbose)
 306         {
 307             WarningIn
 308             (
 309                 "triSurface::checkTriangles(bool verbose)"
 310             )   << "Removing " << size() - newFaceI
 311                 << " illegal faces." << endl;
 312         }
 313         (*this).setSize(newFaceI);
 314
 315         // Topology can change because of renumbering
 316         clearOut();
 317     }
 318 }
 319
 320
 321 // Check/fix edges with more than two triangles
 322 void Foam::triSurface::checkEdges(const bool verbose)
 323 {
 324     const labelListList& eFaces = edgeFaces();
 325
 326     forAll(eFaces, edgeI)
 327     {
 328         const labelList& myFaces = eFaces[edgeI];
 329
 330         if (myFaces.empty())
 331         {
 332             FatalErrorIn("triSurface::checkEdges(bool verbose)")
 333                 << "Edge " << edgeI << " with vertices " << edges()[edgeI]
 334                 << " has no edgeFaces"
 335                 << exit(FatalError);
 336         }
 337         else if (myFaces.size() > 2)
 338         {
 339             WarningIn
 340             (
 341                 "triSurface::checkEdges(bool verbose)"
 342             )   << "Edge " << edgeI << " with vertices " << edges()[edgeI]
 343                 << " has more than 2 faces connected to it : " << myFaces
 344                 << endl;
 345         }
 346     }
 347 }
 348
 349
 350 // Read triangles, points from Istream
 351 bool Foam::triSurface::read(Istream& is)
 352 {
 353     is  >> patches_ >> storedPoints() >> storedFaces();
 354
 355     return true;
 356 }
 357
 358
 359 // Read from file in given format
 360 bool Foam::triSurface::read
 361 (
 362     const fileName& name,
 363     const word& ext,
 364     const bool check
 365 )
 366 {
 367     if (check && !exists(name))
 368     {
 369         FatalErrorIn
 370         (
 371             "triSurface::read(const fileName&, const word&, const bool)"
 372         )   << "Cannnot read " << name << exit(FatalError);
 373     }
 374
 375     if (ext == "gz")
 376     {
 377         fileName unzipName = name.lessExt();
 378
 379         // Do not check for existence. Let IFstream do the unzipping.
 380         return read(unzipName, unzipName.ext(), false);
 381     }
 382     else if (ext == "ftr")
 383     {
 384         return read(IFstream(name)());
 385     }
 386     else if (ext == "stl")
 387     {
 388         return readSTL(name);
 389     }
 390     else if (ext == "stlb")
 391     {
 392         return readSTL(name);
 393     }
 394     else if (ext == "gts")
 395     {
 396         return readGTS(name);
 397     }
 398     else if (ext == "obj")
 399     {
 400         return readOBJ(name);
 401     }
 402     else if (ext == "off")
 403     {
 404         return readOFF(name);
 405     }
 406     else if (ext == "tri")
 407     {
 408         return readTRI(name);
 409     }
 410     else if (ext == "ac")
 411     {
 412         return readAC(name);
 413     }
 414     else if (ext == "nas")
 415     {
 416         return readNAS(name);
 417     }
 418     else
 419     {
 420         FatalErrorIn
 421         (
 422             "triSurface::read(const fileName&, const word&)"
 423         )   << "unknown file extension " << ext
 424             << ". Supported extensions are '.ftr', '.stl', '.stlb', '.gts'"
 425             << ", '.obj', '.ac', '.off', '.nas' and '.tri'"
 426             << exit(FatalError);
 427
 428         return false;
 429     }
 430 }
 431
 432
 433 // Write to file in given format
 434 void Foam::triSurface::write
 435 (
 436     const fileName& name,
 437     const word& ext,
 438     const bool sort
 439 ) const
 440 {
 441     if (ext == "ftr")
 442     {
 443         return write(OFstream(name)());
 444     }
 445     else if (ext == "stl")
 446     {
 447         return writeSTLASCII(OFstream(name)());
 448     }
 449     else if (ext == "stlb")
 450     {
 451         ofstream outFile(name.c_str(), std::ios::binary);
 452
 453         writeSTLBINARY(outFile);
 454     }
 455     else if (ext == "gts")
 456     {
 457         return writeGTS(sort, OFstream(name)());
 458     }
 459     else if (ext == "obj")
 460     {
 461         writeOBJ(sort, OFstream(name)());
 462     }
 463     else if (ext == "off")
 464     {
 465         writeOFF(sort, OFstream(name)());
 466     }
 467     else if (ext == "vtk")
 468     {
 469         writeVTK(sort, OFstream(name)());
 470     }
 471     else if (ext == "tri")
 472     {
 473         writeTRI(sort, OFstream(name)());
 474     }
 475     else if (ext == "dx")
 476     {
 477         writeDX(sort, OFstream(name)());
 478     }
 479     else if (ext == "ac")
 480     {
 481         writeAC(OFstream(name)());
 482     }
 483     else if (ext == "smesh")
 484     {
 485         writeSMESH(sort, OFstream(name)());
 486     }
 487     else
 488     {
 489         FatalErrorIn
 490         (
 491             "triSurface::write(const fileName&, const word&, const bool)"
 492         )   << "unknown file extension " << ext
 493             << ". Supported extensions are '.ftr', '.stl', '.stlb', "
 494             << "'.gts', '.obj', '.vtk'"
 495             << ", '.off', '.dx', '.smesh', '.ac' and '.tri'"
 496             << exit(FatalError);
 497     }
 498 }
 499
 500
 501 // Returns patch info. Sets faceMap to the indexing according to patch
 502 // numbers. Patch numbers start at 0.
 503 Foam::surfacePatchList Foam::triSurface::calcPatches(labelList& faceMap) const
 504 {
 505     // Sort according to region numbers of labelledTri
 506     SortableList<label> sortedRegion(size());
 507
 508     forAll(sortedRegion, faceI)
 509     {
 510         sortedRegion[faceI] = operator[](faceI).region();
 511     }
 512     sortedRegion.sort();
 513
 514     faceMap = sortedRegion.indices();
 515
 516     // Extend regions
 517     label maxRegion = patches_.size()-1;    // for non-compacted regions
 518
 519     if (faceMap.size())
 520     {
 521         maxRegion = max
 522         (
 523             maxRegion,
 524             operator[](faceMap[faceMap.size() - 1]).region()
 525         );
 526     }
 527
 528     // Get new region list
 529     surfacePatchList newPatches(maxRegion + 1);
 530
 531     // Fill patch sizes
 532     forAll(*this, faceI)
 533     {
 534         label region = operator[](faceI).region();
 535
 536         newPatches[region].size()++;
 537     }
 538
 539
 540     // Fill rest of patch info
 541
 542     label startFaceI = 0;
 543     forAll(newPatches, newPatchI)
 544     {
 545         surfacePatch& newPatch = newPatches[newPatchI];
 546
 547         label oldPatchI = newPatchI;
 548
 549         // start of patch
 550         newPatch.start() = startFaceI;
 551
 552
 553         // Take over any information from existing patches
 554         if ((oldPatchI < patches_.size()) && (patches_[oldPatchI].name() != ""))
 555         {
 556             newPatch.name() = patches_[oldPatchI].name();
 557         }
 558         else
 559         {
 560             newPatch.name() = word("patch") + name(newPatchI);
 561         }
 562
 563         if
 564         (
 565             (oldPatchI < patches_.size())
 566          && (patches_[oldPatchI].geometricType() != "")
 567         )
 568         {
 569             newPatch.geometricType() = patches_[oldPatchI].geometricType();
 570         }
 571         else
 572         {
 573             newPatch.geometricType() = "empty";
 574         }
 575
 576         startFaceI += newPatch.size();
 577     }
 578
 579     return newPatches;
 580 }
 581
 582
 583 void Foam::triSurface::setDefaultPatches()
 584 {
 585     labelList faceMap;
 586
 587     // Get names, types and sizes
 588     surfacePatchList newPatches(calcPatches(faceMap));
 589
 590     // Take over names and types (but not size)
 591     patches_.setSize(newPatches.size());
 592
 593     forAll(newPatches, patchI)
 594     {
 595         patches_[patchI].name() = newPatches[patchI].name();
 596         patches_[patchI].geometricType() = newPatches[patchI].geometricType();
 597     }
 598 }
 599
 600
 601 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 602
 603 Foam::triSurface::triSurface()
 604 :
 605     ParentType(List<Face>(), pointField()),
 606     patches_(0),
 607     sortedEdgeFacesPtr_(NULL),
 608     edgeOwnerPtr_(NULL)
 609 {}
 610
 611
 612
 613 Foam::triSurface::triSurface
 614 (
 615     const List<labelledTri>& triangles,
 616     const geometricSurfacePatchList& patches,
 617     const pointField& points
 618 )
 619 :
 620     ParentType(triangles, points),
 621     patches_(patches),
 622     sortedEdgeFacesPtr_(NULL),
 623     edgeOwnerPtr_(NULL)
 624 {}
 625
 626
 627 Foam::triSurface::triSurface
 628 (
 629     List<labelledTri>& triangles,
 630     const geometricSurfacePatchList& patches,
 631     pointField& points,
 632     const bool reUse
 633 )
 634 :
 635     ParentType(triangles, points, reUse),
 636     patches_(patches),
 637     sortedEdgeFacesPtr_(NULL),
 638     edgeOwnerPtr_(NULL)
 639 {}
 640
 641
 642 Foam::triSurface::triSurface
 643 (
 644     const List<labelledTri>& triangles,
 645     const pointField& points
 646 )
 647 :
 648     ParentType(triangles, points),
 649     patches_(),
 650     sortedEdgeFacesPtr_(NULL),
 651     edgeOwnerPtr_(NULL)
 652 {
 653     setDefaultPatches();
 654 }
 655
 656
 657 Foam::triSurface::triSurface
 658 (
 659     const triFaceList& triangles,
 660     const pointField& points
 661 )
 662 :
 663     ParentType(convertToTri(triangles, 0), points),
 664     patches_(),
 665     sortedEdgeFacesPtr_(NULL),
 666     edgeOwnerPtr_(NULL)
 667 {
 668     setDefaultPatches();
 669 }
 670
 671
 672 Foam::triSurface::triSurface(const fileName& name)
 673 :
 674     ParentType(List<Face>(), pointField()),
 675     patches_(),
 676     sortedEdgeFacesPtr_(NULL),
 677     edgeOwnerPtr_(NULL)
 678 {
 679     word ext = name.ext();
 680
 681     read(name, ext);
 682
 683     setDefaultPatches();
 684 }
 685
 686
 687 Foam::triSurface::triSurface(Istream& is)
 688 :
 689     ParentType(List<Face>(), pointField()),
 690     patches_(),
 691     sortedEdgeFacesPtr_(NULL),
 692     edgeOwnerPtr_(NULL)
 693 {
 694     read(is);
 695
 696     setDefaultPatches();
 697 }
 698
 699
 700 Foam::triSurface::triSurface(const Time& d)
 701 :
 702     ParentType(List<Face>(), pointField()),
 703     patches_(),
 704     sortedEdgeFacesPtr_(NULL),
 705     edgeOwnerPtr_(NULL)
 706 {
 707     fileName foamFile(d.caseName() + ".ftr");
 708
 709     fileName foamPath(d.path()/triSurfInstance(d)/typeName/foamFile);
 710
 711     IFstream foamStream(foamPath);
 712
 713     read(foamStream);
 714
 715     setDefaultPatches();
 716 }
 717
 718
 719 Foam::triSurface::triSurface(const triSurface& ts)
 720 :
 721     ParentType(ts, ts.points()),
 722     patches_(ts.patches()),
 723     sortedEdgeFacesPtr_(NULL),
 724     edgeOwnerPtr_(NULL)
 725 {}
 726
 727
 728 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 729
 730 Foam::triSurface::~triSurface()
 731 {
 732     clearOut();
 733 }
 734
 735
 736 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 737
 738 void Foam::triSurface::clearTopology()
 739 {
 740     ParentType::clearTopology();
 741     deleteDemandDrivenData(sortedEdgeFacesPtr_);
 742     deleteDemandDrivenData(edgeOwnerPtr_);
 743 }
 744
 745
 746 void Foam::triSurface::clearPatchMeshAddr()
 747 {
 748     ParentType::clearPatchMeshAddr();
 749 }
 750
 751
 752 void Foam::triSurface::clearOut()
 753 {
 754     ParentType::clearOut();
 755
 756     clearTopology();
 757     clearPatchMeshAddr();
 758 }
 759
 760
 761 const Foam::labelListList& Foam::triSurface::sortedEdgeFaces() const
 762 {
 763     if (!sortedEdgeFacesPtr_)
 764     {
 765         calcSortedEdgeFaces();
 766     }
 767
 768     return *sortedEdgeFacesPtr_;
 769 }
 770
 771
 772 const Foam::labelList& Foam::triSurface::edgeOwner() const
 773 {
 774     if (!edgeOwnerPtr_)
 775     {
 776         calcEdgeOwner();
 777     }
 778
 779     return *edgeOwnerPtr_;
 780 }
 781
 782
 783 //- Move points
 784 void Foam::triSurface::movePoints(const pointField& newPoints)
 785 {
 786     // Remove all geometry dependent data
 787     deleteDemandDrivenData(sortedEdgeFacesPtr_);
 788
 789     // Adapt for new point position
 790     ParentType::movePoints(newPoints);
 791
 792     // Copy new points
 793     storedPoints() = newPoints;
 794 }
 795
 796
 797 // scale points
 798 void Foam::triSurface::scalePoints(const scalar& scaleFactor)
 799 {
 800     // avoid bad scaling
 801     if (scaleFactor > 0 && scaleFactor != 1.0)
 802     {
 803         // Remove all geometry dependent data
 804         clearTopology();
 805
 806         // Adapt for new point position
 807         ParentType::movePoints(pointField());
 808
 809         storedPoints() *= scaleFactor;
 810     }
 811 }
 812
 813
 814 // Remove non-triangles, double triangles.
 815 void Foam::triSurface::cleanup(const bool verbose)
 816 {
 817     // Merge points (already done for STL, TRI)
 818     stitchTriangles(pointField(points()), SMALL, verbose);
 819
 820     // Merging points might have changed geometric factors
 821     clearOut();
 822
 823     checkTriangles(verbose);
 824
 825     checkEdges(verbose);
 826 }
 827
 828
 829 // Finds area, starting at faceI, delimited by borderEdge. Marks all visited
 830 // faces (from face-edge-face walk) with currentZone.
 831 void Foam::triSurface::markZone
 832 (
 833     const boolList& borderEdge,
 834     const label faceI,
 835     const label currentZone,
 836     labelList& faceZone
 837 ) const
 838 {
 839     // List of faces whose faceZone has been set.
 840     labelList changedFaces(1, faceI);
 841
 842     while(true)
 843     {
 844         // Pick up neighbours of changedFaces
 845         DynamicList<label> newChangedFaces(2*changedFaces.size());
 846
 847         forAll(changedFaces, i)
 848         {
 849             label faceI = changedFaces[i];
 850
 851             const labelList& fEdges = faceEdges()[faceI];
 852
 853             forAll(fEdges, i)
 854             {
 855                 label edgeI = fEdges[i];
 856
 857                 if (!borderEdge[edgeI])
 858                 {
 859                     const labelList& eFaces = edgeFaces()[edgeI];
 860
 861                     forAll(eFaces, j)
 862                     {
 863                         label nbrFaceI = eFaces[j];
 864
 865                         if (faceZone[nbrFaceI] == -1)
 866                         {
 867                             faceZone[nbrFaceI] = currentZone;
 868                             newChangedFaces.append(nbrFaceI);
 869                         }
 870                         else if (faceZone[nbrFaceI] != currentZone)
 871                         {
 872                             FatalErrorIn
 873                             (
 874                                 "triSurface::markZone(const boolList&,"
 875                                 "const label, const label, labelList&) const"
 876                             )
 877                                 << "Zones " << faceZone[nbrFaceI]
 878                                 << " at face " << nbrFaceI
 879                                 << " connects to zone " << currentZone
 880                                 << " at face " << faceI
 881                                 << abort(FatalError);
 882                         }
 883                     }
 884                 }
 885             }
 886         }
 887
 888         if (newChangedFaces.empty())
 889         {
 890             break;
 891         }
 892
 893         changedFaces.transfer(newChangedFaces);
 894     }
 895 }
 896
 897
 898 // Finds areas delimited by borderEdge (or 'real' edges).
 899 // Fills faceZone accordingly
 900 Foam::label Foam::triSurface::markZones
 901 (
 902     const boolList& borderEdge,
 903     labelList& faceZone
 904 ) const
 905 {
 906     faceZone.setSize(size());
 907     faceZone = -1;
 908
 909     if (borderEdge.size() != nEdges())
 910     {
 911         FatalErrorIn
 912         (
 913             "triSurface::markZones"
 914             "(const boolList&, labelList&)"
 915         )
 916             << "borderEdge boolList not same size as number of edges" << endl
 917             << "borderEdge:" << borderEdge.size() << endl
 918             << "nEdges    :" << nEdges()
 919             << exit(FatalError);
 920     }
 921
 922     label zoneI = 0;
 923
 924     for (label startFaceI = 0;; zoneI++)
 925     {
 926         // Find first non-coloured face
 927         for (; startFaceI < size(); startFaceI++)
 928         {
 929             if (faceZone[startFaceI] == -1)
 930             {
 931                 break;
 932             }
 933         }
 934
 935         if (startFaceI >= size())
 936         {
 937             break;
 938         }
 939
 940         faceZone[startFaceI] = zoneI;
 941
 942         markZone(borderEdge, startFaceI, zoneI, faceZone);
 943     }
 944
 945     return zoneI;
 946 }
 947
 948
 949 void Foam::triSurface::subsetMeshMap
 950 (
 951     const boolList& include,
 952     labelList& pointMap,
 953     labelList& faceMap
 954 ) const
 955 {
 956     const List<labelledTri>& locFaces = localFaces();
 957
 958     label faceI = 0;
 959     label pointI = 0;
 960
 961     faceMap.setSize(locFaces.size());
 962
 963     pointMap.setSize(nPoints());
 964
 965     boolList pointHad(nPoints(), false);
 966
 967     forAll(include, oldFacei)
 968     {
 969         if (include[oldFacei])
 970         {
 971             // Store new faces compact
 972             faceMap[faceI++] = oldFacei;
 973
 974             // Renumber labels for triangle
 975             const labelledTri& tri = locFaces[oldFacei];
 976
 977             label a = tri[0];
 978             if (!pointHad[a])
 979             {
 980                 pointHad[a] = true;
 981                 pointMap[pointI++] = a;
 982             }
 983
 984             label b = tri[1];
 985             if (!pointHad[b])
 986             {
 987                 pointHad[b] = true;
 988                 pointMap[pointI++] = b;
 989             }
 990
 991             label c = tri[2];
 992             if (!pointHad[c])
 993             {
 994                 pointHad[c] = true;
 995                 pointMap[pointI++] = c;
 996             }
 997         }
 998     }
 999
1000     // Trim
1001     faceMap.setSize(faceI);
1002
1003     pointMap.setSize(pointI);
1004 }
1005
1006
1007 Foam::triSurface Foam::triSurface::subsetMesh
1008 (
1009     const boolList& include,
1010     labelList& pointMap,
1011     labelList& faceMap
1012 ) const
1013 {
1014     const pointField& locPoints = localPoints();
1015     const List<labelledTri>& locFaces = localFaces();
1016
1017     // Fill pointMap, faceMap
1018     subsetMeshMap(include, pointMap, faceMap);
1019
1020
1021     // Create compact coordinate list and forward mapping array
1022     pointField newPoints(pointMap.size());
1023     labelList oldToNew(locPoints.size());
1024     forAll(pointMap, pointi)
1025     {
1026         newPoints[pointi] = locPoints[pointMap[pointi]];
1027         oldToNew[pointMap[pointi]] = pointi;
1028     }
1029
1030     // Renumber triangle node labels and compact
1031     List<labelledTri> newTriangles(faceMap.size());
1032
1033     forAll(faceMap, facei)
1034     {
1035         // Get old vertex labels
1036         const labelledTri& tri = locFaces[faceMap[facei]];
1037
1038         newTriangles[facei][0] = oldToNew[tri[0]];
1039         newTriangles[facei][1] = oldToNew[tri[1]];
1040         newTriangles[facei][2] = oldToNew[tri[2]];
1041         newTriangles[facei].region() = tri.region();
1042     }
1043
1044     // Construct subsurface
1045     return triSurface(newTriangles, patches(), newPoints, true);
1046 }
1047
1048
1049 void Foam::triSurface::write
1050 (
1051     const fileName& name,
1052     const bool sortByRegion
1053 ) const
1054 {
1055     write(name, name.ext(), sortByRegion);
1056 }
1057
1058
1059 void Foam::triSurface::write(Ostream& os) const
1060 {
1061     os  << patches() << endl;
1062
1063     //Note: Write with global point numbering
1064     os  << points() << nl
1065         << static_cast<const List<labelledTri>&>(*this) << endl;
1066
1067     // Check state of Ostream
1068     os.check("triSurface::write(Ostream&)");
1069 }
1070
1071
1072 void Foam::triSurface::write(const Time& d) const
1073 {
1074     fileName foamFile(d.caseName() + ".ftr");
1075
1076     fileName foamPath(d.path()/triSurfInstance(d)/typeName/foamFile);
1077
1078     OFstream foamStream(foamPath);
1079
1080     write(foamStream);
1081 }
1082
1083
1084 void Foam::triSurface::writeStats(Ostream& os) const
1085 {
1086     // Unfortunately nPoints constructs meshPoints() so do compact version
1087     // ourselves.
1088     PackedBoolList pointIsUsed(points().size());
1089
1090     label nPoints = 0;
1091     boundBox bb = boundBox::invertedBox;
1092
1093     forAll(*this, triI)
1094     {
1095         const labelledTri& f = operator[](triI);
1096
1097         forAll(f, fp)
1098         {
1099             label pointI = f[fp];
1100             if (pointIsUsed.set(pointI, 1))
1101             {
1102                 bb.min() = ::Foam::min(bb.min(), points()[pointI]);
1103                 bb.max() = ::Foam::max(bb.max(), points()[pointI]);
1104                 nPoints++;
1105             }
1106         }
1107     }
1108
1109     os  << "Triangles    : " << size() << endl
1110         << "Vertices     : " << nPoints << endl
1111         << "Bounding Box : " << bb << endl;
1112 }
1113
1114
1115 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
1116
1117 void Foam::triSurface::operator=(const triSurface& ts)
1118 {
1119     List<labelledTri>::operator=(ts);
1120     clearOut();
1121     storedPoints() = ts.points();
1122     patches_ = ts.patches();
1123 }
1124
1125
1126 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
1127
1128 Foam::Ostream& Foam::operator<<(Ostream& os, const triSurface& sm)
1129 {
1130     sm.write(os);
1131     return os;
1132 }
1133
1134
1135 // ************************************************************************* //