FreeFOAM The Cross-Platform CFD Toolkit
isoSurfaceCell.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 "isoSurfaceCell.H"
27 #include <OpenFOAM/dictionary.H>
28 #include <OpenFOAM/polyMesh.H>
29 #include <OpenFOAM/mergePoints.H>
30 #include <OpenFOAM/tetMatcher.H>
31 #include <OpenFOAM/syncTools.H>
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39  defineTypeNameAndDebug(isoSurfaceCell, 0);
40 }
41 
42 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
43 
44 Foam::scalar Foam::isoSurfaceCell::isoFraction
45 (
46  const scalar s0,
47  const scalar s1
48 ) const
49 {
50  scalar d = s1-s0;
51 
52  if (mag(d) > VSMALL)
53  {
54  return (iso_-s0)/d;
55  }
56  else
57  {
58  return -1.0;
59  }
60 }
61 
62 
63 //Foam::List<Foam::triFace> Foam::isoSurfaceCell::triangulate(const face& f)
64 // const
65 //{
66 // faceList triFaces(f.nTriangles(mesh_.points()));
67 // label triFaceI = 0;
68 // f.triangles(mesh_.points(), triFaceI, triFaces);
69 //
70 // List<triFace> tris(triFaces.size());
71 // forAll(triFaces, i)
72 // {
73 // tris[i][0] = triFaces[i][0];
74 // tris[i][1] = triFaces[i][1];
75 // tris[i][2] = triFaces[i][2];
76 // }
77 // return tris;
78 //}
79 
80 
81 Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
82 (
83  const PackedBoolList& isTet,
84  const scalarField& cellValues,
85  const scalarField& pointValues,
86  const label cellI
87 ) const
88 {
89  const cell& cFaces = mesh_.cells()[cellI];
90 
91  if (isTet.get(cellI) == 1)
92  {
93  forAll(cFaces, cFaceI)
94  {
95  const face& f = mesh_.faces()[cFaces[cFaceI]];
96 
97  for (label fp = 1; fp < f.size() - 1; fp++)
98  {
99  triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
100 
101  bool aLower = (pointValues[tri[0]] < iso_);
102  bool bLower = (pointValues[tri[1]] < iso_);
103  bool cLower = (pointValues[tri[2]] < iso_);
104 
105  if (aLower == bLower && aLower == cLower)
106  {}
107  else
108  {
109  return CUT;
110  }
111  }
112  }
113  return NOTCUT;
114  }
115  else
116  {
117  bool cellLower = (cellValues[cellI] < iso_);
118 
119  // First check if there is any cut in cell
120  bool edgeCut = false;
121 
122  forAll(cFaces, cFaceI)
123  {
124  const face& f = mesh_.faces()[cFaces[cFaceI]];
125 
126  // Check pyramids cut
127  forAll(f, fp)
128  {
129  if ((pointValues[f[fp]] < iso_) != cellLower)
130  {
131  edgeCut = true;
132  break;
133  }
134  }
135 
136  if (edgeCut)
137  {
138  break;
139  }
140 
141  for (label fp = 1; fp < f.size() - 1; fp++)
142  {
143  triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
144  //List<triFace> tris(triangulate(f));
145  //forAll(tris, i)
146  //{
147  // const triFace& tri = tris[i];
148 
149  bool aLower = (pointValues[tri[0]] < iso_);
150  bool bLower = (pointValues[tri[1]] < iso_);
151  bool cLower = (pointValues[tri[2]] < iso_);
152 
153  if (aLower == bLower && aLower == cLower)
154  {}
155  else
156  {
157  edgeCut = true;
158  break;
159  }
160  }
161 
162  if (edgeCut)
163  {
164  break;
165  }
166  }
167 
168  if (edgeCut)
169  {
170  // Count actual cuts (expensive since addressing needed)
171  // Note: not needed if you don't want to preserve maxima/minima
172  // centred around cellcentre. In that case just always return CUT
173 
174  const labelList& cPoints = mesh_.cellPoints(cellI);
175 
176  label nPyrCuts = 0;
177 
178  forAll(cPoints, i)
179  {
180  if ((pointValues[cPoints[i]] < iso_) != cellLower)
181  {
182  nPyrCuts++;
183  }
184  }
185 
186  if (nPyrCuts == cPoints.size())
187  {
188  return SPHERE;
189  }
190  else
191  {
192  return CUT;
193  }
194  }
195  else
196  {
197  return NOTCUT;
198  }
199  }
200 }
201 
202 
203 void Foam::isoSurfaceCell::calcCutTypes
204 (
205  const PackedBoolList& isTet,
206  const scalarField& cVals,
207  const scalarField& pVals
208 )
209 {
210  cellCutType_.setSize(mesh_.nCells());
211  nCutCells_ = 0;
212  forAll(mesh_.cells(), cellI)
213  {
214  cellCutType_[cellI] = calcCutType(isTet, cVals, pVals, cellI);
215 
216  if (cellCutType_[cellI] == CUT)
217  {
218  nCutCells_++;
219  }
220  }
221 
222  if (debug)
223  {
224  Pout<< "isoSurfaceCell : detected " << nCutCells_
225  << " candidate cut cells." << endl;
226  }
227 }
228 
229 
230 
231 // Return the two common points between two triangles
232 Foam::labelPair Foam::isoSurfaceCell::findCommonPoints
233 (
234  const labelledTri& tri0,
235  const labelledTri& tri1
236 )
237 {
238  labelPair common(-1, -1);
239 
240  label fp0 = 0;
241  label fp1 = findIndex(tri1, tri0[fp0]);
242 
243  if (fp1 == -1)
244  {
245  fp0 = 1;
246  fp1 = findIndex(tri1, tri0[fp0]);
247  }
248 
249  if (fp1 != -1)
250  {
251  // So tri0[fp0] is tri1[fp1]
252 
253  // Find next common point
254  label fp0p1 = tri0.fcIndex(fp0);
255  label fp1p1 = tri1.fcIndex(fp1);
256  label fp1m1 = tri1.rcIndex(fp1);
257 
258  if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
259  {
260  common[0] = tri0[fp0];
261  common[1] = tri0[fp0p1];
262  }
263  }
264  return common;
265 }
266 
267 
268 // Caculate centre of surface.
269 Foam::point Foam::isoSurfaceCell::calcCentre(const triSurface& s)
270 {
272 
273  forAll(s, i)
274  {
275  sum += s[i].centre(s.points());
276  }
277  return sum/s.size();
278 }
279 
280 
281 // Replace surface (localPoints, localTris) with single point. Returns
282 // point. Destructs arguments.
283 Foam::pointIndexHit Foam::isoSurfaceCell::collapseSurface
284 (
285  pointField& localPoints,
286  DynamicList<labelledTri, 64>& localTris
287 )
288 {
289  pointIndexHit info(false, vector::zero, localTris.size());
290 
291  if (localTris.size() == 1)
292  {
293  const labelledTri& tri = localTris[0];
294  info.setPoint(tri.centre(localPoints));
295  info.setHit();
296  }
297  else if (localTris.size() == 2)
298  {
299  // Check if the two triangles share an edge.
300  const labelledTri& tri0 = localTris[0];
301  const labelledTri& tri1 = localTris[0];
302 
303  labelPair shared = findCommonPoints(tri0, tri1);
304 
305  if (shared[0] != -1)
306  {
307  info.setPoint
308  (
309  0.5
310  * (
311  tri0.centre(localPoints)
312  + tri1.centre(localPoints)
313  )
314  );
315  info.setHit();
316  }
317  }
318  else if (localTris.size())
319  {
320  // Check if single region. Rare situation.
321  triSurface surf
322  (
323  localTris,
325  localPoints,
326  true
327  );
328  localTris.clearStorage();
329 
330  labelList faceZone;
331  label nZones = surf.markZones
332  (
333  boolList(surf.nEdges(), false),
334  faceZone
335  );
336 
337  if (nZones == 1)
338  {
339  info.setPoint(calcCentre(surf));
340  info.setHit();
341  }
342  }
343 
344  return info;
345 }
346 
347 
348 void Foam::isoSurfaceCell::calcSnappedCc
349 (
350  const PackedBoolList& isTet,
351  const scalarField& cVals,
352  const scalarField& pVals,
353 
354  DynamicList<point>& snappedPoints,
355  labelList& snappedCc
356 ) const
357 {
358  const pointField& cc = mesh_.cellCentres();
359  const pointField& pts = mesh_.points();
360 
361  snappedCc.setSize(mesh_.nCells());
362  snappedCc = -1;
363 
364  // Work arrays
365  DynamicList<point, 64> localPoints(64);
366  DynamicList<labelledTri, 64> localTris(64);
367  Map<label> pointToLocal(64);
368 
369  forAll(mesh_.cells(), cellI)
370  {
371  if (cellCutType_[cellI] == CUT && isTet.get(cellI) == 0)
372  {
373  scalar cVal = cVals[cellI];
374 
375  const cell& cFaces = mesh_.cells()[cellI];
376 
377  localPoints.clear();
378  localTris.clear();
379  pointToLocal.clear();
380 
381  // Create points for all intersections close to cell centre
382  // (i.e. from pyramid edges)
383 
384  forAll(cFaces, cFaceI)
385  {
386  const face& f = mesh_.faces()[cFaces[cFaceI]];
387 
388  forAll(f, fp)
389  {
390  label pointI = f[fp];
391 
392  scalar s = isoFraction(cVal, pVals[pointI]);
393 
394  if (s >= 0.0 && s <= 0.5)
395  {
396  if (pointToLocal.insert(pointI, localPoints.size()))
397  {
398  localPoints.append((1.0-s)*cc[cellI]+s*pts[pointI]);
399  }
400  }
401  }
402  }
403 
404  if (localPoints.size() == 1)
405  {
406  // No need for any analysis.
407  snappedCc[cellI] = snappedPoints.size();
408  snappedPoints.append(localPoints[0]);
409 
410  //Pout<< "cell:" << cellI
411  // << " at " << mesh_.cellCentres()[cellI]
412  // << " collapsing " << localPoints
413  // << " intersections down to "
414  // << snappedPoints[snappedCc[cellI]] << endl;
415  }
416  else if (localPoints.size() == 2)
417  {
418  //? No need for any analysis.???
419  snappedCc[cellI] = snappedPoints.size();
420  snappedPoints.append(0.5*(localPoints[0]+localPoints[1]));
421 
422  //Pout<< "cell:" << cellI
423  // << " at " << mesh_.cellCentres()[cellI]
424  // << " collapsing " << localPoints
425  // << " intersections down to "
426  // << snappedPoints[snappedCc[cellI]] << endl;
427  }
428  else if (localPoints.size())
429  {
430  // Need to analyse
431  forAll(cFaces, cFaceI)
432  {
433  const face& f = mesh_.faces()[cFaces[cFaceI]];
434 
435  // Do a tetrahedrisation. Each face to cc becomes pyr.
436  // Each pyr gets split into tets by diagonalisation
437  // of face.
438 
439  for (label fp = 1; fp < f.size() - 1; fp++)
440  {
441  triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
442  //List<triFace> tris(triangulate(f));
443  //forAll(tris, i)
444  //{
445  // const triFace& tri = tris[i];
446 
447  // Get fractions for the three edges to cell centre
448  FixedList<scalar, 3> s(3);
449  s[0] = isoFraction(cVal, pVals[tri[0]]);
450  s[1] = isoFraction(cVal, pVals[tri[1]]);
451  s[2] = isoFraction(cVal, pVals[tri[2]]);
452 
453  if
454  (
455  (s[0] >= 0.0 && s[0] <= 0.5)
456  && (s[1] >= 0.0 && s[1] <= 0.5)
457  && (s[2] >= 0.0 && s[2] <= 0.5)
458  )
459  {
460  localTris.append
461  (
462  labelledTri
463  (
464  pointToLocal[tri[0]],
465  pointToLocal[tri[1]],
466  pointToLocal[tri[2]],
467  0
468  )
469  );
470  }
471  }
472  }
473 
474  pointField surfPoints;
475  surfPoints.transfer(localPoints);
476  pointIndexHit info = collapseSurface(surfPoints, localTris);
477 
478  if (info.hit())
479  {
480  snappedCc[cellI] = snappedPoints.size();
481  snappedPoints.append(info.hitPoint());
482 
483  //Pout<< "cell:" << cellI
484  // << " at " << mesh_.cellCentres()[cellI]
485  // << " collapsing " << surfPoints
486  // << " intersections down to "
487  // << snappedPoints[snappedCc[cellI]] << endl;
488  }
489  }
490  }
491  }
492 }
493 
494 
495 // Generate triangles for face connected to pointI
496 void Foam::isoSurfaceCell::genPointTris
497 (
498  const scalarField& cellValues,
499  const scalarField& pointValues,
500  const label pointI,
501  const face& f,
502  const label cellI,
503  DynamicList<point, 64>& localTriPoints
504 ) const
505 {
506  const pointField& cc = mesh_.cellCentres();
507  const pointField& pts = mesh_.points();
508 
509  for (label fp = 1; fp < f.size() - 1; fp++)
510  {
511  triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
512  //List<triFace> tris(triangulate(f));
513  //forAll(tris, i)
514  //{
515  // const triFace& tri = tris[i];
516 
517  label index = findIndex(tri, pointI);
518 
519  if (index == -1)
520  {
521  continue;
522  }
523 
524  // Tet between index..index-1, index..index+1, index..cc
525  label b = tri[tri.fcIndex(index)];
526  label c = tri[tri.rcIndex(index)];
527 
528  // Get fractions for the three edges emanating from point
529  FixedList<scalar, 3> s(3);
530  s[0] = isoFraction(pointValues[pointI], pointValues[b]);
531  s[1] = isoFraction(pointValues[pointI], pointValues[c]);
532  s[2] = isoFraction(pointValues[pointI], cellValues[cellI]);
533 
534  if
535  (
536  (s[0] >= 0.0 && s[0] <= 0.5)
537  && (s[1] >= 0.0 && s[1] <= 0.5)
538  && (s[2] >= 0.0 && s[2] <= 0.5)
539  )
540  {
541  localTriPoints.append((1.0-s[0])*pts[pointI] + s[0]*pts[b]);
542  localTriPoints.append((1.0-s[1])*pts[pointI] + s[1]*pts[c]);
543  localTriPoints.append((1.0-s[2])*pts[pointI] + s[2]*cc[cellI]);
544  }
545  }
546 }
547 
548 
549 // Generate triangle for tet connected to pointI
550 void Foam::isoSurfaceCell::genPointTris
551 (
552  const scalarField& pointValues,
553  const label pointI,
554  const label cellI,
555  DynamicList<point, 64>& localTriPoints
556 ) const
557 {
558  const pointField& pts = mesh_.points();
559  const cell& cFaces = mesh_.cells()[cellI];
560 
561  FixedList<label, 4> tet;
562 
563  label face0 = cFaces[0];
564  const face& f0 = mesh_.faces()[face0];
565 
566  if (mesh_.faceOwner()[face0] != cellI)
567  {
568  tet[0] = f0[0];
569  tet[1] = f0[1];
570  tet[2] = f0[2];
571  }
572  else
573  {
574  tet[0] = f0[0];
575  tet[1] = f0[2];
576  tet[2] = f0[1];
577  }
578 
579  // Find the point on the next face that is not on f0
580  const face& f1 = mesh_.faces()[cFaces[1]];
581 
582  forAll(f1, fp)
583  {
584  label p1 = f1[fp];
585 
586  if (p1 != tet[0] && p1 != tet[1] && p1 != tet[2])
587  {
588  tet[3] = p1;
589  break;
590  }
591  }
592 
593 
594  // Get the index of pointI
595 
596  label i0 = findIndex(tet, pointI);
597  label i1 = tet.fcIndex(i0);
598  label i2 = tet.fcIndex(i1);
599  label i3 = tet.fcIndex(i2);
600 
601  // Get fractions for the three edges emanating from point
602  FixedList<scalar, 3> s(3);
603  s[0] = isoFraction(pointValues[pointI], pointValues[tet[i1]]);
604  s[1] = isoFraction(pointValues[pointI], pointValues[tet[i2]]);
605  s[2] = isoFraction(pointValues[pointI], pointValues[tet[i3]]);
606 
607  if
608  (
609  (s[0] >= 0.0 && s[0] <= 0.5)
610  && (s[1] >= 0.0 && s[1] <= 0.5)
611  && (s[2] >= 0.0 && s[2] <= 0.5)
612  )
613  {
614  localTriPoints.append((1.0-s[0])*pts[pointI] + s[0]*pts[tet[i1]]);
615  localTriPoints.append((1.0-s[1])*pts[pointI] + s[1]*pts[tet[i2]]);
616  localTriPoints.append((1.0-s[2])*pts[pointI] + s[2]*pts[tet[i3]]);
617  }
618 }
619 
620 
621 void Foam::isoSurfaceCell::calcSnappedPoint
622 (
623  const PackedBoolList& isBoundaryPoint,
624  const PackedBoolList& isTet,
625  const scalarField& cVals,
626  const scalarField& pVals,
627 
628  DynamicList<point>& snappedPoints,
629  labelList& snappedPoint
630 ) const
631 {
632  const point greatPoint(VGREAT, VGREAT, VGREAT);
633  pointField collapsedPoint(mesh_.nPoints(), greatPoint);
634 
635 
636  // Work arrays
637  DynamicList<point, 64> localTriPoints(100);
638  labelHashSet localPointCells(100);
639 
640  forAll(mesh_.pointFaces(), pointI)
641  {
642  if (isBoundaryPoint.get(pointI) == 1)
643  {
644  continue;
645  }
646 
647  const labelList& pFaces = mesh_.pointFaces()[pointI];
648 
649  bool anyCut = false;
650 
651  forAll(pFaces, i)
652  {
653  label faceI = pFaces[i];
654 
655  if
656  (
657  cellCutType_[mesh_.faceOwner()[faceI]] == CUT
658  || (
659  mesh_.isInternalFace(faceI)
660  && cellCutType_[mesh_.faceNeighbour()[faceI]] == CUT
661  )
662  )
663  {
664  anyCut = true;
665  break;
666  }
667  }
668 
669  if (!anyCut)
670  {
671  continue;
672  }
673 
674 
675  localPointCells.clear();
676  localTriPoints.clear();
677 
678  forAll(pFaces, pFaceI)
679  {
680  label faceI = pFaces[pFaceI];
681  const face& f = mesh_.faces()[faceI];
682  label own = mesh_.faceOwner()[faceI];
683 
684  // Triangulate around f[0] on owner side
685  if (isTet.get(own) == 1)
686  {
687  if (localPointCells.insert(own))
688  {
689  genPointTris(pVals, pointI, own, localTriPoints);
690  }
691  }
692  else
693  {
694  genPointTris(cVals, pVals, pointI, f, own, localTriPoints);
695  }
696 
697  if (mesh_.isInternalFace(faceI))
698  {
699  label nei = mesh_.faceNeighbour()[faceI];
700 
701  if (isTet.get(nei) == 1)
702  {
703  if (localPointCells.insert(nei))
704  {
705  genPointTris(pVals, pointI, nei, localTriPoints);
706  }
707  }
708  else
709  {
710  genPointTris(cVals, pVals, pointI, f, nei, localTriPoints);
711  }
712  }
713  }
714 
715  if (localTriPoints.size() == 3)
716  {
717  // Single triangle. No need for any analysis. Average points.
719  points.transfer(localTriPoints);
720  collapsedPoint[pointI] = sum(points)/points.size();
721 
722  //Pout<< " point:" << pointI
723  // << " replacing coord:" << mesh_.points()[pointI]
724  // << " by average:" << collapsedPoint[pointI] << endl;
725  }
726  else if (localTriPoints.size())
727  {
728  // Convert points into triSurface.
729 
730  // Merge points and compact out non-valid triangles
731  labelList triMap; // merged to unmerged triangle
732  labelList triPointReverseMap; // unmerged to merged point
733  triSurface surf
734  (
735  stitchTriPoints
736  (
737  false, // do not check for duplicate tris
738  localTriPoints,
739  triPointReverseMap,
740  triMap
741  )
742  );
743 
744  labelList faceZone;
745  label nZones = surf.markZones
746  (
747  boolList(surf.nEdges(), false),
748  faceZone
749  );
750 
751  if (nZones == 1)
752  {
753  collapsedPoint[pointI] = calcCentre(surf);
754  //Pout<< " point:" << pointI << " nZones:" << nZones
755  // << " replacing coord:" << mesh_.points()[pointI]
756  // << " by average:" << collapsedPoint[pointI] << endl;
757  }
758  }
759  }
760 
762  (
763  mesh_,
764  collapsedPoint,
765  minEqOp<point>(),
766  greatPoint,
767  true // are coordinates so separate
768  );
769 
770  snappedPoint.setSize(mesh_.nPoints());
771  snappedPoint = -1;
772 
773  forAll(collapsedPoint, pointI)
774  {
775  if (collapsedPoint[pointI] != greatPoint)
776  {
777  snappedPoint[pointI] = snappedPoints.size();
778  snappedPoints.append(collapsedPoint[pointI]);
779  }
780  }
781 }
782 
783 
784 
785 
786 Foam::triSurface Foam::isoSurfaceCell::stitchTriPoints
787 (
788  const bool checkDuplicates,
789  const List<point>& triPoints,
790  labelList& triPointReverseMap, // unmerged to merged point
791  labelList& triMap // merged to unmerged triangle
792 ) const
793 {
794  label nTris = triPoints.size()/3;
795 
796  if ((triPoints.size() % 3) != 0)
797  {
798  FatalErrorIn("isoSurfaceCell::stitchTriPoints(..)")
799  << "Problem: number of points " << triPoints.size()
800  << " not a multiple of 3." << abort(FatalError);
801  }
802 
803  pointField newPoints;
805  (
806  triPoints,
807  mergeDistance_,
808  false,
809  triPointReverseMap,
810  newPoints
811  );
812 
813  // Check that enough merged.
814  if (debug)
815  {
816  pointField newNewPoints;
817  labelList oldToNew;
818  bool hasMerged = mergePoints
819  (
820  newPoints,
821  mergeDistance_,
822  true,
823  oldToNew,
824  newNewPoints
825  );
826 
827  if (hasMerged)
828  {
829  FatalErrorIn("isoSurfaceCell::stitchTriPoints(..)")
830  << "Merged points contain duplicates"
831  << " when merging with distance " << mergeDistance_ << endl
832  << "merged:" << newPoints.size() << " re-merged:"
833  << newNewPoints.size()
834  << abort(FatalError);
835  }
836  }
837 
838 
839  List<labelledTri> tris;
840  {
841  DynamicList<labelledTri> dynTris(nTris);
842  label rawPointI = 0;
843  DynamicList<label> newToOldTri(nTris);
844 
845  for (label oldTriI = 0; oldTriI < nTris; oldTriI++)
846  {
847  labelledTri tri
848  (
849  triPointReverseMap[rawPointI],
850  triPointReverseMap[rawPointI+1],
851  triPointReverseMap[rawPointI+2],
852  0
853  );
854  rawPointI += 3;
855 
856  if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
857  {
858  newToOldTri.append(oldTriI);
859  dynTris.append(tri);
860  }
861  }
862 
863  triMap.transfer(newToOldTri);
864  tris.transfer(dynTris);
865  }
866 
867 
868  // Use face centres to determine 'flat hole' situation (see RMT paper).
869  // Two unconnected triangles get connected because (some of) the edges
870  // separating them get collapsed. Below only checks for duplicate triangles,
871  // not non-manifold edge connectivity.
872  if (checkDuplicates)
873  {
874  if (debug)
875  {
876  Pout<< "isoSurfaceCell : merged from " << nTris
877  << " down to " << tris.size() << " triangles." << endl;
878  }
879 
880  pointField centres(tris.size());
881  forAll(tris, triI)
882  {
883  centres[triI] = tris[triI].centre(newPoints);
884  }
885 
886  pointField mergedCentres;
887  labelList oldToMerged;
888  bool hasMerged = mergePoints
889  (
890  centres,
891  mergeDistance_,
892  false,
893  oldToMerged,
894  mergedCentres
895  );
896 
897  if (debug)
898  {
899  Pout<< "isoSurfaceCell : detected "
900  << centres.size()-mergedCentres.size()
901  << " duplicate triangles." << endl;
902  }
903 
904  if (hasMerged)
905  {
906  // Filter out duplicates.
907  label newTriI = 0;
908  DynamicList<label> newToOldTri(tris.size());
909  labelList newToMaster(mergedCentres.size(), -1);
910  forAll(tris, triI)
911  {
912  label mergedI = oldToMerged[triI];
913 
914  if (newToMaster[mergedI] == -1)
915  {
916  newToMaster[mergedI] = triI;
917  newToOldTri.append(triMap[triI]);
918  tris[newTriI++] = tris[triI];
919  }
920  }
921 
922  triMap.transfer(newToOldTri);
923  tris.setSize(newTriI);
924  }
925  }
926 
927  return triSurface(tris, geometricSurfacePatchList(0), newPoints, true);
928 }
929 
930 
931 // Does face use valid vertices?
932 bool Foam::isoSurfaceCell::validTri(const triSurface& surf, const label faceI)
933 {
934  // Simple check on indices ok.
935 
936  const labelledTri& f = surf[faceI];
937 
938  if
939  (
940  (f[0] < 0) || (f[0] >= surf.points().size())
941  || (f[1] < 0) || (f[1] >= surf.points().size())
942  || (f[2] < 0) || (f[2] >= surf.points().size())
943  )
944  {
945  WarningIn("validTri(const triSurface&, const label)")
946  << "triangle " << faceI << " vertices " << f
947  << " uses point indices outside point range 0.."
948  << surf.points().size()-1 << endl;
949 
950  return false;
951  }
952 
953  if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
954  {
955  WarningIn("validTri(const triSurface&, const label)")
956  << "triangle " << faceI
957  << " uses non-unique vertices " << f
958  << endl;
959  return false;
960  }
961 
962  // duplicate triangle check
963 
964  const labelList& fFaces = surf.faceFaces()[faceI];
965 
966  // Check if faceNeighbours use same points as this face.
967  // Note: discards normal information - sides of baffle are merged.
968  forAll(fFaces, i)
969  {
970  label nbrFaceI = fFaces[i];
971 
972  if (nbrFaceI <= faceI)
973  {
974  // lower numbered faces already checked
975  continue;
976  }
977 
978  const labelledTri& nbrF = surf[nbrFaceI];
979 
980  if
981  (
982  ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2]))
983  && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2]))
984  && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
985  )
986  {
987  WarningIn("validTri(const triSurface&, const label)")
988  << "triangle " << faceI << " vertices " << f
989  << " has the same vertices as triangle " << nbrFaceI
990  << " vertices " << nbrF
991  << endl;
992 
993  return false;
994  }
995  }
996  return true;
997 }
998 
999 
1000 void Foam::isoSurfaceCell::calcAddressing
1001 (
1002  const triSurface& surf,
1003  List<FixedList<label, 3> >& faceEdges,
1004  labelList& edgeFace0,
1005  labelList& edgeFace1,
1006  Map<labelList>& edgeFacesRest
1007 ) const
1008 {
1009  const pointField& points = surf.points();
1010 
1011  pointField edgeCentres(3*surf.size());
1012  label edgeI = 0;
1013  forAll(surf, triI)
1014  {
1015  const labelledTri& tri = surf[triI];
1016  edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
1017  edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
1018  edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
1019  }
1020 
1021  pointField mergedCentres;
1022  labelList oldToMerged;
1023  bool hasMerged = mergePoints
1024  (
1025  edgeCentres,
1026  mergeDistance_,
1027  false,
1028  oldToMerged,
1029  mergedCentres
1030  );
1031 
1032  if (debug)
1033  {
1034  Pout<< "isoSurfaceCell : detected "
1035  << mergedCentres.size()
1036  << " edges on " << surf.size() << " triangles." << endl;
1037  }
1038 
1039  if (!hasMerged)
1040  {
1041  if (surf.size() == 1)
1042  {
1043  faceEdges.setSize(1);
1044  faceEdges[0][0] = 0;
1045  faceEdges[0][1] = 1;
1046  faceEdges[0][2] = 2;
1047  edgeFace0.setSize(1);
1048  edgeFace0[0] = 0;
1049  edgeFace1.setSize(1);
1050  edgeFace1[0] = -1;
1051  edgeFacesRest.clear();
1052  }
1053  return;
1054  }
1055 
1056 
1057  // Determine faceEdges
1058  faceEdges.setSize(surf.size());
1059  edgeI = 0;
1060  forAll(surf, triI)
1061  {
1062  faceEdges[triI][0] = oldToMerged[edgeI++];
1063  faceEdges[triI][1] = oldToMerged[edgeI++];
1064  faceEdges[triI][2] = oldToMerged[edgeI++];
1065  }
1066 
1067 
1068  // Determine edgeFaces
1069  edgeFace0.setSize(mergedCentres.size());
1070  edgeFace0 = -1;
1071  edgeFace1.setSize(mergedCentres.size());
1072  edgeFace1 = -1;
1073  edgeFacesRest.clear();
1074 
1075  forAll(oldToMerged, oldEdgeI)
1076  {
1077  label triI = oldEdgeI / 3;
1078  label edgeI = oldToMerged[oldEdgeI];
1079 
1080  if (edgeFace0[edgeI] == -1)
1081  {
1082  edgeFace0[edgeI] = triI;
1083  }
1084  else if (edgeFace1[edgeI] == -1)
1085  {
1086  edgeFace1[edgeI] = triI;
1087  }
1088  else
1089  {
1090  //WarningIn("orientSurface(triSurface&)")
1091  // << "Edge " << edgeI << " with centre " << mergedCentres[edgeI]
1092  // << " used by more than two triangles: " << edgeFace0[edgeI]
1093  // << ", "
1094  // << edgeFace1[edgeI] << " and " << triI << endl;
1095  Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);
1096 
1097  if (iter != edgeFacesRest.end())
1098  {
1099  labelList& eFaces = iter();
1100  label sz = eFaces.size();
1101  eFaces.setSize(sz+1);
1102  eFaces[sz] = triI;
1103  }
1104  else
1105  {
1106  edgeFacesRest.insert(edgeI, labelList(1, triI));
1107  }
1108  }
1109  }
1110 }
1111 
1112 
1113 void Foam::isoSurfaceCell::walkOrientation
1114 (
1115  const triSurface& surf,
1116  const List<FixedList<label, 3> >& faceEdges,
1117  const labelList& edgeFace0,
1118  const labelList& edgeFace1,
1119  const label seedTriI,
1120  labelList& flipState
1121 )
1122 {
1123  // Do walk for consistent orientation.
1124  DynamicList<label> changedFaces(surf.size());
1125 
1126  changedFaces.append(seedTriI);
1127 
1128  while (changedFaces.size())
1129  {
1130  DynamicList<label> newChangedFaces(changedFaces.size());
1131 
1132  forAll(changedFaces, i)
1133  {
1134  label triI = changedFaces[i];
1135  const labelledTri& tri = surf[triI];
1136  const FixedList<label, 3>& fEdges = faceEdges[triI];
1137 
1138  forAll(fEdges, fp)
1139  {
1140  label edgeI = fEdges[fp];
1141 
1142  // my points:
1143  label p0 = tri[fp];
1144  label p1 = tri[tri.fcIndex(fp)];
1145 
1146  label nbrI =
1147  (
1148  edgeFace0[edgeI] != triI
1149  ? edgeFace0[edgeI]
1150  : edgeFace1[edgeI]
1151  );
1152 
1153  if (nbrI != -1 && flipState[nbrI] == -1)
1154  {
1155  const labelledTri& nbrTri = surf[nbrI];
1156 
1157  // nbr points
1158  label nbrFp = findIndex(nbrTri, p0);
1159  label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
1160 
1161  bool sameOrientation = (p1 == nbrP1);
1162 
1163  if (flipState[triI] == 0)
1164  {
1165  flipState[nbrI] = (sameOrientation ? 0 : 1);
1166  }
1167  else
1168  {
1169  flipState[nbrI] = (sameOrientation ? 1 : 0);
1170  }
1171  newChangedFaces.append(nbrI);
1172  }
1173  }
1174  }
1175 
1176  changedFaces.transfer(newChangedFaces);
1177  }
1178 }
1179 
1180 
1181 void Foam::isoSurfaceCell::orientSurface
1182 (
1183  triSurface& surf,
1184  const List<FixedList<label, 3> >& faceEdges,
1185  const labelList& edgeFace0,
1186  const labelList& edgeFace1,
1187  const Map<labelList>& edgeFacesRest
1188 )
1189 {
1190  // -1 : unvisited
1191  // 0 : leave as is
1192  // 1 : flip
1193  labelList flipState(surf.size(), -1);
1194 
1195  label seedTriI = 0;
1196 
1197  while (true)
1198  {
1199  // Find first unvisited triangle
1200  for
1201  (
1202  ;
1203  seedTriI < surf.size() && flipState[seedTriI] != -1;
1204  seedTriI++
1205  )
1206  {}
1207 
1208  if (seedTriI == surf.size())
1209  {
1210  break;
1211  }
1212 
1213  // Note: Determine orientation of seedTriI?
1214  // for now assume it is ok
1215  flipState[seedTriI] = 0;
1216 
1217  walkOrientation
1218  (
1219  surf,
1220  faceEdges,
1221  edgeFace0,
1222  edgeFace1,
1223  seedTriI,
1224  flipState
1225  );
1226  }
1227 
1228  // Do actual flipping
1229  surf.clearOut();
1230  forAll(surf, triI)
1231  {
1232  if (flipState[triI] == 1)
1233  {
1234  labelledTri tri(surf[triI]);
1235 
1236  surf[triI][0] = tri[0];
1237  surf[triI][1] = tri[2];
1238  surf[triI][2] = tri[1];
1239  }
1240  else if (flipState[triI] == -1)
1241  {
1242  FatalErrorIn
1243  (
1244  "isoSurfaceCell::orientSurface(triSurface&, const label)"
1245  ) << "problem" << abort(FatalError);
1246  }
1247  }
1248 }
1249 
1250 
1251 // Checks if triangle is connected through edgeI only.
1252 bool Foam::isoSurfaceCell::danglingTriangle
1253 (
1254  const FixedList<label, 3>& fEdges,
1255  const labelList& edgeFace1
1256 )
1257 {
1258  label nOpen = 0;
1259  forAll(fEdges, i)
1260  {
1261  if (edgeFace1[fEdges[i]] == -1)
1262  {
1263  nOpen++;
1264  }
1265  }
1266 
1267  if (nOpen == 1 || nOpen == 2 || nOpen == 3)
1268  {
1269  return true;
1270  }
1271  else
1272  {
1273  return false;
1274  }
1275 }
1276 
1277 
1278 // Mark triangles to keep. Returns number of dangling triangles.
1279 Foam::label Foam::isoSurfaceCell::markDanglingTriangles
1280 (
1281  const List<FixedList<label, 3> >& faceEdges,
1282  const labelList& edgeFace0,
1283  const labelList& edgeFace1,
1284  const Map<labelList>& edgeFacesRest,
1285  boolList& keepTriangles
1286 )
1287 {
1288  keepTriangles.setSize(faceEdges.size());
1289  keepTriangles = true;
1290 
1291  label nDangling = 0;
1292 
1293  // Remove any dangling triangles
1294  forAllConstIter(Map<labelList>, edgeFacesRest, iter)
1295  {
1296  // These are all the non-manifold edges. Filter out all triangles
1297  // with only one connected edge (= this edge)
1298 
1299  label edgeI = iter.key();
1300  const labelList& otherEdgeFaces = iter();
1301 
1302  // Remove all dangling triangles
1303  if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
1304  {
1305  keepTriangles[edgeFace0[edgeI]] = false;
1306  nDangling++;
1307  }
1308  if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
1309  {
1310  keepTriangles[edgeFace1[edgeI]] = false;
1311  nDangling++;
1312  }
1313  forAll(otherEdgeFaces, i)
1314  {
1315  label triI = otherEdgeFaces[i];
1316  if (danglingTriangle(faceEdges[triI], edgeFace1))
1317  {
1318  keepTriangles[triI] = false;
1319  nDangling++;
1320  }
1321  }
1322  }
1323  return nDangling;
1324 }
1325 
1326 
1327 Foam::triSurface Foam::isoSurfaceCell::subsetMesh
1328 (
1329  const triSurface& s,
1330  const labelList& newToOldFaces,
1331  labelList& oldToNewPoints,
1332  labelList& newToOldPoints
1333 )
1334 {
1335  const boolList include
1336  (
1337  createWithValues<boolList>
1338  (
1339  s.size(),
1340  false,
1341  newToOldFaces,
1342  true
1343  )
1344  );
1345 
1346  newToOldPoints.setSize(s.points().size());
1347  oldToNewPoints.setSize(s.points().size());
1348  oldToNewPoints = -1;
1349  {
1350  label pointI = 0;
1351 
1352  forAll(include, oldFacei)
1353  {
1354  if (include[oldFacei])
1355  {
1356  // Renumber labels for triangle
1357  const labelledTri& tri = s[oldFacei];
1358 
1359  forAll(tri, fp)
1360  {
1361  label oldPointI = tri[fp];
1362 
1363  if (oldToNewPoints[oldPointI] == -1)
1364  {
1365  oldToNewPoints[oldPointI] = pointI;
1366  newToOldPoints[pointI++] = oldPointI;
1367  }
1368  }
1369  }
1370  }
1371  newToOldPoints.setSize(pointI);
1372  }
1373 
1374  // Extract points
1375  pointField newPoints(newToOldPoints.size());
1376  forAll(newToOldPoints, i)
1377  {
1378  newPoints[i] = s.points()[newToOldPoints[i]];
1379  }
1380  // Extract faces
1381  List<labelledTri> newTriangles(newToOldFaces.size());
1382 
1383  forAll(newToOldFaces, i)
1384  {
1385  // Get old vertex labels
1386  const labelledTri& tri = s[newToOldFaces[i]];
1387 
1388  newTriangles[i][0] = oldToNewPoints[tri[0]];
1389  newTriangles[i][1] = oldToNewPoints[tri[1]];
1390  newTriangles[i][2] = oldToNewPoints[tri[2]];
1391  newTriangles[i].region() = tri.region();
1392  }
1393 
1394  // Reuse storage.
1395  return triSurface(newTriangles, s.patches(), newPoints, true);
1396 }
1397 
1398 
1399 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1400 
1403  const polyMesh& mesh,
1404  const scalarField& cVals,
1405  const scalarField& pVals,
1406  const scalar iso,
1407  const bool regularise,
1408  const scalar mergeTol
1409 )
1410 :
1411  mesh_(mesh),
1412  iso_(iso),
1413  mergeDistance_(mergeTol*mesh.bounds().mag())
1414 {
1415  // Determine if cell is tet
1416  PackedBoolList isTet(mesh_.nCells());
1417  {
1418  tetMatcher tet;
1419 
1420  forAll(isTet, cellI)
1421  {
1422  if (tet.isA(mesh_, cellI))
1423  {
1424  isTet.set(cellI, 1);
1425  }
1426  }
1427  }
1428 
1429  // Determine if point is on boundary. Points on boundaries are never
1430  // snapped. Coupled boundaries are handled explicitly so not marked here.
1431  PackedBoolList isBoundaryPoint(mesh_.nPoints());
1432  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1433  forAll(patches, patchI)
1434  {
1435  const polyPatch& pp = patches[patchI];
1436 
1437  if (!pp.coupled())
1438  {
1439  label faceI = pp.start();
1440  forAll(pp, i)
1441  {
1442  const face& f = mesh_.faces()[faceI++];
1443 
1444  forAll(f, fp)
1445  {
1446  isBoundaryPoint.set(f[fp], 1);
1447  }
1448  }
1449  }
1450  }
1451 
1452 
1453  // Determine if any cut through cell
1454  calcCutTypes(isTet, cVals, pVals);
1455 
1456  DynamicList<point> snappedPoints(nCutCells_);
1457 
1458  // Per cc -1 or a point inside snappedPoints.
1459  labelList snappedCc;
1460  if (regularise)
1461  {
1462  calcSnappedCc
1463  (
1464  isTet,
1465  cVals,
1466  pVals,
1467  snappedPoints,
1468  snappedCc
1469  );
1470  }
1471  else
1472  {
1473  snappedCc.setSize(mesh_.nCells());
1474  snappedCc = -1;
1475  }
1476 
1477  if (debug)
1478  {
1479  Pout<< "isoSurfaceCell : shifted " << snappedPoints.size()
1480  << " cell centres to intersection." << endl;
1481  }
1482 
1483  snappedPoints.shrink();
1484  label nCellSnaps = snappedPoints.size();
1485 
1486  // Per point -1 or a point inside snappedPoints.
1487  labelList snappedPoint;
1488  if (regularise)
1489  {
1490  calcSnappedPoint
1491  (
1492  isBoundaryPoint,
1493  isTet,
1494  cVals,
1495  pVals,
1496  snappedPoints,
1497  snappedPoint
1498  );
1499  }
1500  else
1501  {
1502  snappedPoint.setSize(mesh_.nPoints());
1503  snappedPoint = -1;
1504  }
1505 
1506  if (debug)
1507  {
1508  Pout<< "isoSurfaceCell : shifted " << snappedPoints.size()-nCellSnaps
1509  << " vertices to intersection." << endl;
1510  }
1511 
1512 
1513 
1514  DynamicList<point> triPoints(nCutCells_);
1515  DynamicList<label> triMeshCells(nCutCells_);
1516 
1517  generateTriPoints
1518  (
1519  cVals,
1520  pVals,
1521 
1522  mesh_.cellCentres(),
1523  mesh_.points(),
1524 
1525  snappedPoints,
1526  snappedCc,
1527  snappedPoint,
1528 
1529  triPoints,
1530  triMeshCells
1531  );
1532 
1533  if (debug)
1534  {
1535  Pout<< "isoSurfaceCell : generated " << triMeshCells.size()
1536  << " unmerged triangles." << endl;
1537  }
1538 
1539  // Merge points and compact out non-valid triangles
1540  labelList triMap; // merged to unmerged triangle
1542  (
1543  stitchTriPoints
1544  (
1545  true, // check for duplicate tris
1546  triPoints,
1547  triPointMergeMap_, // unmerged to merged point
1548  triMap
1549  )
1550  );
1551 
1552  if (debug)
1553  {
1554  Pout<< "isoSurfaceCell : generated " << triMap.size()
1555  << " merged triangles." << endl;
1556  }
1557 
1558  meshCells_.setSize(triMap.size());
1559  forAll(triMap, i)
1560  {
1561  meshCells_[i] = triMeshCells[triMap[i]];
1562  }
1563 
1564  if (debug)
1565  {
1566  Pout<< "isoSurfaceCell : checking " << size()
1567  << " triangles for validity." << endl;
1568 
1569  forAll(*this, triI)
1570  {
1571  // Copied from surfaceCheck
1572  validTri(*this, triI);
1573  }
1574  }
1575 
1576 
1577 //if (false)
1578 {
1579  List<FixedList<label, 3> > faceEdges;
1580  labelList edgeFace0, edgeFace1;
1581  Map<labelList> edgeFacesRest;
1582 
1583 
1584  while (true)
1585  {
1586  // Calculate addressing
1587  calcAddressing(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
1588 
1589  // See if any dangling triangles
1590  boolList keepTriangles;
1591  label nDangling = markDanglingTriangles
1592  (
1593  faceEdges,
1594  edgeFace0,
1595  edgeFace1,
1596  edgeFacesRest,
1597  keepTriangles
1598  );
1599 
1600  if (debug)
1601  {
1602  Pout<< "isoSurfaceCell : detected " << nDangling
1603  << " dangling triangles." << endl;
1604  }
1605 
1606  if (nDangling == 0)
1607  {
1608  break;
1609  }
1610 
1611  // Create face map (new to old)
1612  labelList subsetTriMap(findIndices(keepTriangles, true));
1613 
1614  labelList subsetPointMap;
1615  labelList reversePointMap;
1617  (
1618  subsetMesh
1619  (
1620  *this,
1621  subsetTriMap,
1622  reversePointMap,
1623  subsetPointMap
1624  )
1625  );
1626  meshCells_ = labelField(meshCells_, subsetTriMap);
1627  inplaceRenumber(reversePointMap, triPointMergeMap_);
1628  }
1629 
1630  orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
1631 }
1632 
1633  //combineCellTriangles();
1634 }
1635 
1636 
1641 //void Foam::isoSurfaceCell::combineCellTriangles()
1642 //{
1643 // if (size())
1644 // {
1645 // DynamicList<labelledTri> newTris(size());
1646 // DynamicList<label> newTriToCell(size());
1647 //
1648 // label startTriI = 0;
1649 //
1650 // DynamicList<labelledTri> tris;
1651 //
1652 // for (label triI = 1; triI <= meshCells_.size(); triI++)
1653 // {
1654 // if
1655 // (
1656 // triI == meshCells_.size()
1657 // || meshCells_[triI] != meshCells_[startTriI]
1658 // )
1659 // {
1660 // label nTris = triI-startTriI;
1661 //
1662 // if (nTris == 1)
1663 // {
1664 // newTris.append(operator[](startTriI));
1665 // newTriToCell.append(meshCells_[startTriI]);
1666 // }
1667 // else
1668 // {
1669 // // Collect from startTriI to triI in a triSurface
1670 // tris.clear();
1671 // for (label i = startTriI; i < triI; i++)
1672 // {
1673 // tris.append(operator[](i));
1674 // }
1675 // triSurface cellTris(tris, patches(), points());
1676 // tris.clear();
1677 //
1678 // // Get outside
1679 // const labelListList& loops = cellTris.edgeLoops();
1680 //
1681 // forAll(loops, i)
1682 // {
1683 // // Do proper triangulation of loop
1684 // face loop(renumber(cellTris.meshPoints(), loops[i]));
1685 //
1686 // faceTriangulation faceTris
1687 // (
1688 // points(),
1689 // loop,
1690 // true
1691 // );
1692 //
1693 // // Copy into newTris
1694 // forAll(faceTris, faceTriI)
1695 // {
1696 // const triFace& tri = faceTris[faceTriI];
1697 //
1698 // newTris.append
1699 // (
1700 // labelledTri
1701 // (
1702 // tri[0],
1703 // tri[1],
1704 // tri[2],
1705 // operator[](startTriI).region()
1706 // )
1707 // );
1708 // newTriToCell.append(meshCells_[startTriI]);
1709 // }
1710 // }
1711 // }
1712 //
1713 // startTriI = triI;
1714 // }
1715 // }
1716 // newTris.shrink();
1717 // newTriToCell.shrink();
1718 //
1719 // // Compact
1720 // pointField newPoints(points().size());
1721 // label newPointI = 0;
1722 // labelList oldToNewPoint(points().size(), -1);
1723 //
1724 // forAll(newTris, i)
1725 // {
1726 // labelledTri& tri = newTris[i];
1727 // forAll(tri, j)
1728 // {
1729 // label pointI = tri[j];
1730 //
1731 // if (oldToNewPoint[pointI] == -1)
1732 // {
1733 // oldToNewPoint[pointI] = newPointI;
1734 // newPoints[newPointI++] = points()[pointI];
1735 // }
1736 // tri[j] = oldToNewPoint[pointI];
1737 // }
1738 // }
1739 // newPoints.setSize(newPointI);
1740 //
1741 // triSurface::operator=
1742 // (
1743 // triSurface
1744 // (
1745 // newTris,
1746 // patches(),
1747 // newPoints,
1748 // true
1749 // )
1750 // );
1751 // meshCells_.transfer(newTriToCell);
1752 // }
1753 //}
1755 
1756 // ************************ vim: set sw=4 sts=4 et: ************************ //