FreeFOAM The Cross-Platform CFD Toolkit
collapseBase.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 "collapseBase.H"
28 #include <OpenFOAM/argList.H>
29 #include <OpenFOAM/OFstream.H>
30 #include <OpenFOAM/SubList.H>
31 #include <OpenFOAM/labelPair.H>
32 #include <meshTools/meshTools.H>
33 #include <OpenFOAM/OSspecific.H>
34 
35 using namespace Foam;
36 
37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 
39 // Dump collapse region to .obj file
40 static void writeRegionOBJ
41 (
42  const triSurface& surf,
43  const label regionI,
44  const labelList& collapseRegion,
45  const labelList& outsideVerts
46 )
47 {
48  fileName dir("regions");
49 
50  mkDir(dir);
51  fileName regionName(dir / "region_" + name(regionI) + ".obj");
52 
53  Pout<< "Dumping region " << regionI << " to file " << regionName << endl;
54 
55  boolList include(surf.size(), false);
56 
57  forAll(collapseRegion, faceI)
58  {
59  if (collapseRegion[faceI] == regionI)
60  {
61  include[faceI] = true;
62  }
63  }
64 
65  labelList pointMap, faceMap;
66 
67  triSurface regionSurf(surf.subsetMesh(include, pointMap, faceMap));
68 
69  //Pout<< "Region " << regionI << " surface:" << nl;
70  //regionSurf.writeStats(Pout);
71 
72  regionSurf.write(regionName);
73 
74 
75  // Dump corresponding outside vertices.
76  fileName pointsName(dir / "regionPoints_" + name(regionI) + ".obj");
77 
78  Pout<< "Dumping region " << regionI << " points to file " << pointsName
79  << endl;
80 
81  OFstream str(pointsName);
82 
83  forAll(outsideVerts, i)
84  {
85  meshTools::writeOBJ(str, surf.localPoints()[outsideVerts[i]]);
86  }
87 }
88 
89 
90 // Split triangle into multiple triangles because edge e being split
91 // into multiple edges.
92 static void splitTri
93 (
94  const labelledTri& f,
95  const edge& e,
96  const labelList& splitPoints,
98 )
99 {
100  label oldNTris = tris.size();
101 
102  label fp = findIndex(f, e[0]);
103  label fp1 = (fp+1)%3;
104  label fp2 = (fp1+1)%3;
105 
106  if (f[fp1] == e[1])
107  {
108  // Split triangle along fp to fp1
109  tris.append(labelledTri(f[fp2], f[fp], splitPoints[0], f.region()));
110 
111  for (label i = 1; i < splitPoints.size(); i++)
112  {
113  tris.append
114  (
116  (
117  f[fp2],
118  splitPoints[i-1],
119  splitPoints[i],
120  f.region()
121  )
122  );
123  }
124 
125  tris.append
126  (
128  (
129  f[fp2],
130  splitPoints[splitPoints.size()-1],
131  f[fp1],
132  f.region()
133  )
134  );
135  }
136  else if (f[fp2] == e[1])
137  {
138  // Split triangle along fp2 to fp. (Reverse order of splitPoints)
139 
140  tris.append
141  (
143  (
144  f[fp1],
145  f[fp2],
146  splitPoints[splitPoints.size()-1],
147  f.region()
148  )
149  );
150 
151  for (label i = splitPoints.size()-1; i > 0; --i)
152  {
153  tris.append
154  (
156  (
157  f[fp1],
158  splitPoints[i],
159  splitPoints[i-1],
160  f.region()
161  )
162  );
163  }
164 
165  tris.append
166  (
168  (
169  f[fp1],
170  splitPoints[0],
171  f[fp],
172  f.region()
173  )
174  );
175  }
176  else
177  {
178  FatalErrorIn("splitTri")
179  << "Edge " << e << " not part of triangle " << f
180  << " fp:" << fp
181  << " fp1:" << fp1
182  << " fp2:" << fp2
183  << abort(FatalError);
184  }
185 
186  Pout<< "Split face " << f << " along edge " << e
187  << " into triangles:" << endl;
188 
189  for (label i = oldNTris; i < tris.size(); i++)
190  {
191  Pout<< " " << tris[i] << nl;
192  }
193 }
194 
195 
196 // Insert scalar into sortedVerts/sortedWeights so the weights are in
197 // incrementing order.
198 static bool insertSorted
199 (
200  const label vertI,
201  const scalar weight,
202 
203  labelList& sortedVerts,
204  scalarField& sortedWeights
205 )
206 {
207  if (findIndex(sortedVerts, vertI) != -1)
208  {
209  FatalErrorIn("insertSorted") << "Trying to insert vertex " << vertI
210  << " which is already in list of sorted vertices "
211  << sortedVerts << abort(FatalError);
212  }
213 
214  if (weight <= 0 || weight >= 1)
215  {
216  FatalErrorIn("insertSorted") << "Trying to insert vertex " << vertI
217  << " with illegal weight " << weight
218  << " into list of sorted vertices "
219  << sortedVerts << abort(FatalError);
220  }
221 
222 
223  label insertI = sortedVerts.size();
224 
225  forAll(sortedVerts, sortedI)
226  {
227  scalar w = sortedWeights[sortedI];
228 
229  if (mag(w - weight) < SMALL)
230  {
231  WarningIn("insertSorted")
232  << "Trying to insert weight " << weight << " which is close to"
233  << " existing weight " << w << " in " << sortedWeights
234  << endl;
235  }
236 
237  if (w > weight)
238  {
239  // Start inserting before sortedI.
240  insertI = sortedI;
241  break;
242  }
243  }
244 
245 
246  label sz = sortedWeights.size();
247 
248  sortedWeights.setSize(sz + 1);
249  sortedVerts.setSize(sz + 1);
250 
251  // Leave everything up to (not including) insertI intact.
252 
253  // Make space by copying from insertI up.
254  for (label i = sz-1; i >= insertI; --i)
255  {
256  sortedWeights[i+1] = sortedWeights[i];
257  sortedVerts[i+1] = sortedVerts[i];
258  }
259  sortedWeights[insertI] = weight;
260  sortedVerts[insertI] = vertI;
261 
262  return true;
263 }
264 
265 
266 // Mark all faces that are going to be collapsed.
267 // faceToEdge: per face -1 or the base edge of the face.
268 static void markCollapsedFaces
269 (
270  const triSurface& surf,
271  const scalar minLen,
272  labelList& faceToEdge
273 )
274 {
275  faceToEdge.setSize(surf.size());
276  faceToEdge = -1;
277 
278  const pointField& localPoints = surf.localPoints();
279  const labelListList& edgeFaces = surf.edgeFaces();
280 
281  forAll(edgeFaces, edgeI)
282  {
283  const edge& e = surf.edges()[edgeI];
284 
285  const labelList& eFaces = surf.edgeFaces()[edgeI];
286 
287  forAll(eFaces, i)
288  {
289  label faceI = eFaces[i];
290 
291  // Check distance of vertex to edge.
292  label opposite0 =
294  (
295  surf,
296  faceI,
297  edgeI
298  );
299 
300  pointHit pHit =
301  e.line(localPoints).nearestDist
302  (
303  localPoints[opposite0]
304  );
305 
306  if (pHit.hit() && pHit.distance() < minLen)
307  {
308  // Remove faceI and split all other faces using this
309  // edge. This is done by 'replacing' the edgeI with the
310  // opposite0 vertex
311  Pout<< "Splitting face " << faceI << " since distance "
312  << pHit.distance()
313  << " from vertex " << opposite0
314  << " to edge " << edgeI
315  << " points "
316  << localPoints[e[0]]
317  << localPoints[e[1]]
318  << " is too small" << endl;
319 
320  // Mark face as being collapsed
321  if (faceToEdge[faceI] != -1)
322  {
323  FatalErrorIn("markCollapsedFaces")
324  << "Cannot collapse face " << faceI << " since "
325  << " is marked to be collapsed both to edge "
326  << faceToEdge[faceI] << " and " << edgeI
327  << abort(FatalError);
328  }
329 
330  faceToEdge[faceI] = edgeI;
331  }
332  }
333  }
334 }
335 
336 
337 // Recurse through collapsed faces marking all of them with regionI (in
338 // collapseRegion)
339 static void markRegion
340 (
341  const triSurface& surf,
342  const labelList& faceToEdge,
343  const label regionI,
344  const label faceI,
345  labelList& collapseRegion
346 )
347 {
348  if (faceToEdge[faceI] == -1 || collapseRegion[faceI] != -1)
349  {
350  FatalErrorIn("markRegion")
351  << "Problem : crossed into uncollapsed/regionized face"
352  << abort(FatalError);
353  }
354 
355  collapseRegion[faceI] = regionI;
356 
357  // Recurse across edges to collapsed neighbours
358 
359  const labelList& fEdges = surf.faceEdges()[faceI];
360 
361  forAll(fEdges, fEdgeI)
362  {
363  label edgeI = fEdges[fEdgeI];
364 
365  const labelList& eFaces = surf.edgeFaces()[edgeI];
366 
367  forAll(eFaces, i)
368  {
369  label nbrFaceI = eFaces[i];
370 
371  if (faceToEdge[nbrFaceI] != -1)
372  {
373  if (collapseRegion[nbrFaceI] == -1)
374  {
375  markRegion
376  (
377  surf,
378  faceToEdge,
379  regionI,
380  nbrFaceI,
381  collapseRegion
382  );
383  }
384  else if (collapseRegion[nbrFaceI] != regionI)
385  {
386  FatalErrorIn("markRegion")
387  << "Edge:" << edgeI << " between face " << faceI
388  << " with region " << regionI
389  << " and face " << nbrFaceI
390  << " with region " << collapseRegion[nbrFaceI]
391  << endl;
392  }
393  }
394  }
395  }
396 }
397 
398 
399 // Mark every face with region (in collapseRegion) (or -1).
400 // Return number of regions.
401 static label markRegions
402 (
403  const triSurface& surf,
404  const labelList& faceToEdge,
405  labelList& collapseRegion
406 )
407 {
408  label regionI = 0;
409 
410  forAll(faceToEdge, faceI)
411  {
412  if (collapseRegion[faceI] == -1 && faceToEdge[faceI] != -1)
413  {
414  Pout<< "markRegions : Marking region:" << regionI
415  << " starting from face " << faceI << endl;
416 
417  // Collapsed face. Mark connected region with current region number
418  markRegion(surf, faceToEdge, regionI++, faceI, collapseRegion);
419  }
420  }
421  return regionI;
422 }
423 
424 
425 // Type of region.
426 // -1 : edge inbetween uncollapsed faces.
427 // -2 : edge inbetween collapsed faces
428 // >=0 : edge inbetween uncollapsed and collapsed region. Returns region.
429 static label edgeType
430 (
431  const triSurface& surf,
432  const labelList& collapseRegion,
433  const label edgeI
434 )
435 {
436  const labelList& eFaces = surf.edgeFaces()[edgeI];
437 
438  // Detect if edge is inbetween collapseRegion and non-collapse face
439  bool usesUncollapsed = false;
440  label usesRegion = -1;
441 
442  forAll(eFaces, i)
443  {
444  label faceI = eFaces[i];
445 
446  label region = collapseRegion[faceI];
447 
448  if (region == -1)
449  {
450  usesUncollapsed = true;
451  }
452  else if (usesRegion == -1)
453  {
454  usesRegion = region;
455  }
456  else if (usesRegion != region)
457  {
458  FatalErrorIn("edgeRegion") << abort(FatalError);
459  }
460  else
461  {
462  // Equal regions.
463  }
464  }
465 
466  if (usesUncollapsed)
467  {
468  if (usesRegion == -1)
469  {
470  // uncollapsed faces only.
471  return -1;
472  }
473  else
474  {
475  // between uncollapsed and collapsed.
476  return usesRegion;
477  }
478  }
479  else
480  {
481  if (usesRegion == -1)
482  {
483  FatalErrorIn("edgeRegion") << abort(FatalError);
484  return -2;
485  }
486  else
487  {
488  return -2;
489  }
490  }
491 }
492 
493 
494 // Get points on outside edge of region (= outside points)
495 static labelListList getOutsideVerts
496 (
497  const triSurface& surf,
498  const labelList& collapseRegion,
499  const label nRegions
500 )
501 {
502  const labelListList& edgeFaces = surf.edgeFaces();
503 
504  // Per region all the outside vertices.
505  labelListList outsideVerts(nRegions);
506 
507  forAll(edgeFaces, edgeI)
508  {
509  // Detect if edge is inbetween collapseRegion and non-collapse face
510  label regionI = edgeType(surf, collapseRegion, edgeI);
511 
512  if (regionI >= 0)
513  {
514  // Edge borders both uncollapsed face and collapsed face on region
515  // usesRegion.
516 
517  const edge& e = surf.edges()[edgeI];
518 
519  labelList& regionVerts = outsideVerts[regionI];
520 
521  // Add both edge points to regionVerts.
522  forAll(e, eI)
523  {
524  label v = e[eI];
525 
526  if (findIndex(regionVerts, v) == -1)
527  {
528  label sz = regionVerts.size();
529  regionVerts.setSize(sz+1);
530  regionVerts[sz] = v;
531  }
532  }
533  }
534  }
535 
536  return outsideVerts;
537 }
538 
539 
540 // n^2 search for furthest removed point pair.
541 static labelPair getSpanPoints
542 (
543  const triSurface& surf,
544  const labelList& outsideVerts
545 )
546 {
547  const pointField& localPoints = surf.localPoints();
548 
549  scalar maxDist = -GREAT;
550  labelPair maxPair;
551 
552  forAll(outsideVerts, i)
553  {
554  label v0 = outsideVerts[i];
555 
556  for (label j = i+1; j < outsideVerts.size(); j++)
557  {
558  label v1 = outsideVerts[j];
559 
560  scalar d = mag(localPoints[v0] - localPoints[v1]);
561 
562  if (d > maxDist)
563  {
564  maxDist = d;
565  maxPair[0] = v0;
566  maxPair[1] = v1;
567  }
568  }
569  }
570 
571  return maxPair;
572 }
573 
574 
575 // Project all non-span points onto the span edge.
576 static void projectNonSpanPoints
577 (
578  const triSurface& surf,
579  const labelList& outsideVerts,
580  const labelPair& spanPair,
581  labelList& sortedVertices,
582  scalarField& sortedWeights
583 )
584 {
585  const point& p0 = surf.localPoints()[spanPair[0]];
586  const point& p1 = surf.localPoints()[spanPair[1]];
587 
588  forAll(outsideVerts, i)
589  {
590  label v = outsideVerts[i];
591 
592  if (v != spanPair[0] && v != spanPair[1])
593  {
594  // Is a non-span point. Project onto spanning edge.
595 
596  pointHit pHit =
597  linePointRef(p0, p1).nearestDist
598  (
599  surf.localPoints()[v]
600  );
601 
602  if (!pHit.hit())
603  {
604  FatalErrorIn("projectNonSpanPoints")
605  << abort(FatalError);
606  }
607 
608  scalar w = mag(pHit.hitPoint() - p0) / mag(p1 - p0);
609 
610  insertSorted(v, w, sortedVertices, sortedWeights);
611  }
612  }
613 }
614 
615 
616 // Slice part of the orderVertices (and optionally reverse) for this edge.
617 static void getSplitVerts
618 (
619  const triSurface& surf,
620  const label regionI,
621  const labelPair& spanPoints,
622  const labelList& orderedVerts,
623  const scalarField& orderedWeights,
624  const label edgeI,
625 
626  labelList& splitVerts,
627  scalarField& splitWeights
628 )
629 {
630  const edge& e = surf.edges()[edgeI];
631  const label sz = orderedVerts.size();
632 
633  if (e[0] == spanPoints[0])
634  {
635  // Edge in same order as spanPoints&orderedVerts. Keep order.
636 
637  if (e[1] == spanPoints[1])
638  {
639  // Copy all.
640  splitVerts = orderedVerts;
641  splitWeights = orderedWeights;
642  }
643  else
644  {
645  // Copy upto (but not including) e[1]
646  label i1 = findIndex(orderedVerts, e[1]);
647  splitVerts = SubList<label>(orderedVerts, i1, 0);
648  splitWeights = SubList<scalar>(orderedWeights, i1, 0);
649  }
650  }
651  else if (e[0] == spanPoints[1])
652  {
653  // Reverse.
654 
655  if (e[1] == spanPoints[0])
656  {
657  // Copy all.
658  splitVerts = orderedVerts;
659  reverse(splitVerts);
660  splitWeights = orderedWeights;
661  reverse(splitWeights);
662  }
663  else
664  {
665  // Copy downto (but not including) e[1]
666 
667  label i1 = findIndex(orderedVerts, e[1]);
668  splitVerts = SubList<label>(orderedVerts, sz-(i1+1), i1+1);
669  reverse(splitVerts);
670  splitWeights = SubList<scalar>(orderedWeights, sz-(i1+1), i1+1);
671  reverse(splitWeights);
672  }
673  }
674  else if (e[1] == spanPoints[0])
675  {
676  // Reverse.
677 
678  // Copy upto (but not including) e[0]
679 
680  label i0 = findIndex(orderedVerts, e[0]);
681  splitVerts = SubList<label>(orderedVerts, i0, 0);
682  reverse(splitVerts);
683  splitWeights = SubList<scalar>(orderedWeights, i0, 0);
684  reverse(splitWeights);
685  }
686  else if (e[1] == spanPoints[1])
687  {
688  // Copy from (but not including) e[0] to end
689 
690  label i0 = findIndex(orderedVerts, e[0]);
691  splitVerts = SubList<label>(orderedVerts, sz-(i0+1), i0+1);
692  splitWeights = SubList<scalar>(orderedWeights, sz-(i0+1), i0+1);
693  }
694  else
695  {
696  label i0 = findIndex(orderedVerts, e[0]);
697  label i1 = findIndex(orderedVerts, e[1]);
698 
699  if (i0 == -1 || i1 == -1)
700  {
701  FatalErrorIn("getSplitVerts")
702  << "Did not find edge in projected vertices." << nl
703  << "region:" << regionI << nl
704  << "spanPoints:" << spanPoints
705  << " coords:" << surf.localPoints()[spanPoints[0]]
706  << surf.localPoints()[spanPoints[1]] << nl
707  << "edge:" << edgeI
708  << " verts:" << e
709  << " coords:" << surf.localPoints()[e[0]]
710  << surf.localPoints()[e[1]] << nl
711  << "orderedVerts:" << orderedVerts << nl
712  << abort(FatalError);
713  }
714 
715  if (i0 < i1)
716  {
717  splitVerts = SubList<label>(orderedVerts, i1-i0-1, i0+1);
718  splitWeights = SubList<scalar>(orderedWeights, i1-i0-1, i0+1);
719  }
720  else
721  {
722  splitVerts = SubList<label>(orderedVerts, i0-i1-1, i1+1);
723  reverse(splitVerts);
724  splitWeights = SubList<scalar>(orderedWeights, i0-i1-1, i1+1);
725  reverse(splitWeights);
726  }
727  }
728 }
729 
730 
731 label collapseBase(triSurface& surf, const scalar minLen)
732 {
733  label nTotalSplit = 0;
734 
735  label iter = 0;
736 
737  while (true)
738  {
739  // Detect faces to collapse
740  // ~~~~~~~~~~~~~~~~~~~~~~~~
741 
742  // -1 or edge the face is collapsed onto.
743  labelList faceToEdge(surf.size(), -1);
744 
745  // Calculate faceToEdge (face collapses)
746  markCollapsedFaces(surf, minLen, faceToEdge);
747 
748 
749  // Find regions of connected collapsed faces
750  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
751 
752  // per face -1 or region
753  labelList collapseRegion(surf.size(), -1);
754 
755  label nRegions = markRegions(surf, faceToEdge, collapseRegion);
756 
757  Pout<< "Detected " << nRegions << " regions of faces to be collapsed"
758  << nl << endl;
759 
760  // Pick up all vertices on outside of region
761  labelListList outsideVerts
762  (
763  getOutsideVerts(surf, collapseRegion, nRegions)
764  );
765 
766  // For all regions determine maximum distance between points
767  List<labelPair> spanPoints(nRegions);
768  labelListList orderedVertices(nRegions);
769  List<scalarField> orderedWeights(nRegions);
770 
771  forAll(spanPoints, regionI)
772  {
773  spanPoints[regionI] = getSpanPoints(surf, outsideVerts[regionI]);
774 
775  Pout<< "For region " << regionI << " found extrema at points "
776  << surf.localPoints()[spanPoints[regionI][0]]
777  << surf.localPoints()[spanPoints[regionI][1]]
778  << endl;
779 
780  // Project all non-span points onto the span edge.
781  projectNonSpanPoints
782  (
783  surf,
784  outsideVerts[regionI],
785  spanPoints[regionI],
786  orderedVertices[regionI],
787  orderedWeights[regionI]
788  );
789 
790  Pout<< "For region:" << regionI
791  << " span:" << spanPoints[regionI]
792  << " orderedVerts:" << orderedVertices[regionI]
793  << " orderedWeights:" << orderedWeights[regionI]
794  << endl;
795 
796  writeRegionOBJ
797  (
798  surf,
799  regionI,
800  collapseRegion,
801  outsideVerts[regionI]
802  );
803 
804  Pout<< endl;
805  }
806 
807 
808 
809  // Actually split the edges
810  // ~~~~~~~~~~~~~~~~~~~~~~~~
811 
812 
813  const List<labelledTri>& localFaces = surf.localFaces();
814  const edgeList& edges = surf.edges();
815 
816  label nSplit = 0;
817 
818  // Storage for new triangles.
819  DynamicList<labelledTri> newTris(surf.size());
820 
821  // Whether face has been dealt with (either copied/split or deleted)
822  boolList faceHandled(surf.size(), false);
823 
824 
825  forAll(edges, edgeI)
826  {
827  const edge& e = edges[edgeI];
828 
829  // Detect if edge is inbetween collapseRegion and non-collapse face
830  label regionI = edgeType(surf, collapseRegion, edgeI);
831 
832  if (regionI == -2)
833  {
834  // inbetween collapsed faces. nothing needs to be done.
835  }
836  else if (regionI == -1)
837  {
838  // edge inbetween uncollapsed faces. Handle these later on.
839  }
840  else
841  {
842  // some faces around edge are collapsed.
843 
844  // Find additional set of points on edge to be used to split
845  // the remaining faces.
846 
847  labelList splitVerts;
848  scalarField splitWeights;
849  getSplitVerts
850  (
851  surf,
852  regionI,
853  spanPoints[regionI],
854  orderedVertices[regionI],
855  orderedWeights[regionI],
856  edgeI,
857 
858  splitVerts,
859  splitWeights
860  );
861 
862  if (splitVerts.size())
863  {
864  // Split edge using splitVerts. All non-collapsed triangles
865  // using edge will get split.
866 
867 
868  {
869  const pointField& localPoints = surf.localPoints();
870  Pout<< "edge " << edgeI << ' ' << e
871  << " points "
872  << localPoints[e[0]] << ' ' << localPoints[e[1]]
873  << " split into edges with extra points:"
874  << endl;
875  forAll(splitVerts, i)
876  {
877  Pout<< " " << splitVerts[i] << " weight "
878  << splitWeights[i] << nl;
879  }
880  }
881 
882  const labelList& eFaces = surf.edgeFaces()[edgeI];
883 
884  forAll(eFaces, i)
885  {
886  label faceI = eFaces[i];
887 
888  if (!faceHandled[faceI] && faceToEdge[faceI] == -1)
889  {
890  // Split face to use vertices.
891  splitTri
892  (
893  localFaces[faceI],
894  e,
895  splitVerts,
896  newTris
897  );
898 
899  faceHandled[faceI] = true;
900 
901  nSplit++;
902  }
903  }
904  }
905  }
906  }
907 
908  // Copy all unsplit faces
909  forAll(faceHandled, faceI)
910  {
911  if (!faceHandled[faceI] && faceToEdge[faceI] == -1)
912  {
913  newTris.append(localFaces[faceI]);
914  }
915  }
916 
917  Pout<< "collapseBase : splitting " << nSplit << " triangles"
918  << endl;
919 
920  nTotalSplit += nSplit;
921 
922  if (nSplit == 0)
923  {
924  break;
925  }
926 
927  // Pack the triangles
928  newTris.shrink();
929 
930  Pout<< "Resetting surface from " << surf.size() << " to "
931  << newTris.size() << " triangles" << endl;
932  surf = triSurface(newTris, surf.patches(), surf.localPoints());
933 
934  {
935  fileName fName("bla" + name(iter) + ".obj");
936  Pout<< "Writing surf to " << fName << endl;
937  surf.write(fName);
938  }
939 
940  iter++;
941  }
942 
943  // Remove any unused vertices
944  surf = triSurface(surf.localFaces(), surf.patches(), surf.localPoints());
945 
946  return nTotalSplit;
947 }
948 
949 
950 // ************************ vim: set sw=4 sts=4 et: ************************ //