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