isoSurfacePoint.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-2017 OpenFOAM Foundation
9  Copyright (C) 2015-2021 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 "isoSurfacePoint.H"
30 #include "mergePoints.H"
31 #include "slicedVolFields.H"
32 #include "volFields.H"
33 #include "triSurfaceTools.H"
34 #include "triSurface.H"
35 #include "triPoints.H"
36 
37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 
39 #include "isoSurfaceBaseMethods.H"
41 
42 
43 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
44 
45 namespace Foam
46 {
48 }
49 
50 
51 // * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
52 
53 namespace Foam
54 {
55 
56 // Helper class for slicing triangles
57 struct storeOp
58 {
60 
62  :
63  list(tris)
64  {}
65 
66  void operator()(const triPoints& tri)
67  {
68  list.append(tri);
69  }
70 
71  void operator()(triPoints&& tri)
72  {
73  list.append(std::move(tri));
74  }
75 };
76 
77 
78 // Is patch co-located (i.e. non-separated) coupled patch?
79 static inline bool collocatedPatch(const polyPatch& pp)
80 {
81  const auto* cpp = isA<coupledPolyPatch>(pp);
82  return cpp && cpp->parallel() && !cpp->separated();
83 }
84 
85 
86 // Collocated patch, with extra checks
87 #if 0
88 static bool isCollocatedPatch(const coupledPolyPatch& pp)
89 {
90  if (isA<processorPolyPatch>(pp) || isA<cyclicPolyPatch>(pp))
91  {
92  return collocatedPatch(pp);
93  }
94 
96  << "Unhandled coupledPolyPatch type " << pp.type() << nl
97  << abort(FatalError);
98 
99  return false;
100 }
101 #endif
102 
103 } // End namespace Foam
104 
105 
106 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
107 
108 bool Foam::isoSurfacePoint::noTransform(const tensor& tt) const
109 {
110  return
111  (mag(tt.xx()-1) < mergeDistance_)
112  && (mag(tt.yy()-1) < mergeDistance_)
113  && (mag(tt.zz()-1) < mergeDistance_)
114  && (mag(tt.xy()) < mergeDistance_)
115  && (mag(tt.xz()) < mergeDistance_)
116  && (mag(tt.yx()) < mergeDistance_)
117  && (mag(tt.yz()) < mergeDistance_)
118  && (mag(tt.zx()) < mergeDistance_)
119  && (mag(tt.zy()) < mergeDistance_);
120 }
121 
122 
123 Foam::bitSet Foam::isoSurfacePoint::collocatedFaces
124 (
125  const coupledPolyPatch& pp
126 )
127 {
128  // Initialise to false
129  bitSet collocated(pp.size());
130 
131  if (isA<processorPolyPatch>(pp) || isA<cyclicPolyPatch>(pp))
132  {
133  if (collocatedPatch(pp))
134  {
135  // All on
136  collocated = true;
137  }
138  }
139  else
140  {
142  << "Unhandled coupledPolyPatch type " << pp.type()
143  << abort(FatalError);
144  }
145  return collocated;
146 }
147 
148 
149 void Foam::isoSurfacePoint::syncUnseparatedPoints
150 (
151  pointField& pointValues,
152  const point& nullValue
153 ) const
154 {
155  // Until syncPointList handles separated coupled patches with multiple
156  // transforms do our own synchronisation of non-separated patches only
157  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
158 
159  if (Pstream::parRun())
160  {
161  const labelList& procPatches = mesh_.globalData().processorPatches();
162 
163  // Send
164  for (const label patchi : procPatches)
165  {
166  const polyPatch& pp = patches[patchi];
167  const auto& procPatch = refCast<const processorPolyPatch>(pp);
168 
169  if (pp.nPoints() && collocatedPatch(pp))
170  {
171  const labelList& meshPts = procPatch.meshPoints();
172  const labelList& nbrPts = procPatch.neighbPoints();
173 
174  pointField patchInfo(meshPts.size());
175 
176  forAll(nbrPts, pointi)
177  {
178  const label nbrPointi = nbrPts[pointi];
179  patchInfo[nbrPointi] = pointValues[meshPts[pointi]];
180  }
181 
182  OPstream toNbr
183  (
185  procPatch.neighbProcNo()
186  );
187  toNbr << patchInfo;
188  }
189  }
190 
191  // Receive and combine.
192  for (const label patchi : procPatches)
193  {
194  const polyPatch& pp = patches[patchi];
195  const auto& procPatch = refCast<const processorPolyPatch>(pp);
196 
197  if (pp.nPoints() && collocatedPatch(pp))
198  {
199  pointField nbrPatchInfo(procPatch.nPoints());
200  {
201  // We do not know the number of points on the other side
202  // so cannot use Pstream::read.
203  IPstream fromNbr
204  (
206  procPatch.neighbProcNo()
207  );
208  fromNbr >> nbrPatchInfo;
209  }
210 
211  const labelList& meshPts = procPatch.meshPoints();
212 
213  forAll(meshPts, pointi)
214  {
215  const label meshPointi = meshPts[pointi];
216  minEqOp<point>()
217  (
218  pointValues[meshPointi],
219  nbrPatchInfo[pointi]
220  );
221  }
222  }
223  }
224  }
225 
226  // Do the cyclics.
227  for (const polyPatch& pp : patches)
228  {
229  const cyclicPolyPatch* cpp = isA<cyclicPolyPatch>(pp);
230 
231  if (cpp && cpp->owner() && collocatedPatch(*cpp))
232  {
233  // Owner does all.
234 
235  const auto& cycPatch = *cpp;
236  const auto& nbrPatch = cycPatch.neighbPatch();
237 
238  const edgeList& coupledPoints = cycPatch.coupledPoints();
239  const labelList& meshPts = cycPatch.meshPoints();
240  const labelList& nbrMeshPoints = nbrPatch.meshPoints();
241 
242  pointField half0Values(coupledPoints.size());
243  pointField half1Values(coupledPoints.size());
244 
245  forAll(coupledPoints, i)
246  {
247  const edge& e = coupledPoints[i];
248  half0Values[i] = pointValues[meshPts[e[0]]];
249  half1Values[i] = pointValues[nbrMeshPoints[e[1]]];
250  }
251 
252  forAll(coupledPoints, i)
253  {
254  const edge& e = coupledPoints[i];
255  const label p0 = meshPts[e[0]];
256  const label p1 = nbrMeshPoints[e[1]];
257 
258  minEqOp<point>()(pointValues[p0], half1Values[i]);
259  minEqOp<point>()(pointValues[p1], half0Values[i]);
260  }
261  }
262  }
263 
264  // Synchronize multiple shared points.
265  const globalMeshData& pd = mesh_.globalData();
266 
267  if (pd.nGlobalPoints() > 0)
268  {
269  // Values on shared points.
270  pointField sharedPts(pd.nGlobalPoints(), nullValue);
271 
272  forAll(pd.sharedPointLabels(), i)
273  {
274  const label meshPointi = pd.sharedPointLabels()[i];
275  // Fill my entries in the shared points
276  sharedPts[pd.sharedPointAddr()[i]] = pointValues[meshPointi];
277  }
278 
279  // Combine on master.
280  Pstream::listCombineGather(sharedPts, minEqOp<point>());
281  Pstream::listCombineScatter(sharedPts);
282 
283  // Now we will all have the same information. Merge it back with
284  // my local information.
285  forAll(pd.sharedPointLabels(), i)
286  {
287  const label meshPointi = pd.sharedPointLabels()[i];
288  pointValues[meshPointi] = sharedPts[pd.sharedPointAddr()[i]];
289  }
290  }
291 }
292 
293 
294 Foam::scalar Foam::isoSurfacePoint::isoFraction
295 (
296  const scalar s0,
297  const scalar s1
298 ) const
299 {
300  const scalar d = s1-s0;
301 
302  if (mag(d) > VSMALL)
303  {
304  return (iso_-s0)/d;
305  }
306 
307  return -1.0;
308 }
309 
310 
311 bool Foam::isoSurfacePoint::isEdgeOfFaceCut
312 (
313  const scalarField& pVals,
314  const face& f,
315  const bool ownLower,
316  const bool neiLower
317 ) const
318 {
319  forAll(f, fp)
320  {
321  const bool fpLower = (pVals[f[fp]] < iso_);
322 
323  if
324  (
325  fpLower != ownLower
326  || fpLower != neiLower
327  || fpLower != (pVals[f[f.fcIndex(fp)]] < iso_)
328  )
329  {
330  return true;
331  }
332  }
333 
334  return false;
335 }
336 
337 
338 void Foam::isoSurfacePoint::getNeighbour
339 (
340  const labelList& boundaryRegion,
341  const volVectorField& meshC,
342  const volScalarField& cVals,
343  const label celli,
344  const label facei,
345  scalar& nbrValue,
346  point& nbrPoint
347 ) const
348 {
349  const labelList& own = mesh_.faceOwner();
350  const labelList& nei = mesh_.faceNeighbour();
351 
352  if (mesh_.isInternalFace(facei))
353  {
354  const label nbr = (own[facei] == celli ? nei[facei] : own[facei]);
355  nbrValue = cVals[nbr];
356  nbrPoint = meshC[nbr];
357  }
358  else
359  {
360  const label bFacei = facei-mesh_.nInternalFaces();
361  const label patchi = boundaryRegion[bFacei];
362  const polyPatch& pp = mesh_.boundaryMesh()[patchi];
363  const label patchFacei = facei-pp.start();
364 
365  nbrValue = cVals.boundaryField()[patchi][patchFacei];
366  nbrPoint = meshC.boundaryField()[patchi][patchFacei];
367  }
368 }
369 
370 
371 void Foam::isoSurfacePoint::calcCutTypes
372 (
373  const labelList& boundaryRegion,
374  const volVectorField& meshC,
375  const volScalarField& cVals,
376  const scalarField& pVals
377 )
378 {
379  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
380  const labelList& own = mesh_.faceOwner();
381  const labelList& nei = mesh_.faceNeighbour();
382 
383  faceCutType_.resize(mesh_.nFaces());
384  faceCutType_ = cutType::NOTCUT;
385 
386  for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
387  {
388  const scalar ownValue = cVals[own[facei]];
389  const bool ownLower = (ownValue < iso_);
390 
391  scalar nbrValue;
392  point nbrPoint;
393  getNeighbour
394  (
395  boundaryRegion,
396  meshC,
397  cVals,
398  own[facei],
399  facei,
400  nbrValue,
401  nbrPoint
402  );
403 
404  const bool neiLower = (nbrValue < iso_);
405 
406  if (ownLower != neiLower)
407  {
408  faceCutType_[facei] = cutType::CUT;
409  }
410  else
411  {
412  // Is mesh edge cut?
413  // - by looping over all the edges of the face.
414  const face& f = mesh_.faces()[facei];
415 
416  if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
417  {
418  faceCutType_[facei] = cutType::CUT;
419  }
420  }
421  }
422 
423  for (const polyPatch& pp : patches)
424  {
425  for (const label facei : pp.range())
426  {
427  const scalar ownValue = cVals[own[facei]];
428  const bool ownLower = (ownValue < iso_);
429 
430  scalar nbrValue;
431  point nbrPoint;
432  getNeighbour
433  (
434  boundaryRegion,
435  meshC,
436  cVals,
437  own[facei],
438  facei,
439  nbrValue,
440  nbrPoint
441  );
442 
443  const bool neiLower = (nbrValue < iso_);
444 
445  if (ownLower != neiLower)
446  {
447  faceCutType_[facei] = cutType::CUT;
448  }
449  else
450  {
451  // Is mesh edge cut?
452  // - by looping over all the edges of the face.
453  const face& f = mesh_.faces()[facei];
454 
455  if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
456  {
457  faceCutType_[facei] = cutType::CUT;
458  }
459  }
460  }
461  }
462 
463  nCutCells_ = 0;
464  cellCutType_.resize(mesh_.nCells());
465  cellCutType_ = cutType::NOTCUT;
466 
467 
468  // Propagate internal face cuts into the cells.
469 
470  for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
471  {
472  if (faceCutType_[facei] == cutType::NOTCUT)
473  {
474  continue;
475  }
476 
477  if (cellCutType_[own[facei]] == cutType::NOTCUT)
478  {
479  cellCutType_[own[facei]] = cutType::CUT;
480  ++nCutCells_;
481  }
482  if (cellCutType_[nei[facei]] == cutType::NOTCUT)
483  {
484  cellCutType_[nei[facei]] = cutType::CUT;
485  ++nCutCells_;
486  }
487  }
488 
489 
490  // Propagate boundary face cuts into the cells.
491 
492  for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); ++facei)
493  {
494  if (faceCutType_[facei] == cutType::NOTCUT)
495  {
496  continue;
497  }
498 
499  if (cellCutType_[own[facei]] == cutType::NOTCUT)
500  {
501  cellCutType_[own[facei]] = cutType::CUT;
502  ++nCutCells_;
503  }
504  }
505 
506  if (debug)
507  {
508  Pout<< "isoSurfacePoint : candidate cut cells "
509  << nCutCells_ << " / " << mesh_.nCells() << endl;
510  }
511 }
512 
513 
514 Foam::point Foam::isoSurfacePoint::calcCentre(const triSurface& s)
515 {
516  // Calculate centre of surface.
517 
518  vector sum = Zero;
519 
520  forAll(s, i)
521  {
522  sum += s[i].centre(s.points());
523  }
524  return sum/s.size();
525 }
526 
527 
528 void Foam::isoSurfacePoint::calcSnappedCc
529 (
530  const labelList& boundaryRegion,
531  const volVectorField& meshC,
532  const volScalarField& cVals,
533  const scalarField& pVals,
534 
535  DynamicList<point>& snappedPoints,
536  labelList& snappedCc
537 ) const
538 {
539  const pointField& pts = mesh_.points();
540  const pointField& cc = mesh_.cellCentres();
541 
542  snappedCc.setSize(mesh_.nCells());
543  snappedCc = -1;
544 
545  // Work arrays
546  DynamicList<point, 64> localTriPoints(64);
547 
548  forAll(mesh_.cells(), celli)
549  {
550  if (cellCutType_[celli] == cutType::CUT)
551  {
552  const scalar cVal = cVals[celli];
553 
554  localTriPoints.clear();
555  label nOther = 0;
556  point otherPointSum = Zero;
557 
558  // Create points for all intersections close to cell centre
559  // (i.e. from pyramid edges)
560 
561  for (const label facei : mesh_.cells()[celli])
562  {
563  const face& f = mesh_.faces()[facei];
564 
565  scalar nbrValue;
566  point nbrPoint;
567  getNeighbour
568  (
569  boundaryRegion,
570  meshC,
571  cVals,
572  celli,
573  facei,
574  nbrValue,
575  nbrPoint
576  );
577 
578  // Calculate intersection points of edges to cell centre
579  FixedList<scalar, 3> s;
580  FixedList<point, 3> pt;
581 
582  // From cc to neighbour cc.
583  s[2] = isoFraction(cVal, nbrValue);
584  pt[2] = (1.0-s[2])*cc[celli] + s[2]*nbrPoint;
585 
586  forAll(f, fp)
587  {
588  // From cc to fp
589  label p0 = f[fp];
590  s[0] = isoFraction(cVal, pVals[p0]);
591  pt[0] = (1.0-s[0])*cc[celli] + s[0]*pts[p0];
592 
593  // From cc to fp+1
594  label p1 = f[f.fcIndex(fp)];
595  s[1] = isoFraction(cVal, pVals[p1]);
596  pt[1] = (1.0-s[1])*cc[celli] + s[1]*pts[p1];
597 
598  if
599  (
600  (s[0] >= 0.0 && s[0] <= 0.5)
601  && (s[1] >= 0.0 && s[1] <= 0.5)
602  && (s[2] >= 0.0 && s[2] <= 0.5)
603  )
604  {
605  localTriPoints.append(pt[0]);
606  localTriPoints.append(pt[1]);
607  localTriPoints.append(pt[2]);
608  }
609  else
610  {
611  // Get average of all other points
612  forAll(s, i)
613  {
614  if (s[i] >= 0.0 && s[i] <= 0.5)
615  {
616  otherPointSum += pt[i];
617  ++nOther;
618  }
619  }
620  }
621  }
622  }
623 
624  if (localTriPoints.size() == 0)
625  {
626  // No complete triangles. Get average of all intersection
627  // points.
628  if (nOther > 0)
629  {
630  snappedCc[celli] = snappedPoints.size();
631  snappedPoints.append(otherPointSum/nOther);
632 
633  //Pout<< " point:" << pointi
634  // << " replacing coord:" << mesh_.points()[pointi]
635  // << " by average:" << collapsedPoint[pointi] << endl;
636  }
637  }
638  else if (localTriPoints.size() == 3)
639  {
640  // Single triangle. No need for any analysis. Average points.
641  snappedCc[celli] = snappedPoints.size();
642  snappedPoints.append(sum(localTriPoints)/3);
643  localTriPoints.clear();
644 
645  //Pout<< " point:" << pointi
646  // << " replacing coord:" << mesh_.points()[pointi]
647  // << " by average:" << collapsedPoint[pointi] << endl;
648  }
649  else
650  {
651  // Convert points into triSurface.
652 
653  // Merge points and compact out non-valid triangles
654  labelList triMap; // merged to unmerged triangle
655  labelList triPointReverseMap; // unmerged to merged point
656  triSurface surf
657  (
658  stitchTriPoints
659  (
660  false, // do not check for duplicate tris
661  localTriPoints,
662  triPointReverseMap,
663  triMap
664  )
665  );
666 
667  labelList faceZone;
668  label nZones = surf.markZones
669  (
670  boolList(surf.nEdges(), false),
671  faceZone
672  );
673 
674  if (nZones == 1)
675  {
676  snappedCc[celli] = snappedPoints.size();
677  snappedPoints.append(calcCentre(surf));
678  //Pout<< " point:" << pointi << " nZones:" << nZones
679  // << " replacing coord:" << mesh_.points()[pointi]
680  // << " by average:" << collapsedPoint[pointi] << endl;
681  }
682  }
683  }
684  }
685 }
686 
687 
688 void Foam::isoSurfacePoint::calcSnappedPoint
689 (
690  const bitSet& isBoundaryPoint,
691  const labelList& boundaryRegion,
692  const volVectorField& meshC,
693  const volScalarField& cVals,
694  const scalarField& pVals,
695 
696  DynamicList<point>& snappedPoints,
697  labelList& snappedPoint
698 ) const
699 {
700  const pointField& pts = mesh_.points();
701  const pointField& cc = mesh_.cellCentres();
702 
703  pointField collapsedPoint(mesh_.nPoints(), point::max);
704 
705  // Work arrays
706  DynamicList<point, 64> localTriPoints(100);
707 
708  forAll(mesh_.pointFaces(), pointi)
709  {
710  if (isBoundaryPoint.test(pointi))
711  {
712  continue;
713  }
714 
715  const labelList& pFaces = mesh_.pointFaces()[pointi];
716 
717  bool anyCut = false;
718 
719  for (const label facei : pFaces)
720  {
721  if (faceCutType_[facei] == cutType::CUT)
722  {
723  anyCut = true;
724  break;
725  }
726  }
727 
728  if (!anyCut)
729  {
730  continue;
731  }
732 
733 
734  localTriPoints.clear();
735  label nOther = 0;
736  point otherPointSum = Zero;
737 
738  for (const label facei : pFaces)
739  {
740  // Create points for all intersections close to point
741  // (i.e. from pyramid edges)
742  const face& f = mesh_.faces()[facei];
743  const label own = mesh_.faceOwner()[facei];
744 
745  // Get neighbour value
746  scalar nbrValue;
747  point nbrPoint;
748  getNeighbour
749  (
750  boundaryRegion,
751  meshC,
752  cVals,
753  own,
754  facei,
755  nbrValue,
756  nbrPoint
757  );
758 
759  // Calculate intersection points of edges emanating from point
760  FixedList<scalar, 4> s;
761  FixedList<point, 4> pt;
762 
763  label fp = f.find(pointi);
764  s[0] = isoFraction(pVals[pointi], cVals[own]);
765  pt[0] = (1.0-s[0])*pts[pointi] + s[0]*cc[own];
766 
767  s[1] = isoFraction(pVals[pointi], nbrValue);
768  pt[1] = (1.0-s[1])*pts[pointi] + s[1]*nbrPoint;
769 
770  label nextPointi = f[f.fcIndex(fp)];
771  s[2] = isoFraction(pVals[pointi], pVals[nextPointi]);
772  pt[2] = (1.0-s[2])*pts[pointi] + s[2]*pts[nextPointi];
773 
774  label prevPointi = f[f.rcIndex(fp)];
775  s[3] = isoFraction(pVals[pointi], pVals[prevPointi]);
776  pt[3] = (1.0-s[3])*pts[pointi] + s[3]*pts[prevPointi];
777 
778  if
779  (
780  (s[0] >= 0.0 && s[0] <= 0.5)
781  && (s[1] >= 0.0 && s[1] <= 0.5)
782  && (s[2] >= 0.0 && s[2] <= 0.5)
783  )
784  {
785  localTriPoints.append(pt[0]);
786  localTriPoints.append(pt[1]);
787  localTriPoints.append(pt[2]);
788  }
789  if
790  (
791  (s[0] >= 0.0 && s[0] <= 0.5)
792  && (s[1] >= 0.0 && s[1] <= 0.5)
793  && (s[3] >= 0.0 && s[3] <= 0.5)
794  )
795  {
796  localTriPoints.append(pt[3]);
797  localTriPoints.append(pt[0]);
798  localTriPoints.append(pt[1]);
799  }
800 
801  // Get average of points as fallback
802  forAll(s, i)
803  {
804  if (s[i] >= 0.0 && s[i] <= 0.5)
805  {
806  otherPointSum += pt[i];
807  ++nOther;
808  }
809  }
810  }
811 
812  if (localTriPoints.size() == 0)
813  {
814  // No complete triangles. Get average of all intersection
815  // points.
816  if (nOther > 0)
817  {
818  collapsedPoint[pointi] = otherPointSum/nOther;
819  }
820  }
821  else if (localTriPoints.size() == 3)
822  {
823  // Single triangle. No need for any analysis. Average points.
825  points.transfer(localTriPoints);
826  collapsedPoint[pointi] = sum(points)/points.size();
827  }
828  else
829  {
830  // Convert points into triSurface.
831 
832  // Merge points and compact out non-valid triangles
833  labelList triMap; // merged to unmerged triangle
834  labelList triPointReverseMap; // unmerged to merged point
835  triSurface surf
836  (
837  stitchTriPoints
838  (
839  false, // do not check for duplicate tris
840  localTriPoints,
841  triPointReverseMap,
842  triMap
843  )
844  );
845 
846  labelList faceZone;
847  label nZones = surf.markZones
848  (
849  boolList(surf.nEdges(), false),
850  faceZone
851  );
852 
853  if (nZones == 1)
854  {
855  collapsedPoint[pointi] = calcCentre(surf);
856  }
857  }
858  }
859 
860 
861  // Synchronise snap location
862  syncUnseparatedPoints(collapsedPoint, point::max);
863 
864 
865  snappedPoint.setSize(mesh_.nPoints());
866  snappedPoint = -1;
867 
868  forAll(collapsedPoint, pointi)
869  {
870  if (collapsedPoint[pointi] != point::max)
871  {
872  snappedPoint[pointi] = snappedPoints.size();
873  snappedPoints.append(collapsedPoint[pointi]);
874  }
875  }
876 }
877 
878 
879 Foam::triSurface Foam::isoSurfacePoint::stitchTriPoints
880 (
881  const bool checkDuplicates,
882  const List<point>& triPoints,
883  labelList& triPointReverseMap, // unmerged to merged point
884  labelList& triMap // merged to unmerged triangle
885 ) const
886 {
887  label nTris = triPoints.size()/3;
888 
889  if ((triPoints.size() % 3) != 0)
890  {
892  << "Problem: number of points " << triPoints.size()
893  << " not a multiple of 3." << abort(FatalError);
894  }
895 
896  pointField newPoints;
898  (
899  triPoints,
900  mergeDistance_,
901  false,
902  triPointReverseMap,
903  newPoints
904  );
905 
906  // Check that enough merged.
907  if (debug)
908  {
909  Pout<< "isoSurfacePoint : merged from " << triPoints.size()
910  << " down to " << newPoints.size() << " unique points." << endl;
911 
912  pointField newNewPoints;
913  labelList oldToNew;
914  bool hasMerged = mergePoints
915  (
916  newPoints,
917  mergeDistance_,
918  true,
919  oldToNew,
920  newNewPoints
921  );
922 
923  if (hasMerged)
924  {
926  << "Merged points contain duplicates"
927  << " when merging with distance " << mergeDistance_ << endl
928  << "merged:" << newPoints.size() << " re-merged:"
929  << newNewPoints.size()
930  << abort(FatalError);
931  }
932  }
933 
934 
935  List<labelledTri> tris;
936  {
937  DynamicList<labelledTri> dynTris(nTris);
938  label rawPointi = 0;
939  DynamicList<label> newToOldTri(nTris);
940 
941  for (label oldTriI = 0; oldTriI < nTris; oldTriI++)
942  {
943  labelledTri tri
944  (
945  triPointReverseMap[rawPointi],
946  triPointReverseMap[rawPointi+1],
947  triPointReverseMap[rawPointi+2],
948  0
949  );
950  rawPointi += 3;
951 
952  if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
953  {
954  newToOldTri.append(oldTriI);
955  dynTris.append(tri);
956  }
957  }
958 
959  triMap.transfer(newToOldTri);
960  tris.transfer(dynTris);
961  }
962 
963 
964 
965  // Determine 'flat hole' situation (see RMT paper).
966  // Two unconnected triangles get connected because (some of) the edges
967  // separating them get collapsed. Below only checks for duplicate triangles,
968  // not non-manifold edge connectivity.
969  if (checkDuplicates)
970  {
971  labelListList pointFaces;
972  invertManyToMany(newPoints.size(), tris, pointFaces);
973 
974  // Filter out duplicates.
975  DynamicList<label> newToOldTri(tris.size());
976 
977  forAll(tris, triI)
978  {
979  const labelledTri& tri = tris[triI];
980  const labelList& pFaces = pointFaces[tri[0]];
981 
982  // Find the maximum of any duplicates. Maximum since the tris
983  // below triI
984  // get overwritten so we cannot use them in a comparison.
985  label dupTriI = -1;
986  forAll(pFaces, i)
987  {
988  label nbrTriI = pFaces[i];
989 
990  if (nbrTriI > triI && (tris[nbrTriI] == tri))
991  {
992  //Pout<< "Duplicate : " << triI << " verts:" << tri
993  // << " to " << nbrTriI << " verts:" << tris[nbrTriI]
994  // << endl;
995  dupTriI = nbrTriI;
996  break;
997  }
998  }
999 
1000  if (dupTriI == -1)
1001  {
1002  // There is no (higher numbered) duplicate triangle
1003  label newTriI = newToOldTri.size();
1004  newToOldTri.append(triMap[triI]);
1005  tris[newTriI] = tris[triI];
1006  }
1007  }
1008 
1009  triMap.transfer(newToOldTri);
1010  tris.setSize(triMap.size());
1011 
1012  if (debug)
1013  {
1014  Pout<< "isoSurfacePoint : merged from " << nTris
1015  << " down to " << tris.size() << " unique triangles." << endl;
1016  }
1017 
1018 
1019  if (debug)
1020  {
1021  triSurface surf(tris, geometricSurfacePatchList(), newPoints);
1022 
1023  forAll(surf, facei)
1024  {
1025  const labelledTri& f = surf[facei];
1026  const labelList& fFaces = surf.faceFaces()[facei];
1027 
1028  forAll(fFaces, i)
1029  {
1030  label nbrFacei = fFaces[i];
1031 
1032  if (nbrFacei <= facei)
1033  {
1034  // lower numbered faces already checked
1035  continue;
1036  }
1037 
1038  const labelledTri& nbrF = surf[nbrFacei];
1039 
1040  if (f == nbrF)
1041  {
1043  << "Check : "
1044  << " triangle " << facei << " vertices " << f
1045  << " fc:" << f.centre(surf.points())
1046  << " has the same vertices as triangle " << nbrFacei
1047  << " vertices " << nbrF
1048  << " fc:" << nbrF.centre(surf.points())
1049  << abort(FatalError);
1050  }
1051  }
1052  }
1053  }
1054  }
1055 
1056  return triSurface(tris, geometricSurfacePatchList(), newPoints, true);
1057 }
1058 
1059 
1060 void Foam::isoSurfacePoint::trimToPlanes
1061 (
1062  const PtrList<plane>& planes,
1063  const triPointRef& tri,
1064  DynamicList<point>& newTriPoints
1065 )
1066 {
1067  // Buffer for generated triangles
1068  DynamicList<triPoints> insideTrisA;
1069  storeOp insideOpA(insideTrisA);
1070 
1071  // Buffer for generated triangles
1072  DynamicList<triPoints> insideTrisB;
1073  storeOp insideOpB(insideTrisB);
1074 
1075  triPointRef::dummyOp dop;
1076 
1077  // Store starting triangle in insideTrisA
1078  insideOpA(triPoints(tri.a(), tri.b(), tri.c()));
1079 
1080 
1081  bool useA = true;
1082 
1083  forAll(planes, faceI)
1084  {
1085  const plane& pln = planes[faceI];
1086 
1087  if (useA)
1088  {
1089  insideTrisB.clear();
1090  for (const triPoints& tri : insideTrisA)
1091  {
1092  triPointRef(tri).sliceWithPlane(pln, insideOpB, dop);
1093  }
1094  }
1095  else
1096  {
1097  insideTrisA.clear();
1098  for (const triPoints& tri : insideTrisB)
1099  {
1100  triPointRef(tri).sliceWithPlane(pln, insideOpA, dop);
1101  }
1102  }
1103  useA = !useA;
1104  }
1105 
1106 
1107  DynamicList<triPoints>& newTris = (useA ? insideTrisA : insideTrisB);
1108 
1109  newTriPoints.reserve(newTriPoints.size() + 3*newTris.size());
1110 
1111  // Transfer
1112  for (const triPoints& tri : newTris)
1113  {
1114  newTriPoints.append(tri[0]);
1115  newTriPoints.append(tri[1]);
1116  newTriPoints.append(tri[2]);
1117  }
1118 }
1119 
1120 
1121 void Foam::isoSurfacePoint::trimToBox
1122 (
1123  const treeBoundBox& bb,
1124  DynamicList<point>& triPoints,
1125  DynamicList<label>& triMap // trimmed to original tris
1126 )
1127 {
1128  if (debug)
1129  {
1130  Pout<< "isoSurfacePoint : trimming to " << bb << endl;
1131  }
1132 
1133  // Generate inwards pointing planes
1134  PtrList<plane> planes(treeBoundBox::faceNormals.size());
1136  {
1137  const vector& n = treeBoundBox::faceNormals[faceI];
1138  planes.set(faceI, new plane(bb.faceCentre(faceI), -n));
1139  }
1140 
1141  label nTris = triPoints.size()/3;
1142 
1143  DynamicList<point> newTriPoints(triPoints.size()/16);
1144  triMap.setCapacity(nTris/16);
1145 
1146  label vertI = 0;
1147  for (label triI = 0; triI < nTris; triI++)
1148  {
1149  const point& p0 = triPoints[vertI++];
1150  const point& p1 = triPoints[vertI++];
1151  const point& p2 = triPoints[vertI++];
1152 
1153  label oldNPoints = newTriPoints.size();
1154  trimToPlanes
1155  (
1156  planes,
1157  triPointRef(p0, p1, p2),
1158  newTriPoints
1159  );
1160 
1161  label nCells = (newTriPoints.size()-oldNPoints)/3;
1162  for (label i = 0; i < nCells; i++)
1163  {
1164  triMap.append(triI);
1165  }
1166  }
1167 
1168  if (debug)
1169  {
1170  Pout<< "isoSurfacePoint : trimmed from " << nTris
1171  << " down to " << triMap.size()
1172  << " triangles." << endl;
1173  }
1174 
1175  triPoints.transfer(newTriPoints);
1176 }
1177 
1178 
1179 void Foam::isoSurfacePoint::trimToBox
1180 (
1181  const treeBoundBox& bb,
1182  DynamicList<point>& triPoints, // new points
1183  DynamicList<label>& triMap, // map from (new) triangle to original
1184  labelList& triPointMap, // map from (new) point to original
1185  labelList& interpolatedPoints, // labels of newly introduced points
1186  List<FixedList<label, 3>>& interpolatedOldPoints,// and their interpolation
1187  List<FixedList<scalar, 3>>& interpolationWeights
1188 )
1189 {
1190  const List<point> oldTriPoints(triPoints);
1191 
1192  // Trim triPoints, return map
1193  trimToBox(bb, triPoints, triMap);
1194 
1195 
1196  // Find point correspondence:
1197  // - one-to-one for preserved original points (triPointMap)
1198  // - interpolation for newly introduced points
1199  // (interpolatedOldPoints)
1200  label sz = oldTriPoints.size()/100;
1201  DynamicList<label> dynInterpolatedPoints(sz);
1202  DynamicList<FixedList<label, 3>> dynInterpolatedOldPoints(sz);
1203  DynamicList<FixedList<scalar, 3>> dynInterpolationWeights(sz);
1204 
1205 
1206  triPointMap.setSize(triPoints.size());
1207  forAll(triMap, triI)
1208  {
1209  label oldTriI = triMap[triI];
1210 
1211  // Find point correspondence. Assumes coordinate is bit-copy.
1212  for (label i = 0; i < 3; i++)
1213  {
1214  label pointI = 3*triI+i;
1215  const point& pt = triPoints[pointI];
1216 
1217  // Compare to old-triangle's points
1218  label matchPointI = -1;
1219  for (label j = 0; j < 3; j++)
1220  {
1221  label oldPointI = 3*oldTriI+j;
1222  if (pt == oldTriPoints[oldPointI])
1223  {
1224  matchPointI = oldPointI;
1225  break;
1226  }
1227  }
1228 
1229  triPointMap[pointI] = matchPointI;
1230 
1231  // If new point: calculate and store interpolation
1232  if (matchPointI == -1)
1233  {
1234  dynInterpolatedPoints.append(pointI);
1235 
1236  FixedList<label, 3> oldPoints
1237  (
1238  {3*oldTriI, 3*oldTriI+1, 3*oldTriI+2}
1239  );
1240  dynInterpolatedOldPoints.append(oldPoints);
1241 
1242  triPointRef tri(oldTriPoints, oldPoints);
1243  barycentric2D bary = tri.pointToBarycentric(pt);
1244  FixedList<scalar, 3> weights({bary.a(), bary.b(), bary.c()});
1245 
1246  dynInterpolationWeights.append(weights);
1247  }
1248  }
1249  }
1250 
1251  interpolatedPoints.transfer(dynInterpolatedPoints);
1252  interpolatedOldPoints.transfer(dynInterpolatedOldPoints);
1253  interpolationWeights.transfer(dynInterpolationWeights);
1254 }
1255 
1256 
1257 Foam::triSurface Foam::isoSurfacePoint::subsetMesh
1258 (
1259  const triSurface& s,
1260  const labelList& newToOldFaces,
1261  labelList& oldToNewPoints,
1262  labelList& newToOldPoints
1263 )
1264 {
1265  const boolList include
1266  (
1267  ListOps::createWithValue<bool>(s.size(), newToOldFaces, true, false)
1268  );
1269 
1270  newToOldPoints.setSize(s.points().size());
1271  oldToNewPoints.setSize(s.points().size());
1272  oldToNewPoints = -1;
1273  {
1274  label pointi = 0;
1275 
1276  forAll(include, oldFacei)
1277  {
1278  if (include[oldFacei])
1279  {
1280  // Renumber labels for triangle
1281  const labelledTri& tri = s[oldFacei];
1282 
1283  forAll(tri, fp)
1284  {
1285  label oldPointi = tri[fp];
1286 
1287  if (oldToNewPoints[oldPointi] == -1)
1288  {
1289  oldToNewPoints[oldPointi] = pointi;
1290  newToOldPoints[pointi++] = oldPointi;
1291  }
1292  }
1293  }
1294  }
1295  newToOldPoints.setSize(pointi);
1296  }
1297 
1298  // Extract points
1299  pointField newPoints(newToOldPoints.size());
1300  forAll(newToOldPoints, i)
1301  {
1302  newPoints[i] = s.points()[newToOldPoints[i]];
1303  }
1304  // Extract faces
1305  List<labelledTri> newTriangles(newToOldFaces.size());
1306 
1307  forAll(newToOldFaces, i)
1308  {
1309  // Get old vertex labels
1310  const labelledTri& tri = s[newToOldFaces[i]];
1311 
1312  newTriangles[i][0] = oldToNewPoints[tri[0]];
1313  newTriangles[i][1] = oldToNewPoints[tri[1]];
1314  newTriangles[i][2] = oldToNewPoints[tri[2]];
1315  newTriangles[i].region() = tri.region();
1316  }
1317 
1318  // Reuse storage.
1319  return triSurface(newTriangles, s.patches(), newPoints, true);
1320 }
1321 
1322 
1323 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1324 
1327  const volScalarField& cellValues,
1328  const scalarField& pointValues,
1329  const scalar iso,
1330  const isoSurfaceParams& params,
1331  const bitSet& /*unused*/
1332 )
1333 :
1335  (
1336  cellValues.mesh(),
1337  cellValues.primitiveField(),
1338  pointValues,
1339  iso,
1340  params
1341  ),
1342  cValsPtr_(nullptr),
1343  mergeDistance_(params.mergeTol()*mesh_.bounds().mag())
1344 {
1345  const bool regularise = (params.filter() != filterType::NONE);
1346  const fvMesh& fvmesh = cellValues.mesh();
1347 
1348  if (debug)
1349  {
1350  Pout<< "isoSurfacePoint:" << nl
1351  << " isoField : " << cellValues.name() << nl
1352  << " cell min/max : " << minMax(cVals_) << nl
1353  << " point min/max : " << minMax(pVals_) << nl
1354  << " isoValue : " << iso << nl
1355  << " filter : " << Switch(regularise) << nl
1356  << " mergeTol : " << params.mergeTol() << nl
1357  << endl;
1358  }
1359 
1360  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1361 
1362  // Rewrite input field
1363  // ~~~~~~~~~~~~~~~~~~~
1364  // Rewrite input volScalarField to have interpolated values
1365  // on separated patches.
1366 
1367  cValsPtr_.reset(adaptPatchFields(cellValues).ptr());
1368 
1369 
1370  // Construct cell centres field consistent with cVals
1371  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1372  // Generate field to interpolate. This is identical to the mesh.C()
1373  // except on separated coupled patches and on empty patches.
1374 
1375  slicedVolVectorField meshC
1376  (
1377  IOobject
1378  (
1379  "C",
1380  fvmesh.pointsInstance(),
1381  fvmesh.meshSubDir,
1382  fvmesh,
1385  false
1386  ),
1387  fvmesh,
1388  dimLength,
1389  fvmesh.cellCentres(),
1390  fvmesh.faceCentres()
1391  );
1392  forAll(patches, patchi)
1393  {
1394  const polyPatch& pp = patches[patchi];
1395 
1396  // Adapt separated coupled (proc and cyclic) patches
1397  if (pp.coupled())
1398  {
1399  auto& pfld = const_cast<fvPatchVectorField&>
1400  (
1401  meshC.boundaryField()[patchi]
1402  );
1403 
1404  bitSet isCollocated
1405  (
1406  collocatedFaces(refCast<const coupledPolyPatch>(pp))
1407  );
1408 
1409  forAll(isCollocated, i)
1410  {
1411  if (!isCollocated[i])
1412  {
1413  pfld[i] = mesh_.faceCentres()[pp.start()+i];
1414  }
1415  }
1416  }
1417  else if (isA<emptyPolyPatch>(pp))
1418  {
1419  auto& bfld = const_cast<slicedVolVectorField::Boundary&>
1420  (
1421  meshC.boundaryField()
1422  );
1423 
1424  // Clear old value. Cannot resize it since is a slice.
1425  bfld.set(patchi, nullptr);
1426 
1427  // Set new value we can change
1428  bfld.set
1429  (
1430  patchi,
1432  (
1433  fvmesh.boundary()[patchi],
1434  meshC
1435  )
1436  );
1437 
1438  // Change to face centres
1439  bfld[patchi] = pp.patchSlice(mesh_.faceCentres());
1440  }
1441  }
1442 
1443 
1444  // Pre-calculate patch-per-face to avoid whichPatch call.
1445  labelList boundaryRegion(mesh_.nBoundaryFaces());
1446 
1447  for (const polyPatch& pp : patches)
1448  {
1449  SubList<label>(boundaryRegion, pp.size(), pp.offset()) = pp.index();
1450  }
1451 
1452 
1453  // Determine if any cut through face/cell
1454  calcCutTypes(boundaryRegion, meshC, cValsPtr_(), pVals_);
1455 
1456  if (debug)
1457  {
1458  volScalarField debugField
1459  (
1460  IOobject
1461  (
1462  "isoSurfacePoint.cutType",
1463  fvmesh.time().timeName(),
1464  fvmesh.time(),
1467  false
1468  ),
1469  fvmesh,
1471  );
1472 
1473  auto& debugFld = debugField.primitiveFieldRef();
1474 
1475  forAll(cellCutType_, celli)
1476  {
1477  debugFld[celli] = cellCutType_[celli];
1478  }
1479 
1480  Pout<< "Writing cut types:"
1481  << debugField.objectPath() << endl;
1482 
1483  debugField.write();
1484  }
1485 
1486 
1487  DynamicList<point> snappedPoints(nCutCells_);
1488 
1489  // Per cc -1 or a point inside snappedPoints.
1490  labelList snappedCc;
1491  if (regularise)
1492  {
1493  calcSnappedCc
1494  (
1496  meshC,
1497  cValsPtr_(),
1498  pVals_,
1499 
1500  snappedPoints,
1501  snappedCc
1502  );
1503  }
1504  else
1505  {
1506  snappedCc.setSize(mesh_.nCells());
1507  snappedCc = -1;
1508  }
1509 
1510 
1511 
1512  if (debug)
1513  {
1514  Pout<< "isoSurfacePoint : shifted " << snappedPoints.size()
1515  << " cell centres to intersection." << endl;
1516  }
1517 
1518  label nCellSnaps = snappedPoints.size();
1519 
1520 
1521  // Per point -1 or a point inside snappedPoints.
1522  labelList snappedPoint;
1523  if (regularise)
1524  {
1525  // Determine if point is on boundary.
1526  bitSet isBoundaryPoint(mesh_.nPoints());
1527 
1528  for (const polyPatch& pp : patches)
1529  {
1530  // Mark all boundary points that are not physically coupled
1531  // (so anything but collocated coupled patches)
1532 
1533  if (pp.coupled())
1534  {
1535  const coupledPolyPatch& cpp =
1536  refCast<const coupledPolyPatch>(pp);
1537 
1538  bitSet isCollocated(collocatedFaces(cpp));
1539 
1540  forAll(isCollocated, i)
1541  {
1542  if (!isCollocated[i])
1543  {
1544  const face& f = mesh_.faces()[cpp.start()+i];
1545 
1546  isBoundaryPoint.set(f);
1547  }
1548  }
1549  }
1550  else
1551  {
1552  forAll(pp, i)
1553  {
1554  const face& f = mesh_.faces()[pp.start()+i];
1555 
1556  isBoundaryPoint.set(f);
1557  }
1558  }
1559  }
1560 
1561  calcSnappedPoint
1562  (
1563  isBoundaryPoint,
1565  meshC,
1566  cValsPtr_(),
1567  pVals_,
1568 
1569  snappedPoints,
1570  snappedPoint
1571  );
1572  }
1573  else
1574  {
1575  snappedPoint.setSize(mesh_.nPoints());
1576  snappedPoint = -1;
1577  }
1578 
1579  if (debug)
1580  {
1581  Pout<< "isoSurfacePoint : shifted " << snappedPoints.size()-nCellSnaps
1582  << " vertices to intersection." << endl;
1583  }
1584 
1585 
1586  // Use a triSurface as a temporary for various operations
1587  triSurface tmpsurf;
1588 
1589  {
1590  DynamicList<point> triPoints(3*nCutCells_);
1591  DynamicList<label> triMeshCells(nCutCells_);
1592 
1593  generateTriPoints
1594  (
1595  cValsPtr_(),
1596  pVals_,
1597 
1598  meshC,
1599  mesh_.points(),
1600 
1601  snappedPoints,
1602  snappedCc,
1603  snappedPoint,
1604 
1605  triPoints, // 3 points of the triangle
1606  triMeshCells // per triangle the originating cell
1607  );
1608 
1609  if (debug)
1610  {
1611  Pout<< "isoSurfacePoint : generated " << triMeshCells.size()
1612  << " unmerged triangles from " << triPoints.size()
1613  << " unmerged points." << endl;
1614  }
1615 
1616  label nOldPoints = triPoints.size();
1617 
1618  // Trimmed to original triangle
1619  DynamicList<label> trimTriMap;
1620  // Trimmed to original point
1621  labelList trimTriPointMap;
1622  if (getClipBounds().valid())
1623  {
1624  trimToBox
1625  (
1626  treeBoundBox(getClipBounds()),
1627  triPoints, // new points
1628  trimTriMap, // map from (new) triangle to original
1629  trimTriPointMap, // map from (new) point to original
1630  interpolatedPoints_, // labels of newly introduced points
1631  interpolatedOldPoints_, // and their interpolation
1632  interpolationWeights_
1633  );
1634  triMeshCells = labelField(triMeshCells, trimTriMap);
1635  }
1636 
1637 
1638  // Merge points and compact out non-valid triangles
1639  labelList triMap; // merged to unmerged triangle
1640  tmpsurf = stitchTriPoints
1641  (
1642  true, // check for duplicate tris
1643  triPoints,
1644  triPointMergeMap_, // unmerged to merged point
1645  triMap
1646  );
1647 
1648  if (debug)
1649  {
1650  Pout<< "isoSurfacePoint : generated " << triMap.size()
1651  << " merged triangles." << endl;
1652  }
1653 
1654 
1655  if (getClipBounds().valid())
1656  {
1657  // Adjust interpolatedPoints_
1658  inplaceRenumber(triPointMergeMap_, interpolatedPoints_);
1659 
1660  // Adjust triPointMergeMap_
1661  labelList newTriPointMergeMap(nOldPoints, -1);
1662  forAll(trimTriPointMap, trimPointI)
1663  {
1664  label oldPointI = trimTriPointMap[trimPointI];
1665  if (oldPointI >= 0)
1666  {
1667  label pointI = triPointMergeMap_[trimPointI];
1668  if (pointI >= 0)
1669  {
1670  newTriPointMergeMap[oldPointI] = pointI;
1671  }
1672  }
1673  }
1674  triPointMergeMap_.transfer(newTriPointMergeMap);
1675  }
1676 
1677  meshCells_.setSize(triMap.size());
1678  forAll(triMap, i)
1679  {
1680  meshCells_[i] = triMeshCells[triMap[i]];
1681  }
1682  }
1683 
1684  if (debug)
1685  {
1686  Pout<< "isoSurfacePoint : checking " << tmpsurf.size()
1687  << " triangles for validity." << endl;
1688 
1689  forAll(tmpsurf, facei)
1690  {
1691  triSurfaceTools::validTri(tmpsurf, facei);
1692  }
1693 
1694  fileName stlFile = mesh_.time().path() + ".stl";
1695  Pout<< "Dumping surface to " << stlFile << endl;
1696  tmpsurf.write(stlFile);
1697  }
1698 
1699 
1700  // Transfer to mesh storage. Note, an iso-surface has no zones
1701  {
1702  // Recover the pointField
1703  pointField pts;
1704  tmpsurf.swapPoints(pts);
1705 
1706  // Transcribe from triFace to face
1707  faceList faces;
1708  tmpsurf.triFaceFaces(faces);
1709 
1710  tmpsurf.clearOut();
1711 
1712  Mesh updated(std::move(pts), std::move(faces), surfZoneList());
1713 
1714  this->Mesh::transfer(updated);
1715  }
1716 }
1717 
1718 
1719 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::fvPatchField< vector >
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
volFields.H
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::storeOp
Definition: isoSurfacePoint.C:57
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::isoSurfacePoint::isoSurfacePoint
isoSurfacePoint(const volScalarField &cellValues, const scalarField &pointValues, const scalar iso, const isoSurfaceParams &params=isoSurfaceParams(), const bitSet &ignoreCells=bitSet())
Construct from cell values and point values.
Definition: isoSurfacePoint.C:1326
Foam::triSurface::clearOut
void clearOut()
Definition: triSurface.C:566
Foam::edgeList
List< edge > edgeList
A List of edges.
Definition: edgeList.H:63
Foam::isoSurfaceBase
Low-level components common to various iso-surface algorithms.
Definition: isoSurfaceBase.H:66
Foam::fileName
A class for handling file names.
Definition: fileName.H:73
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::List::resize
void resize(const label len)
Adjust allocated size of list.
Definition: ListI.H:139
Foam::polyBoundaryMesh
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Definition: polyBoundaryMesh.H:63
Foam::dimLength
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
Foam::fileName::path
static std::string path(const std::string &str)
Return directory path name (part before last /)
Definition: fileNameI.H:176
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::DynamicList
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:55
slicedVolFields.H
Foam::triSurface::triFaceFaces
void triFaceFaces(List< face > &plainFaceList) const
Create a list of faces from the triFaces.
Definition: triSurface.C:723
Foam::SubList
A List obtained as a section of another List.
Definition: SubList.H:54
Foam::GeometricField::primitiveField
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
Definition: GeometricFieldI.H:53
Foam::treeBoundBox
Standard boundBox with extra functionality for use in octree.
Definition: treeBoundBox.H:86
Foam::polyPatch::coupled
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:377
Foam::polyMesh::meshSubDir
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
Definition: polyMesh.H:321
Foam::Time::timeName
static word timeName(const scalar t, const int precision=precision_)
Definition: Time.C:780
Foam::boundaryRegion
The boundaryRegion persistent data saved as a Map<dictionary>.
Definition: boundaryRegion.H:72
Foam::FixedList::size
static constexpr label size() noexcept
Return the number of elements in the FixedList.
Definition: FixedList.H:416
Foam::barycentric2D
Barycentric2D< scalar > barycentric2D
A scalar version of the templated Barycentric2D.
Definition: barycentric2D.H:47
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::coupledPolyPatch
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
Definition: coupledPolyPatch.H:54
defineIsoSurfaceInterpolateMethods
#define defineIsoSurfaceInterpolateMethods(ThisClass)
Definition: isoSurfaceBaseMethods.H:54
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
triSurface.H
Foam::inplaceRenumber
void inplaceRenumber(const labelUList &oldToNew, IntListType &input)
Inplace renumber the values (not the indices) of a list.
Definition: ListOpsTemplates.C:61
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::invertManyToMany
void invertManyToMany(const label len, const UList< InputIntListType > &input, List< OutputIntListType > &output)
Invert many-to-many.
Definition: ListOpsTemplates.C:720
isoSurfaceBaseMethods.H
Convenience macros for instantiating iso-surface interpolate methods.
Foam::geometricSurfacePatchList
List< geometricSurfacePatch > geometricSurfacePatchList
A List of geometricSurfacePatch.
Definition: geometricSurfacePatchList.H:47
Foam::isoSurfacePoint
A surface formed by the iso value. After "Regularised Marching Tetrahedra: improved iso-surface extra...
Definition: isoSurfacePoint.H:87
Foam::isoSurfaceParams::filter
filterType filter() const noexcept
Get current filter type.
Definition: isoSurfaceParams.H:225
Foam::polyMesh::pointsInstance
const fileName & pointsInstance() const
Return the current instance directory for points.
Definition: polyMesh.C:846
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::triPoints
Triangle storage. Null constructable (unfortunately triangle<point, point> is not)
Definition: triPoints.H:52
Foam::Field< scalar >
isoSurfacePoint.H
Foam::labelField
Field< label > labelField
Specialisation of Field<T> for label.
Definition: labelField.H:52
Foam::triSurface::write
void write(Ostream &os) const
Write to Ostream in simple OpenFOAM format.
Definition: triSurfaceIO.C:336
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::collocatedPatch
static bool collocatedPatch(const polyPatch &pp)
Definition: isoSurfacePoint.C:79
Foam::polyPatch
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:68
Foam::storeOp::list
DynamicList< triPoints > & list
Definition: isoSurfacePoint.C:59
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::volScalarField
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
Foam::storeOp::operator()
void operator()(const triPoints &tri)
Definition: isoSurfacePoint.C:66
Foam::FatalError
error FatalError
Foam::isoSurfaceParams
Preferences for controlling iso-surface algorithms.
Definition: isoSurfaceParams.H:107
Foam::Ostream::write
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
Foam::fvMesh
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:85
Foam::calculatedFvPatchField
This boundary condition is not designed to be evaluated; it is assmued that the value is assigned via...
Definition: calculatedFvPatchField.H:66
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::polyPatch::start
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:361
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::volVectorField
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:62
Foam::Pstream::listCombineGather
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
Definition: combineGatherScatter.C:290
Foam::fvMesh::boundary
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:685
Foam::labelListList
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
Foam::primitiveMesh::cellCentres
const vectorField & cellCentres() const
Definition: primitiveMeshCellCentresAndVols.C:84
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Foam::polyPatch::patchSlice
const List< T >::subList patchSlice(const UList< T > &l) const
Slice list to patch.
Definition: polyPatch.H:400
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::VectorSpace< Vector< Cmpt >, Cmpt, 3 >::max
static const Vector< Cmpt > max
Definition: VectorSpace.H:117
f
labelList f(nPoints)
Foam::Vector
Templated 3D Vector derived from VectorSpace adding construction from 3 components,...
Definition: Vector.H:62
Foam::storeOp::storeOp
storeOp(DynamicList< triPoints > &tris)
Definition: isoSurfacePoint.C:61
Foam::UPstream::parRun
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:433
Foam::List< label >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::primitiveMesh::faceCentres
const vectorField & faceCentres() const
Definition: primitiveMeshFaceCentresAndAreas.C:77
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::storeOp::operator()
void operator()(triPoints &&tri)
Definition: isoSurfacePoint.C:71
Foam::sum
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Definition: DimensionedFieldFunctions.C:327
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::boundBox::faceNormals
static const FixedList< vector, 6 > faceNormals
The unit normal per face.
Definition: boundBox.H:92
Foam::SlicedGeometricField
Specialization of GeometricField which holds slices of given complete fields in a form that they act ...
Definition: slicedSurfaceFieldsFwd.H:58
Foam::fvMesh::time
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:280
Foam::face
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:72
Foam::minMax
MinMax< label > minMax(const labelHashSet &set)
Find the min/max values of labelHashSet.
Definition: hashSets.C:61
p0
const volScalarField & p0
Definition: EEqn.H:36
Foam::VectorSpace::size
static constexpr direction size() noexcept
The number of elements in the VectorSpace = Ncmpts.
Definition: VectorSpace.H:176
Foam::Pstream::listCombineScatter
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
Definition: combineGatherScatter.C:432
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::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)
triSurfaceTools.H
Foam::tensor
Tensor< scalar > tensor
Tensor of scalars, i.e. Tensor<scalar>.
Definition: symmTensor.H:61
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
Foam::triPointRef
triangle< point, const point & > triPointRef
A triangle using referred points.
Definition: triangle.H:71
Foam::UPstream::commsTypes::blocking