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