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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2016 OpenFOAM Foundation
9  Copyright (C) 2016-2020 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "isoSurfaceCell.H"
30 #include "isoSurfacePoint.H"
31 #include "dictionary.H"
32 #include "polyMesh.H"
33 #include "mergePoints.H"
34 #include "tetMatcher.H"
35 #include "syncTools.H"
36 #include "triSurface.H"
37 #include "triSurfaceTools.H"
38 #include "Time.H"
39 #include "triPoints.H"
40 
41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
42 
43 #include "isoSurfaceBaseMethods.H"
45 
46 
47 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
48 
49 namespace Foam
50 {
52 }
53 
54 
55 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
56 
57 Foam::scalar Foam::isoSurfaceCell::isoFraction
58 (
59  const scalar s0,
60  const scalar s1
61 ) const
62 {
63  const scalar d = s1-s0;
64 
65  if (mag(d) > VSMALL)
66  {
67  return (iso_-s0)/d;
68  }
69 
70  return -1.0;
71 }
72 
73 
74 Foam::labelPair Foam::isoSurfaceCell::findCommonPoints
75 (
76  const labelledTri& tri0,
77  const labelledTri& tri1
78 )
79 {
80  labelPair common(-1, -1);
81 
82  label fp0 = 0;
83  label fp1 = tri1.find(tri0[fp0]);
84 
85  if (fp1 == -1)
86  {
87  fp0 = 1;
88  fp1 = tri1.find(tri0[fp0]);
89  }
90 
91  if (fp1 != -1)
92  {
93  // So tri0[fp0] is tri1[fp1]
94 
95  // Find next common point
96  label fp0p1 = tri0.fcIndex(fp0);
97  label fp1p1 = tri1.fcIndex(fp1);
98  label fp1m1 = tri1.rcIndex(fp1);
99 
100  if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
101  {
102  common[0] = tri0[fp0];
103  common[1] = tri0[fp0p1];
104  }
105  }
106  return common;
107 }
108 
109 
110 Foam::point Foam::isoSurfaceCell::calcCentre(const triSurface& s)
111 {
112  vector sum = Zero;
113 
114  forAll(s, i)
115  {
116  sum += s[i].centre(s.points());
117  }
118  return sum/s.size();
119 }
120 
121 
122 Foam::pointIndexHit Foam::isoSurfaceCell::collapseSurface
123 (
124  const label celli,
125  pointField& localPoints,
126  DynamicList<labelledTri, 64>& localTris
127 ) const
128 {
129  pointIndexHit info(false, Zero, localTris.size());
130 
131  if (localTris.size() == 1)
132  {
133  const labelledTri& tri = localTris[0];
134  info.setPoint(tri.centre(localPoints));
135  info.setHit();
136  }
137  else if (localTris.size() == 2)
138  {
139  // Check if the two triangles share an edge.
140  const labelledTri& tri0 = localTris[0];
141  const labelledTri& tri1 = localTris[1];
142 
143  labelPair shared = findCommonPoints(tri0, tri1);
144 
145  if (shared[0] != -1)
146  {
147  const vector n0 = tri0.areaNormal(localPoints);
148  const vector n1 = tri1.areaNormal(localPoints);
149 
150  // Merge any zero-sized triangles,
151  // or if they point in the same direction.
152 
153  if
154  (
155  mag(n0) <= ROOTVSMALL
156  || mag(n1) <= ROOTVSMALL
157  || (n0 & n1) >= 0
158  )
159  {
160  info.setPoint
161  (
162  0.5
163  * (
164  tri0.centre(localPoints)
165  + tri1.centre(localPoints)
166  )
167  );
168  info.setHit();
169  }
170  }
171  }
172  else if (localTris.size())
173  {
174  // Check if single region. Rare situation.
175  triSurface surf
176  (
177  localTris,
179  localPoints,
180  true
181  );
182  localTris.clearStorage();
183 
184  labelList faceZone;
185  label nZones = surf.markZones
186  (
187  boolList(surf.nEdges(), false),
188  faceZone
189  );
190 
191  if (nZones == 1)
192  {
193  // Check that all normals make a decent angle
194  scalar minCos = GREAT;
195  const vector& n0 = surf.faceNormals()[0];
196  for (label i = 1; i < surf.size(); ++i)
197  {
198  scalar cosAngle = (n0 & surf.faceNormals()[i]);
199  if (cosAngle < minCos)
200  {
201  minCos = cosAngle;
202  }
203  }
204 
205  if (minCos > 0)
206  {
207  info.setPoint(calcCentre(surf));
208  info.setHit();
209  }
210  }
211  }
212 
213  return info;
214 }
215 
216 
217 void Foam::isoSurfaceCell::calcSnappedCc
218 (
219  const bitSet& isTet,
220  const scalarField& cVals,
221  const scalarField& pVals,
222 
223  DynamicList<point>& snappedPoints,
224  labelList& snappedCc
225 ) const
226 {
227  const pointField& cc = mesh_.cellCentres();
228  const pointField& pts = mesh_.points();
229 
230  snappedCc.setSize(mesh_.nCells());
231  snappedCc = -1;
232 
233  // Work arrays
234  DynamicList<point, 64> localPoints(64);
235  DynamicList<labelledTri, 64> localTris(64);
236  Map<label> pointToLocal(64);
237 
238  forAll(mesh_.cells(), celli)
239  {
240  if (cellCutType_[celli] == cutType::CUT) // No tet cuts
241  {
242  const scalar cVal = cVals[celli];
243 
244  const cell& cFaces = mesh_.cells()[celli];
245 
246  localPoints.clear();
247  localTris.clear();
248  pointToLocal.clear();
249 
250  // Create points for all intersections close to cell centre
251  // (i.e. from pyramid edges)
252 
253  for (const label facei : cFaces)
254  {
255  const face& f = mesh_.faces()[facei];
256 
257  for (const label pointi : f)
258  {
259  scalar s = isoFraction(cVal, pVals[pointi]);
260 
261  if (s >= 0.0 && s <= 0.5)
262  {
263  if (pointToLocal.insert(pointi, localPoints.size()))
264  {
265  localPoints.append((1.0-s)*cc[celli]+s*pts[pointi]);
266  }
267  }
268  }
269  }
270 
271  if (localPoints.size() == 1)
272  {
273  // No need for any analysis.
274  snappedCc[celli] = snappedPoints.size();
275  snappedPoints.append(localPoints[0]);
276 
277  //Pout<< "cell:" << celli
278  // << " at " << mesh_.cellCentres()[celli]
279  // << " collapsing " << localPoints
280  // << " intersections down to "
281  // << snappedPoints[snappedCc[celli]] << endl;
282  }
283  else if (localPoints.size() == 2)
284  {
285  //? No need for any analysis.???
286  snappedCc[celli] = snappedPoints.size();
287  snappedPoints.append(0.5*(localPoints[0]+localPoints[1]));
288 
289  //Pout<< "cell:" << celli
290  // << " at " << mesh_.cellCentres()[celli]
291  // << " collapsing " << localPoints
292  // << " intersections down to "
293  // << snappedPoints[snappedCc[celli]] << endl;
294  }
295  else if (localPoints.size())
296  {
297  // Need to analyse
298  for (const label facei : cFaces)
299  {
300  const face& f = mesh_.faces()[facei];
301 
302  // Do a tetrahedralisation. Each face to cc becomes pyr.
303  // Each pyr gets split into tets by diagonalisation
304  // of face.
305 
306  const label fp0 = mesh_.tetBasePtIs()[facei];
307  label fp = f.fcIndex(fp0);
308  for (label i = 2; i < f.size(); ++i)
309  {
310  label nextFp = f.fcIndex(fp);
311  triFace tri(f[fp0], f[fp], f[nextFp]);
312 
313  // Get fractions for the three edges to cell centre
314  FixedList<scalar, 3> s(3);
315  s[0] = isoFraction(cVal, pVals[tri[0]]);
316  s[1] = isoFraction(cVal, pVals[tri[1]]);
317  s[2] = isoFraction(cVal, pVals[tri[2]]);
318 
319  if
320  (
321  (s[0] >= 0.0 && s[0] <= 0.5)
322  && (s[1] >= 0.0 && s[1] <= 0.5)
323  && (s[2] >= 0.0 && s[2] <= 0.5)
324  )
325  {
326  if
327  (
328  (mesh_.faceOwner()[facei] == celli)
329  == (cVal >= pVals[tri[0]])
330  )
331  {
332  localTris.append
333  (
334  labelledTri
335  (
336  pointToLocal[tri[1]],
337  pointToLocal[tri[0]],
338  pointToLocal[tri[2]],
339  0
340  )
341  );
342  }
343  else
344  {
345  localTris.append
346  (
347  labelledTri
348  (
349  pointToLocal[tri[0]],
350  pointToLocal[tri[1]],
351  pointToLocal[tri[2]],
352  0
353  )
354  );
355  }
356  }
357 
358  fp = nextFp;
359  }
360  }
361 
362  pointField surfPoints;
363  surfPoints.transfer(localPoints);
364  pointIndexHit info = collapseSurface
365  (
366  celli,
367  surfPoints,
368  localTris
369  );
370 
371  if (info.hit())
372  {
373  snappedCc[celli] = snappedPoints.size();
374  snappedPoints.append(info.hitPoint());
375 
376  //Pout<< "cell:" << celli
377  // << " at " << mesh_.cellCentres()[celli]
378  // << " collapsing " << surfPoints
379  // << " intersections down to "
380  // << snappedPoints[snappedCc[celli]] << endl;
381  }
382  }
383  }
384  }
385 }
386 
387 
388 void Foam::isoSurfaceCell::genPointTris
389 (
390  const scalarField& cellValues,
391  const scalarField& pointValues,
392  const label pointi,
393  const label facei,
394  const label celli,
395  DynamicList<point, 64>& localTriPoints
396 ) const
397 {
398  const pointField& cc = mesh_.cellCentres();
399  const pointField& pts = mesh_.points();
400  const face& f = mesh_.faces()[facei];
401 
402  const label fp0 = mesh_.tetBasePtIs()[facei];
403  label fp = f.fcIndex(fp0);
404  for (label i = 2; i < f.size(); ++i)
405  {
406  label nextFp = f.fcIndex(fp);
407  triFace tri(f[fp0], f[fp], f[nextFp]);
408 
409  label index = tri.find(pointi);
410 
411  if (index == -1)
412  {
413  continue;
414  }
415 
416  // Tet between index..index-1, index..index+1, index..cc
417  label b = tri[tri.fcIndex(index)];
418  label c = tri[tri.rcIndex(index)];
419 
420  // Get fractions for the three edges emanating from point
421  FixedList<scalar, 3> s(3);
422  s[0] = isoFraction(pointValues[pointi], pointValues[b]);
423  s[1] = isoFraction(pointValues[pointi], pointValues[c]);
424  s[2] = isoFraction(pointValues[pointi], cellValues[celli]);
425 
426  if
427  (
428  (s[0] >= 0.0 && s[0] <= 0.5)
429  && (s[1] >= 0.0 && s[1] <= 0.5)
430  && (s[2] >= 0.0 && s[2] <= 0.5)
431  )
432  {
433  point p0 = (1.0-s[0])*pts[pointi] + s[0]*pts[b];
434  point p1 = (1.0-s[1])*pts[pointi] + s[1]*pts[c];
435  point p2 = (1.0-s[2])*pts[pointi] + s[2]*cc[celli];
436 
437  if
438  (
439  (mesh_.faceOwner()[facei] == celli)
440  == (pointValues[pointi] > cellValues[celli])
441  )
442  {
443  localTriPoints.append(p0);
444  localTriPoints.append(p1);
445  localTriPoints.append(p2);
446  }
447  else
448  {
449  localTriPoints.append(p1);
450  localTriPoints.append(p0);
451  localTriPoints.append(p2);
452  }
453  }
454 
455  fp = nextFp;
456  }
457 }
458 
459 
460 void Foam::isoSurfaceCell::genPointTris
461 (
462  const scalarField& pointValues,
463  const label pointi,
464  const label facei,
465  const label celli,
466  DynamicList<point, 64>& localTriPoints
467 ) const
468 {
469  const pointField& pts = mesh_.points();
470  const cell& cFaces = mesh_.cells()[celli];
471 
472  // Make tet from this face to the 4th point (same as cellcentre in
473  // non-tet cells)
474  const face& f = mesh_.faces()[facei];
475 
476  // Find 4th point
477  label ccPointi = -1;
478  for (const label cfacei : cFaces)
479  {
480  const face& f1 = mesh_.faces()[cfacei];
481  for (const label p1 : f1)
482  {
483  if (!f.found(p1))
484  {
485  ccPointi = p1;
486  break;
487  }
488  }
489  if (ccPointi != -1)
490  {
491  break;
492  }
493  }
494 
495 
496  // Tet between index..index-1, index..index+1, index..cc
497  label index = f.find(pointi);
498  label b = f[f.fcIndex(index)];
499  label c = f[f.rcIndex(index)];
500 
501  //Pout<< " p0:" << pointi << " b:" << b << " c:" << c
502  //<< " d:" << ccPointi << endl;
503 
504  // Get fractions for the three edges emanating from point
505  FixedList<scalar, 3> s(3);
506  s[0] = isoFraction(pointValues[pointi], pointValues[b]);
507  s[1] = isoFraction(pointValues[pointi], pointValues[c]);
508  s[2] = isoFraction(pointValues[pointi], pointValues[ccPointi]);
509 
510  if
511  (
512  (s[0] >= 0.0 && s[0] <= 0.5)
513  && (s[1] >= 0.0 && s[1] <= 0.5)
514  && (s[2] >= 0.0 && s[2] <= 0.5)
515  )
516  {
517  point p0 = (1.0-s[0])*pts[pointi] + s[0]*pts[b];
518  point p1 = (1.0-s[1])*pts[pointi] + s[1]*pts[c];
519  point p2 = (1.0-s[2])*pts[pointi] + s[2]*pts[ccPointi];
520 
521  if (mesh_.faceOwner()[facei] != celli)
522  {
523  localTriPoints.append(p0);
524  localTriPoints.append(p1);
525  localTriPoints.append(p2);
526  }
527  else
528  {
529  localTriPoints.append(p1);
530  localTriPoints.append(p0);
531  localTriPoints.append(p2);
532  }
533  }
534 }
535 
536 
537 void Foam::isoSurfaceCell::calcSnappedPoint
538 (
539  const bitSet& isTet,
540  const scalarField& cVals,
541  const scalarField& pVals,
542 
543  DynamicList<point>& snappedPoints,
544  labelList& snappedPoint
545 ) const
546 {
547  const labelList& faceOwn = mesh_.faceOwner();
548  const labelList& faceNei = mesh_.faceNeighbour();
549 
550  // Determine if point is on boundary. Points on boundaries are never
551  // snapped. Coupled boundaries are handled explicitly so not marked here.
552  bitSet isBoundaryPoint(mesh_.nPoints());
553  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
554 
555  for (const polyPatch& pp : patches)
556  {
557  if (!pp.coupled())
558  {
559  for (const label facei : pp.range())
560  {
561  const face& f = mesh_.faces()[facei];
562 
563  isBoundaryPoint.set(f);
564  }
565  }
566  }
567 
568 
569  const point greatPoint(GREAT, GREAT, GREAT);
570 
571  pointField collapsedPoint(mesh_.nPoints(), greatPoint);
572 
573 
574  // Work arrays
575  DynamicList<point, 64> localTriPoints(100);
576  labelHashSet localPointCells(100);
577 
578  forAll(mesh_.pointFaces(), pointi)
579  {
580  constexpr uint8_t realCut(cutType::CUT | cutType::TETCUT);
581 
582  if (isBoundaryPoint.test(pointi))
583  {
584  continue;
585  }
586 
587  const labelList& pFaces = mesh_.pointFaces()[pointi];
588 
589  bool anyCut = false;
590 
591  for (const label facei : pFaces)
592  {
593  if
594  (
595  (cellCutType_[faceOwn[facei]] & realCut) != 0
596  || (
597  mesh_.isInternalFace(facei)
598  && (cellCutType_[faceNei[facei]] & realCut) != 0
599  )
600  )
601  {
602  anyCut = true;
603  break;
604  }
605  }
606 
607  if (!anyCut)
608  {
609  continue;
610  }
611 
612 
613  // Do a pointCells walk (by using pointFaces)
614 
615  localPointCells.clear();
616  localTriPoints.clear();
617 
618  for (const label facei : pFaces)
619  {
620  const label own = faceOwn[facei];
621 
622  if (isTet.test(own))
623  {
624  // Since tets have no cell centre to include make sure
625  // we only generate a triangle once per point.
626  if (localPointCells.insert(own))
627  {
628  genPointTris(pVals, pointi, facei, own, localTriPoints);
629  }
630  }
631  else
632  {
633  genPointTris
634  (
635  cVals,
636  pVals,
637  pointi,
638  facei,
639  own,
640  localTriPoints
641  );
642  }
643 
644  if (mesh_.isInternalFace(facei))
645  {
646  const label nei = faceNei[facei];
647 
648  if (isTet.test(nei))
649  {
650  if (localPointCells.insert(nei))
651  {
652  genPointTris(pVals, pointi, facei, nei, localTriPoints);
653  }
654  }
655  else
656  {
657  genPointTris
658  (
659  cVals,
660  pVals,
661  pointi,
662  facei,
663  nei,
664  localTriPoints
665  );
666  }
667  }
668  }
669 
670  if (localTriPoints.size() == 3)
671  {
672  // Single triangle. No need for any analysis. Average points.
674  points.transfer(localTriPoints);
675  collapsedPoint[pointi] = sum(points)/points.size();
676 
677  //Pout<< " point:" << pointi
678  // << " replacing coord:" << mesh_.points()[pointi]
679  // << " by average:" << collapsedPoint[pointi] << endl;
680  }
681  else if (localTriPoints.size())
682  {
683  // Convert points into triSurface.
684 
685  // Merge points and compact out non-valid triangles
686  labelList triMap; // merged to unmerged triangle
687  labelList triPointReverseMap; // unmerged to merged point
688  triSurface surf
689  (
690  stitchTriPoints
691  (
692  false, // do not check for duplicate tris
693  localTriPoints,
694  triPointReverseMap,
695  triMap
696  )
697  );
698 
699  labelList faceZone;
700  label nZones = surf.markZones
701  (
702  boolList(surf.nEdges(), false),
703  faceZone
704  );
705 
706  if (nZones == 1)
707  {
708  // Check that all normals make a decent angle
709  scalar minCos = GREAT;
710  const vector& n0 = surf.faceNormals()[0];
711  for (label i = 1; i < surf.size(); ++i)
712  {
713  const vector& n = surf.faceNormals()[i];
714  scalar cosAngle = (n0 & n);
715  if (cosAngle < minCos)
716  {
717  minCos = cosAngle;
718  }
719  }
720  if (minCos > 0)
721  {
722  collapsedPoint[pointi] = calcCentre(surf);
723  }
724  }
725  }
726  }
727 
729  (
730  mesh_,
731  collapsedPoint,
732  minMagSqrEqOp<point>(),
733  greatPoint
734  );
735 
736  snappedPoint.setSize(mesh_.nPoints());
737  snappedPoint = -1;
738 
739  forAll(collapsedPoint, pointi)
740  {
741  // Cannot do == comparison since might be transformed so have
742  // truncation errors.
743  if (magSqr(collapsedPoint[pointi]) < 0.5*magSqr(greatPoint))
744  {
745  snappedPoint[pointi] = snappedPoints.size();
746  snappedPoints.append(collapsedPoint[pointi]);
747  }
748  }
749 }
750 
751 
752 Foam::triSurface Foam::isoSurfaceCell::stitchTriPoints
753 (
754  const bool checkDuplicates,
755  const List<point>& triPoints,
756  labelList& triPointReverseMap, // unmerged to merged point
757  labelList& triMap // merged to unmerged triangle
758 ) const
759 {
760  label nTris = triPoints.size()/3;
761 
762  if ((triPoints.size() % 3) != 0)
763  {
765  << "Problem: number of points " << triPoints.size()
766  << " not a multiple of 3." << abort(FatalError);
767  }
768 
769  pointField newPoints;
771  (
772  triPoints,
773  mergeDistance_,
774  false,
775  triPointReverseMap,
776  newPoints
777  );
778 
779  // Check that enough merged.
780  if (debug)
781  {
782  Pout<< "isoSurfaceCell : merged from " << triPoints.size()
783  << " points down to " << newPoints.size() << endl;
784 
785  pointField newNewPoints;
786  labelList oldToNew;
787  bool hasMerged = mergePoints
788  (
789  newPoints,
790  mergeDistance_,
791  true,
792  oldToNew,
793  newNewPoints
794  );
795 
796  if (hasMerged)
797  {
799  << "Merged points contain duplicates"
800  << " when merging with distance " << mergeDistance_ << endl
801  << "merged:" << newPoints.size() << " re-merged:"
802  << newNewPoints.size()
803  << abort(FatalError);
804  }
805  }
806 
807 
808  List<labelledTri> tris;
809  {
810  DynamicList<labelledTri> dynTris(nTris);
811  label rawPointi = 0;
812  DynamicList<label> newToOldTri(nTris);
813 
814  for (label oldTriI = 0; oldTriI < nTris; ++oldTriI)
815  {
816  labelledTri tri
817  (
818  triPointReverseMap[rawPointi],
819  triPointReverseMap[rawPointi+1],
820  triPointReverseMap[rawPointi+2],
821  0
822  );
823  if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
824  {
825  newToOldTri.append(oldTriI);
826  dynTris.append(tri);
827  }
828 
829  rawPointi += 3;
830  }
831 
832  triMap.transfer(newToOldTri);
833  tris.transfer(dynTris);
834  }
835 
836 
837  // Use face centres to determine 'flat hole' situation (see RMT paper).
838  // Two unconnected triangles get connected because (some of) the edges
839  // separating them get collapsed. Below only checks for duplicate triangles,
840  // not non-manifold edge connectivity.
841  if (checkDuplicates)
842  {
843  if (debug)
844  {
845  Pout<< "isoSurfaceCell : merged from " << nTris
846  << " down to " << tris.size() << " triangles." << endl;
847  }
848 
849  pointField centres(tris.size());
850  forAll(tris, triI)
851  {
852  centres[triI] = tris[triI].centre(newPoints);
853  }
854 
855  pointField mergedCentres;
856  labelList oldToMerged;
857  bool hasMerged = mergePoints
858  (
859  centres,
860  mergeDistance_,
861  false,
862  oldToMerged,
863  mergedCentres
864  );
865 
866  if (debug)
867  {
868  Pout<< "isoSurfaceCell : detected "
869  << centres.size()-mergedCentres.size()
870  << " duplicate triangles." << endl;
871  }
872 
873  if (hasMerged)
874  {
875  // Filter out duplicates.
876  label newTriI = 0;
877  DynamicList<label> newToOldTri(tris.size());
878  labelList newToMaster(mergedCentres.size(), -1);
879  forAll(tris, triI)
880  {
881  label mergedI = oldToMerged[triI];
882 
883  if (newToMaster[mergedI] == -1)
884  {
885  newToMaster[mergedI] = triI;
886  newToOldTri.append(triMap[triI]);
887  tris[newTriI++] = tris[triI];
888  }
889  }
890 
891  triMap.transfer(newToOldTri);
892  tris.setSize(newTriI);
893  }
894  }
895 
896  return triSurface(tris, geometricSurfacePatchList(), newPoints, true);
897 }
898 
899 
900 void Foam::isoSurfaceCell::calcAddressing
901 (
902  const triSurface& surf,
903  List<FixedList<label, 3>>& faceEdges,
904  labelList& edgeFace0,
905  labelList& edgeFace1,
906  Map<labelList>& edgeFacesRest
907 ) const
908 {
909  const pointField& points = surf.points();
910 
911  pointField edgeCentres(3*surf.size());
912  label edgeI = 0;
913  forAll(surf, triI)
914  {
915  const labelledTri& tri = surf[triI];
916  edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
917  edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
918  edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
919  }
920 
921  pointField mergedCentres;
922  labelList oldToMerged;
923  bool hasMerged = mergePoints
924  (
925  edgeCentres,
926  mergeDistance_,
927  false,
928  oldToMerged,
929  mergedCentres
930  );
931 
932  if (debug)
933  {
934  Pout<< "isoSurfaceCell : detected "
935  << mergedCentres.size()
936  << " edges on " << surf.size() << " triangles." << endl;
937  }
938 
939  if (!hasMerged)
940  {
941  return;
942  }
943 
944 
945  // Determine faceEdges
946  faceEdges.setSize(surf.size());
947  edgeI = 0;
948  forAll(surf, triI)
949  {
950  faceEdges[triI][0] = oldToMerged[edgeI++];
951  faceEdges[triI][1] = oldToMerged[edgeI++];
952  faceEdges[triI][2] = oldToMerged[edgeI++];
953  }
954 
955 
956  // Determine edgeFaces
957  edgeFace0.setSize(mergedCentres.size());
958  edgeFace0 = -1;
959  edgeFace1.setSize(mergedCentres.size());
960  edgeFace1 = -1;
961  edgeFacesRest.clear();
962 
963  forAll(oldToMerged, oldEdgeI)
964  {
965  label triI = oldEdgeI / 3;
966  label edgeI = oldToMerged[oldEdgeI];
967 
968  if (edgeFace0[edgeI] == -1)
969  {
970  edgeFace0[edgeI] = triI;
971  }
972  else if (edgeFace1[edgeI] == -1)
973  {
974  edgeFace1[edgeI] = triI;
975  }
976  else
977  {
978  //WarningInFunction
979  // << "Edge " << edgeI << " with centre " << mergedCentres[edgeI]
980  // << " used by more than two triangles: " << edgeFace0[edgeI]
981  // << ", "
982  // << edgeFace1[edgeI] << " and " << triI << endl;
983  Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);
984 
985  if (iter != edgeFacesRest.end())
986  {
987  labelList& eFaces = iter();
988  label sz = eFaces.size();
989  eFaces.setSize(sz+1);
990  eFaces[sz] = triI;
991  }
992  else
993  {
994  edgeFacesRest.insert(edgeI, labelList(1, triI));
995  }
996  }
997  }
998 }
999 
1000 
1001 bool Foam::isoSurfaceCell::danglingTriangle
1002 (
1003  const FixedList<label, 3>& fEdges,
1004  const labelList& edgeFace1
1005 )
1006 {
1007  label nOpen = 0;
1008  for (const label edgei : fEdges)
1009  {
1010  if (edgeFace1[edgei] == -1)
1011  {
1012  ++nOpen;
1013  }
1014  }
1015 
1016  return (nOpen == 1 || nOpen == 2 || nOpen == 3);
1017 }
1018 
1019 
1020 Foam::label Foam::isoSurfaceCell::markDanglingTriangles
1021 (
1022  const List<FixedList<label, 3>>& faceEdges,
1023  const labelList& edgeFace0,
1024  const labelList& edgeFace1,
1025  const Map<labelList>& edgeFacesRest,
1026  boolList& keepTriangles
1027 )
1028 {
1029  keepTriangles.setSize(faceEdges.size());
1030  keepTriangles = true;
1031 
1032  label nDangling = 0;
1033 
1034  // Remove any dangling triangles
1035  forAllConstIters(edgeFacesRest, iter)
1036  {
1037  // These are all the non-manifold edges. Filter out all triangles
1038  // with only one connected edge (= this edge)
1039 
1040  const label edgeI = iter.key();
1041  const labelList& otherEdgeFaces = iter.val();
1042 
1043  // Remove all dangling triangles
1044  if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
1045  {
1046  keepTriangles[edgeFace0[edgeI]] = false;
1047  ++nDangling;
1048  }
1049  if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
1050  {
1051  keepTriangles[edgeFace1[edgeI]] = false;
1052  ++nDangling;
1053  }
1054  for (const label triI : otherEdgeFaces)
1055  {
1056  if (danglingTriangle(faceEdges[triI], edgeFace1))
1057  {
1058  keepTriangles[triI] = false;
1059  ++nDangling;
1060  }
1061  }
1062  }
1063  return nDangling;
1064 }
1065 
1066 
1067 Foam::triSurface Foam::isoSurfaceCell::subsetMesh
1068 (
1069  const triSurface& s,
1070  const labelList& newToOldFaces,
1071  labelList& oldToNewPoints,
1072  labelList& newToOldPoints
1073 )
1074 {
1075  const boolList include
1076  (
1077  ListOps::createWithValue<bool>(s.size(), newToOldFaces, true, false)
1078  );
1079 
1080  newToOldPoints.setSize(s.points().size());
1081  oldToNewPoints.setSize(s.points().size());
1082  oldToNewPoints = -1;
1083  {
1084  label pointi = 0;
1085 
1086  forAll(include, oldFacei)
1087  {
1088  if (include[oldFacei])
1089  {
1090  // Renumber labels for face
1091  for (const label oldPointi : s[oldFacei])
1092  {
1093  if (oldToNewPoints[oldPointi] == -1)
1094  {
1095  oldToNewPoints[oldPointi] = pointi;
1096  newToOldPoints[pointi++] = oldPointi;
1097  }
1098  }
1099  }
1100  }
1101  newToOldPoints.setSize(pointi);
1102  }
1103 
1104  // Extract points
1105  pointField newPoints(newToOldPoints.size());
1106  forAll(newToOldPoints, i)
1107  {
1108  newPoints[i] = s.points()[newToOldPoints[i]];
1109  }
1110  // Extract faces
1111  List<labelledTri> newTriangles(newToOldFaces.size());
1112 
1113  forAll(newToOldFaces, i)
1114  {
1115  // Get old vertex labels
1116  const labelledTri& tri = s[newToOldFaces[i]];
1117 
1118  newTriangles[i][0] = oldToNewPoints[tri[0]];
1119  newTriangles[i][1] = oldToNewPoints[tri[1]];
1120  newTriangles[i][2] = oldToNewPoints[tri[2]];
1121  newTriangles[i].region() = tri.region();
1122  }
1123 
1124  // Reuse storage.
1125  return triSurface(newTriangles, s.patches(), newPoints, true);
1126 }
1127 
1128 
1129 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1130 
1133  const polyMesh& mesh,
1134  const scalarField& cellValues,
1135  const scalarField& pointValues,
1136  const scalar iso,
1137  const isoSurfaceParams& params,
1138  const bitSet& ignoreCells
1139 )
1140 :
1141  isoSurfaceBase(mesh, cellValues, pointValues, iso, params),
1142  mergeDistance_(params.mergeTol()*mesh.bounds().mag()),
1143  cellCutType_(mesh.nCells(), cutType::UNVISITED)
1144 {
1145  const bool regularise = (params.filter() != filterType::NONE);
1146 
1147  if (debug)
1148  {
1149  Pout<< "isoSurfaceCell:" << nl
1150  << " cell min/max : " << minMax(cVals_) << nl
1151  << " point min/max : " << minMax(pVals_) << nl
1152  << " isoValue : " << iso << nl
1153  << " filter : " << Switch(regularise) << nl
1154  << " mergeTol : " << params.mergeTol() << nl
1155  << " mesh span : " << mesh.bounds().mag() << nl
1156  << " mergeDistance : " << mergeDistance_ << nl
1157  << " ignoreCells : " << ignoreCells.count()
1158  << " / " << cVals_.size() << nl
1159  << endl;
1160  }
1161 
1162 
1163  label nBlockedCells = 0;
1164 
1165  // Mark ignoreCells as blocked
1166  nBlockedCells += blockCells(cellCutType_, ignoreCells);
1167 
1168 
1169  // Some processor domains may require tetBasePtIs and others do not.
1170  // Do now to ensure things stay synchronized.
1171  (void)mesh_.tetBasePtIs();
1172 
1173 
1174  // Calculate a tet/non-tet filter
1175  bitSet isTet(mesh_.nCells());
1176  {
1177  for (label celli = 0; celli < mesh_.nCells(); ++celli)
1178  {
1179  if (tetMatcher::test(mesh_, celli))
1180  {
1181  isTet.set(celli);
1182  }
1183  }
1184  }
1185 
1186  // Determine cell cuts
1187  nCutCells_ = calcCellCuts(cellCutType_);
1188 
1189  if (debug)
1190  {
1191  Pout<< "isoSurfaceCell : candidate cells cut "
1192  << nCutCells_
1193  << " blocked " << nBlockedCells
1194  << " total " << mesh_.nCells() << endl;
1195  }
1196 
1197  if (debug && isA<fvMesh>(mesh))
1198  {
1199  const auto& fvmesh = dynamicCast<const fvMesh>(mesh);
1200 
1201  volScalarField debugField
1202  (
1203  IOobject
1204  (
1205  "isoSurfaceCell.cutType",
1206  fvmesh.time().timeName(),
1207  fvmesh.time(),
1210  false
1211  ),
1212  fvmesh,
1214  );
1215 
1216  auto& debugFld = debugField.primitiveFieldRef();
1217 
1218  forAll(cellCutType_, celli)
1219  {
1220  debugFld[celli] = cellCutType_[celli];
1221  }
1222 
1223  Pout<< "Writing cut types:"
1224  << debugField.objectPath() << endl;
1225 
1226  debugField.write();
1227  }
1228 
1229 
1230  DynamicList<point> snappedPoints(nCutCells_);
1231 
1232  // Per cc -1 or a point inside snappedPoints.
1233  labelList snappedCc;
1234  if (regularise)
1235  {
1236  calcSnappedCc
1237  (
1238  isTet,
1239  cellValues,
1240  pointValues,
1241  snappedPoints,
1242  snappedCc
1243  );
1244  }
1245  else
1246  {
1247  snappedCc.setSize(mesh_.nCells());
1248  snappedCc = -1;
1249  }
1250 
1251  if (debug)
1252  {
1253  Pout<< "isoSurfaceCell : shifted " << snappedPoints.size()
1254  << " cell centres to intersection." << endl;
1255  }
1256 
1257  snappedPoints.shrink();
1258  label nCellSnaps = snappedPoints.size();
1259 
1260  // Per point -1 or a point inside snappedPoints.
1261  labelList snappedPoint;
1262  if (regularise)
1263  {
1264  calcSnappedPoint
1265  (
1266  isTet,
1267  cellValues,
1268  pointValues,
1269  snappedPoints,
1270  snappedPoint
1271  );
1272  }
1273  else
1274  {
1275  snappedPoint.setSize(mesh_.nPoints());
1276  snappedPoint = -1;
1277  }
1278 
1279  if (debug)
1280  {
1281  Pout<< "isoSurfaceCell : shifted " << snappedPoints.size()-nCellSnaps
1282  << " vertices to intersection." << endl;
1283  }
1284 
1285 
1286  // Use a triSurface as a temporary for various operations
1287  triSurface tmpsurf;
1288 
1289  {
1290  DynamicList<point> triPoints(nCutCells_);
1291  DynamicList<label> triMeshCells(nCutCells_);
1292 
1293  generateTriPoints
1294  (
1295  cellValues,
1296  pointValues,
1297 
1298  mesh_.cellCentres(),
1299  mesh_.points(),
1300 
1301  snappedPoints,
1302  snappedCc,
1303  snappedPoint,
1304 
1305  triPoints,
1306  triMeshCells
1307  );
1308 
1309  if (debug)
1310  {
1311  Pout<< "isoSurfaceCell : generated " << triMeshCells.size()
1312  << " unmerged triangles." << endl;
1313  }
1314 
1315 
1316  label nOldPoints = triPoints.size();
1317 
1318  // Trimmed to original triangle
1319  DynamicList<label> trimTriMap;
1320  // Trimmed to original point
1321  labelList trimTriPointMap;
1322  if (getClipBounds().valid())
1323  {
1324  isoSurfacePoint::trimToBox
1325  (
1326  treeBoundBox(getClipBounds()),
1327  triPoints, // new points
1328  trimTriMap, // map from (new) triangle to original
1329  trimTriPointMap, // map from (new) point to original
1330  interpolatedPoints_, // labels of newly introduced points
1331  interpolatedOldPoints_, // and their interpolation
1332  interpolationWeights_
1333  );
1334  triMeshCells = labelField(triMeshCells, trimTriMap);
1335  }
1336 
1337 
1338  // Merge points and compact out non-valid triangles
1339  labelList triMap; // merged to unmerged triangle
1340  tmpsurf = stitchTriPoints
1341  (
1342  regularise, // check for duplicate tris
1343  triPoints,
1344  triPointMergeMap_, // unmerged to merged point
1345  triMap // merged to unmerged triangle
1346  );
1347 
1348  if (debug)
1349  {
1350  Pout<< "isoSurfaceCell : generated " << triMap.size()
1351  << " merged triangles." << endl;
1352  }
1353 
1354  if (getClipBounds().valid())
1355  {
1356  // Adjust interpolatedPoints_
1357  inplaceRenumber(triPointMergeMap_, interpolatedPoints_);
1358 
1359  // Adjust triPointMergeMap_
1360  labelList newTriPointMergeMap(nOldPoints, -1);
1361  forAll(trimTriPointMap, trimPointI)
1362  {
1363  label oldPointI = trimTriPointMap[trimPointI];
1364  if (oldPointI >= 0)
1365  {
1366  label pointI = triPointMergeMap_[trimPointI];
1367  if (pointI >= 0)
1368  {
1369  newTriPointMergeMap[oldPointI] = pointI;
1370  }
1371  }
1372  }
1373  triPointMergeMap_.transfer(newTriPointMergeMap);
1374  }
1375 
1376  meshCells_.setSize(triMap.size());
1377  forAll(triMap, i)
1378  {
1379  meshCells_[i] = triMeshCells[triMap[i]];
1380  }
1381  }
1382 
1383 
1384  if (debug)
1385  {
1386  Pout<< "isoSurfaceCell : checking " << tmpsurf.size()
1387  << " triangles for validity." << endl;
1388 
1389  forAll(tmpsurf, triI)
1390  {
1391  triSurfaceTools::validTri(tmpsurf, triI);
1392  }
1393  }
1394 
1395 
1396  if (regularise)
1397  {
1398  List<FixedList<label, 3>> faceEdges;
1399  labelList edgeFace0, edgeFace1;
1400  Map<labelList> edgeFacesRest;
1401 
1402 
1403  while (true)
1404  {
1405  // Calculate addressing
1406  calcAddressing
1407  (
1408  tmpsurf,
1409  faceEdges,
1410  edgeFace0,
1411  edgeFace1,
1412  edgeFacesRest
1413  );
1414 
1415  // See if any dangling triangles
1416  boolList keepTriangles;
1417  label nDangling = markDanglingTriangles
1418  (
1419  faceEdges,
1420  edgeFace0,
1421  edgeFace1,
1422  edgeFacesRest,
1423  keepTriangles
1424  );
1425 
1426  if (debug)
1427  {
1428  Pout<< "isoSurfaceCell : detected " << nDangling
1429  << " dangling triangles." << endl;
1430  }
1431 
1432  if (nDangling == 0)
1433  {
1434  break;
1435  }
1436 
1437  // Create face map (new to old)
1438  labelList subsetTriMap(findIndices(keepTriangles, true));
1439 
1440  labelList subsetPointMap;
1441  labelList reversePointMap;
1442  tmpsurf = subsetMesh
1443  (
1444  tmpsurf,
1445  subsetTriMap,
1446  reversePointMap,
1447  subsetPointMap
1448  );
1449  meshCells_ = labelField(meshCells_, subsetTriMap);
1450  inplaceRenumber(reversePointMap, triPointMergeMap_);
1451  }
1452  }
1453 
1454 
1455  // Transfer to mesh storage. Note, an iso-surface has no zones
1456  {
1457  // Recover the pointField
1458  pointField pts;
1459  tmpsurf.swapPoints(pts);
1460 
1461  // Transcribe from triFace to face
1462  faceList faces;
1463  tmpsurf.triFaceFaces(faces);
1464 
1465  tmpsurf.clearOut();
1466 
1467  Mesh updated(std::move(pts), std::move(faces), surfZoneList());
1468 
1469  this->Mesh::transfer(updated);
1470  }
1471 }
1472 
1473 
1474 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
nZones
label nZones
Definition: interpolatedFaces.H:24
Foam::IOobject::NO_WRITE
Definition: IOobject.H:195
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:67
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:169
Foam::Switch
A simple wrapper around bool so that it can be read as a word: true/false, on/off,...
Definition: Switch.H:77
Foam::boundBox::mag
scalar mag() const
The magnitude of the bounding box span.
Definition: boundBoxI.H:133
Foam::triSurface::clearOut
void clearOut()
Definition: triSurface.C:566
Foam::isoSurfaceBase
Low-level components common to various iso-surface algorithms.
Definition: isoSurfaceBase.H:66
s
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
Definition: gmvOutputSpray.H:25
Foam::bitSet
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:63
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::DynamicList< point >
isoSurfaceCell.H
Foam::triSurface::triFaceFaces
void triFaceFaces(List< face > &plainFaceList) const
Create a list of faces from the triFaces.
Definition: triSurface.C:723
Foam::treeBoundBox
Standard boundBox with extra functionality for use in octree.
Definition: treeBoundBox.H:86
Foam::bitSet::set
void set(const bitSet &bitset)
Set specified bits from another bitset.
Definition: bitSetI.H:574
Foam::Map
A HashTable to objects of type <T> with a label key.
Definition: lumpedPointController.H:69
Foam::FixedList::size
static constexpr label size() noexcept
Return the number of elements in the FixedList.
Definition: FixedList.H:416
Foam::boolList
List< bool > boolList
A List of bools.
Definition: List.H:65
triPoints.H
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::syncTools::syncPointPositions
static void syncPointPositions(const polyMesh &mesh, List< point > &positions, const CombineOp &cop, const point &nullValue)
Synchronize locations on all mesh points.
Definition: syncTools.H:201
defineIsoSurfaceInterpolateMethods
#define defineIsoSurfaceInterpolateMethods(ThisClass)
Definition: isoSurfaceBaseMethods.H:54
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::Map::iterator
typename parent_type::iterator iterator
Definition: Map.H:69
triSurface.H
polyMesh.H
Foam::DynamicList::shrink
DynamicList< T, SizeMin > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:434
syncTools.H
Foam::bitSet::count
unsigned int count(const bool on=true) const
Count number of bits set.
Definition: bitSetI.H:499
Foam::inplaceRenumber
void inplaceRenumber(const labelUList &oldToNew, IntListType &input)
Inplace renumber the values (not the indices) of a list.
Definition: ListOpsTemplates.C:61
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:77
Foam::isoSurfaceCell
A surface formed by the iso value. After "Polygonising A Scalar Field Using Tetrahedrons",...
Definition: isoSurfaceCell.H:66
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
isoSurfaceBaseMethods.H
Convenience macros for instantiating iso-surface interpolate methods.
Foam::geometricSurfacePatchList
List< geometricSurfacePatch > geometricSurfacePatchList
A List of geometricSurfacePatch.
Definition: geometricSurfacePatchList.H:47
Foam::labelPair
Pair< label > labelPair
A pair of labels.
Definition: Pair.H:54
Foam::isoSurfaceParams::filter
filterType filter() const noexcept
Get current filter type.
Definition: isoSurfaceParams.H:225
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::triPoints
Triangle storage. Null constructable (unfortunately triangle<point, point> is not)
Definition: triPoints.H:52
Foam::primitiveMesh::nCells
label nCells() const noexcept
Number of mesh cells.
Definition: primitiveMeshI.H:96
Foam::PointIndexHit
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:52
Foam::constant::physicoChemical::b
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
Foam::Field< scalar >
isoSurfacePoint.H
Foam::labelField
Field< label > labelField
Specialisation of Field<T> for label.
Definition: labelField.H:52
Foam::triSurface
Triangulated surface description with patch information.
Definition: triSurface.H:76
Foam::triSurface::swapPoints
virtual void swapPoints(pointField &pts)
Swap points. Similar to movePoints, but returns the old points.
Definition: triSurface.C:625
Foam::List::setSize
void setSize(const label n)
Alias for resize()
Definition: List.H:222
Foam::List::transfer
void transfer(List< T > &list)
Definition: List.C:456
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition: dimensionedScalarFwd.H:42
Foam::isoSurfaceBase::iso_
const scalar iso_
Iso value.
Definition: isoSurfaceBase.H:111
Foam::FatalError
error FatalError
Foam::isoSurfaceParams
Preferences for controlling iso-surface algorithms.
Definition: isoSurfaceParams.H:107
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::Ostream::write
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::mergePoints
label mergePoints(const PointList &points, const scalar mergeTol, const bool verbose, labelList &pointMap, typename PointList::const_reference origin=PointList::value_type::zero)
Sorts and merges points. All points closer than/equal mergeTol get merged.
Foam::GeometricField::primitiveFieldRef
Internal::FieldType & primitiveFieldRef(const bool updateAccessTime=true)
Return a reference to the internal field.
Definition: GeometricField.C:766
Foam::polyMesh::bounds
const boundBox & bounds() const
Return mesh bounding box.
Definition: polyMesh.H:450
Time.H
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Foam::pointIndexHit
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::Pair
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: Pair.H:54
forAllConstIters
forAllConstIters(mixture.phases(), phase)
Definition: pEqn.H:28
f
labelList f(nPoints)
Foam::Vector< scalar >
Foam::List< label >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
points
const pointField & points
Definition: gmvOutputHeader.H:1
dictionary.H
Foam::sum
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Definition: DimensionedFieldFunctions.C:327
tetMatcher.H
patches
const polyBoundaryMesh & patches
Definition: convertProcessorPatches.H:65
Foam::surfZoneList
List< surfZone > surfZoneList
Definition: surfZoneList.H:47
Foam::List::set
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
Definition: List.H:341
mergePoints.H
Merge points. See below.
Foam::tetMatcher::test
static bool test(const UList< face > &faces)
Definition: tetMatcher.C:89
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
Foam::findIndices
labelList findIndices(const ListType &input, typename ListType::const_reference val, label start=0)
Linear search to find all occurrences of given element.
Foam::minMax
MinMax< label > minMax(const labelHashSet &set)
Find the min/max values of labelHashSet.
Definition: hashSets.C:61
Foam::isoSurfaceCell::isoSurfaceCell
isoSurfaceCell(const polyMesh &mesh, const scalarField &cellValues, const scalarField &pointValues, const scalar iso, const isoSurfaceParams &params=isoSurfaceParams(), const bitSet &ignoreCells=bitSet())
Construct from cell and point values.
Definition: isoSurfaceCell.C:1132
p0
const volScalarField & p0
Definition: EEqn.H:36
Foam::point
vector point
Point is a vector.
Definition: point.H:43
Foam::triSurfaceTools::validTri
static bool validTri(const triSurface &, const label facei, const bool verbose=true)
Check single triangle for (topological) validity.
Definition: triSurfaceTools.C:2696
Foam::labelHashSet
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
Definition: HashSet.H:85
Foam::isoSurfaceParams::mergeTol
scalar mergeTol() const noexcept
Get current merge tolerance.
Definition: isoSurfaceParams.H:249
Foam::GeometricField< scalar, fvPatchField, volMesh >
Foam::MeshedSurface< face >
Foam::IOobject::NO_READ
Definition: IOobject.H:188
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(combustionModel, 0)
triFace
face triFace(3)
triSurfaceTools.H
pFaces
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
Foam::dimless
const dimensionSet dimless
Dimensionless.
Definition: dimensionSets.C:189