isoSurface.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "isoSurface.H"
27 #include "fvMesh.H"
28 #include "mergePoints.H"
30 #include "slicedVolFields.H"
31 #include "volFields.H"
32 #include "surfaceFields.H"
33 #include "OFstream.H"
34 #include "meshTools.H"
35 
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 
38 namespace Foam
39 {
40 defineTypeNameAndDebug(isoSurface, 0);
41 }
42 
43 
44 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
45 
46 bool Foam::isoSurface::noTransform(const tensor& tt) const
47 {
48  return
49  (mag(tt.xx()-1) < mergeDistance_)
50  && (mag(tt.yy()-1) < mergeDistance_)
51  && (mag(tt.zz()-1) < mergeDistance_)
52  && (mag(tt.xy()) < mergeDistance_)
53  && (mag(tt.xz()) < mergeDistance_)
54  && (mag(tt.yx()) < mergeDistance_)
55  && (mag(tt.yz()) < mergeDistance_)
56  && (mag(tt.zx()) < mergeDistance_)
57  && (mag(tt.zy()) < mergeDistance_);
58 }
59 
60 
61 // Calculates per face whether couple is collocated.
62 bool Foam::isoSurface::collocatedPatch(const polyPatch& pp)
63 {
64  const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>(pp);
65  return cpp.parallel() && !cpp.separated();
66 }
67 
68 
69 // Calculates per face whether couple is collocated.
70 Foam::PackedBoolList Foam::isoSurface::collocatedFaces
71 (
72  const coupledPolyPatch& pp
73 ) const
74 {
75  // Initialise to false
76  PackedBoolList collocated(pp.size());
77 
78  if (isA<processorPolyPatch>(pp))
79  {
80  if (collocatedPatch(pp))
81  {
82  forAll(pp, i)
83  {
84  collocated[i] = 1u;
85  }
86  }
87  }
88  else if (isA<cyclicPolyPatch>(pp))
89  {
90  if (collocatedPatch(pp))
91  {
92  forAll(pp, i)
93  {
94  collocated[i] = 1u;
95  }
96  }
97  }
98  else
99  {
101  (
102  "isoSurface::collocatedFaces(const coupledPolyPatch&) const"
103  ) << "Unhandled coupledPolyPatch type " << pp.type()
104  << abort(FatalError);
105  }
106  return collocated;
107 }
108 
109 
110 void Foam::isoSurface::syncUnseparatedPoints
111 (
112  pointField& pointValues,
113  const point& nullValue
114 ) const
115 {
116  // Until syncPointList handles separated coupled patches with multiple
117  // transforms do our own synchronisation of non-separated patches only
118  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
119 
120  if (Pstream::parRun())
121  {
122  // Send
123  forAll(patches, patchI)
124  {
125  if
126  (
127  isA<processorPolyPatch>(patches[patchI])
128  && patches[patchI].nPoints() > 0
129  && collocatedPatch(patches[patchI])
130  )
131  {
132  const processorPolyPatch& pp =
133  refCast<const processorPolyPatch>(patches[patchI]);
134 
135  const labelList& meshPts = pp.meshPoints();
136  const labelList& nbrPts = pp.neighbPoints();
137 
138  pointField patchInfo(meshPts.size());
139 
140  forAll(nbrPts, pointI)
141  {
142  label nbrPointI = nbrPts[pointI];
143  patchInfo[nbrPointI] = pointValues[meshPts[pointI]];
144  }
145 
146  OPstream toNbr(Pstream::blocking, pp.neighbProcNo());
147  toNbr << patchInfo;
148  }
149  }
150 
151  // Receive and combine.
152 
153  forAll(patches, patchI)
154  {
155  if
156  (
157  isA<processorPolyPatch>(patches[patchI])
158  && patches[patchI].nPoints() > 0
159  && collocatedPatch(patches[patchI])
160  )
161  {
162  const processorPolyPatch& pp =
163  refCast<const processorPolyPatch>(patches[patchI]);
164 
165  pointField nbrPatchInfo(pp.nPoints());
166  {
167  // We do not know the number of points on the other side
168  // so cannot use Pstream::read.
169  IPstream fromNbr(Pstream::blocking, pp.neighbProcNo());
170  fromNbr >> nbrPatchInfo;
171  }
172 
173  const labelList& meshPts = pp.meshPoints();
174 
175  forAll(meshPts, pointI)
176  {
177  label meshPointI = meshPts[pointI];
178  minEqOp<point>()
179  (
180  pointValues[meshPointI],
181  nbrPatchInfo[pointI]
182  );
183  }
184  }
185  }
186  }
187 
188  // Do the cyclics.
189  forAll(patches, patchI)
190  {
191  if (isA<cyclicPolyPatch>(patches[patchI]))
192  {
193  const cyclicPolyPatch& cycPatch =
194  refCast<const cyclicPolyPatch>(patches[patchI]);
195 
196  if (cycPatch.owner() && collocatedPatch(cycPatch))
197  {
198  // Owner does all.
199 
200  const edgeList& coupledPoints = cycPatch.coupledPoints();
201  const labelList& meshPts = cycPatch.meshPoints();
202  const cyclicPolyPatch& nbrPatch = cycPatch.neighbPatch();
203  const labelList& nbrMeshPoints = nbrPatch.meshPoints();
204 
205  pointField half0Values(coupledPoints.size());
206  pointField half1Values(coupledPoints.size());
207 
208  forAll(coupledPoints, i)
209  {
210  const edge& e = coupledPoints[i];
211  half0Values[i] = pointValues[meshPts[e[0]]];
212  half1Values[i] = pointValues[nbrMeshPoints[e[1]]];
213  }
214 
215  forAll(coupledPoints, i)
216  {
217  const edge& e = coupledPoints[i];
218  label p0 = meshPts[e[0]];
219  label p1 = nbrMeshPoints[e[1]];
220 
221  minEqOp<point>()(pointValues[p0], half1Values[i]);
222  minEqOp<point>()(pointValues[p1], half0Values[i]);
223  }
224  }
225  }
226  }
227 
228  // Synchronize multiple shared points.
229  const globalMeshData& pd = mesh_.globalData();
230 
231  if (pd.nGlobalPoints() > 0)
232  {
233  // Values on shared points.
234  pointField sharedPts(pd.nGlobalPoints(), nullValue);
235 
236  forAll(pd.sharedPointLabels(), i)
237  {
238  label meshPointI = pd.sharedPointLabels()[i];
239  // Fill my entries in the shared points
240  sharedPts[pd.sharedPointAddr()[i]] = pointValues[meshPointI];
241  }
242 
243  // Combine on master.
244  Pstream::listCombineGather(sharedPts, minEqOp<point>());
245  Pstream::listCombineScatter(sharedPts);
246 
247  // Now we will all have the same information. Merge it back with
248  // my local information.
249  forAll(pd.sharedPointLabels(), i)
250  {
251  label meshPointI = pd.sharedPointLabels()[i];
252  pointValues[meshPointI] = sharedPts[pd.sharedPointAddr()[i]];
253  }
254  }
255 }
256 
257 
258 Foam::scalar Foam::isoSurface::isoFraction
259 (
260  const scalar s0,
261  const scalar s1
262 ) const
263 {
264  scalar d = s1-s0;
265 
266  if (mag(d) > VSMALL)
267  {
268  return (iso_-s0)/d;
269  }
270  else
271  {
272  return -1.0;
273  }
274 }
275 
276 
277 bool Foam::isoSurface::isEdgeOfFaceCut
278 (
279  const scalarField& pVals,
280  const face& f,
281  const bool ownLower,
282  const bool neiLower
283 ) const
284 {
285  forAll(f, fp)
286  {
287  bool fpLower = (pVals[f[fp]] < iso_);
288  if
289  (
290  (fpLower != ownLower)
291  || (fpLower != neiLower)
292  || (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
293  )
294  {
295  return true;
296  }
297  }
298  return false;
299 }
300 
301 
302 // Get neighbour value and position.
303 void Foam::isoSurface::getNeighbour
304 (
305  const labelList& boundaryRegion,
306  const volVectorField& meshC,
307  const volScalarField& cVals,
308  const label cellI,
309  const label faceI,
310  scalar& nbrValue,
311  point& nbrPoint
312 ) const
313 {
314  const labelList& own = mesh_.faceOwner();
315  const labelList& nei = mesh_.faceNeighbour();
316 
317  if (mesh_.isInternalFace(faceI))
318  {
319  label nbr = (own[faceI] == cellI ? nei[faceI] : own[faceI]);
320  nbrValue = cVals[nbr];
321  nbrPoint = meshC[nbr];
322  }
323  else
324  {
325  label bFaceI = faceI-mesh_.nInternalFaces();
326  label patchI = boundaryRegion[bFaceI];
327  const polyPatch& pp = mesh_.boundaryMesh()[patchI];
328  label patchFaceI = faceI-pp.start();
329 
330  nbrValue = cVals.boundaryField()[patchI][patchFaceI];
331  nbrPoint = meshC.boundaryField()[patchI][patchFaceI];
332  }
333 }
334 
335 
336 // Determine for every face/cell whether it (possibly) generates triangles.
337 void Foam::isoSurface::calcCutTypes
338 (
339  const labelList& boundaryRegion,
340  const volVectorField& meshC,
341  const volScalarField& cVals,
342  const scalarField& pVals
343 )
344 {
345  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
346  const labelList& own = mesh_.faceOwner();
347  const labelList& nei = mesh_.faceNeighbour();
348 
349  faceCutType_.setSize(mesh_.nFaces());
350  faceCutType_ = NOTCUT;
351 
352  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
353  {
354  // CC edge.
355  bool ownLower = (cVals[own[faceI]] < iso_);
356 
357  scalar nbrValue;
358  point nbrPoint;
359  getNeighbour
360  (
361  boundaryRegion,
362  meshC,
363  cVals,
364  own[faceI],
365  faceI,
366  nbrValue,
367  nbrPoint
368  );
369 
370  bool neiLower = (nbrValue < iso_);
371 
372  if (ownLower != neiLower)
373  {
374  faceCutType_[faceI] = CUT;
375  }
376  else
377  {
378  // See if any mesh edge is cut by looping over all the edges of the
379  // face.
380  const face f = mesh_.faces()[faceI];
381 
382  if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
383  {
384  faceCutType_[faceI] = CUT;
385  }
386  }
387  }
388 
389  forAll(patches, patchI)
390  {
391  const polyPatch& pp = patches[patchI];
392 
393  label faceI = pp.start();
394 
395  forAll(pp, i)
396  {
397  bool ownLower = (cVals[own[faceI]] < iso_);
398 
399  scalar nbrValue;
400  point nbrPoint;
401  getNeighbour
402  (
403  boundaryRegion,
404  meshC,
405  cVals,
406  own[faceI],
407  faceI,
408  nbrValue,
409  nbrPoint
410  );
411 
412  bool neiLower = (nbrValue < iso_);
413 
414  if (ownLower != neiLower)
415  {
416  faceCutType_[faceI] = CUT;
417  }
418  else
419  {
420  // Mesh edge.
421  const face f = mesh_.faces()[faceI];
422 
423  if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
424  {
425  faceCutType_[faceI] = CUT;
426  }
427  }
428 
429  faceI++;
430  }
431  }
432 
433 
434 
435  nCutCells_ = 0;
436  cellCutType_.setSize(mesh_.nCells());
437  cellCutType_ = NOTCUT;
438 
439  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
440  {
441  if (faceCutType_[faceI] != NOTCUT)
442  {
443  if (cellCutType_[own[faceI]] == NOTCUT)
444  {
445  cellCutType_[own[faceI]] = CUT;
446  nCutCells_++;
447  }
448  if (cellCutType_[nei[faceI]] == NOTCUT)
449  {
450  cellCutType_[nei[faceI]] = CUT;
451  nCutCells_++;
452  }
453  }
454  }
455  for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
456  {
457  if (faceCutType_[faceI] != NOTCUT)
458  {
459  if (cellCutType_[own[faceI]] == NOTCUT)
460  {
461  cellCutType_[own[faceI]] = CUT;
462  nCutCells_++;
463  }
464  }
465  }
466 
467  if (debug)
468  {
469  Pout<< "isoSurface : detected " << nCutCells_
470  << " candidate cut cells (out of " << mesh_.nCells()
471  << ")." << endl;
472  }
473 }
474 
475 
476 // Return the two common points between two triangles
477 Foam::labelPair Foam::isoSurface::findCommonPoints
478 (
479  const labelledTri& tri0,
480  const labelledTri& tri1
481 )
482 {
483  labelPair common(-1, -1);
484 
485  label fp0 = 0;
486  label fp1 = findIndex(tri1, tri0[fp0]);
487 
488  if (fp1 == -1)
489  {
490  fp0 = 1;
491  fp1 = findIndex(tri1, tri0[fp0]);
492  }
493 
494  if (fp1 != -1)
495  {
496  // So tri0[fp0] is tri1[fp1]
497 
498  // Find next common point
499  label fp0p1 = tri0.fcIndex(fp0);
500  label fp1p1 = tri1.fcIndex(fp1);
501  label fp1m1 = tri1.rcIndex(fp1);
502 
503  if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
504  {
505  common[0] = tri0[fp0];
506  common[1] = tri0[fp0p1];
507  }
508  }
509  return common;
510 }
511 
512 
513 // Caculate centre of surface.
514 Foam::point Foam::isoSurface::calcCentre(const triSurface& s)
515 {
516  vector sum = vector::zero;
517 
518  forAll(s, i)
519  {
520  sum += s[i].centre(s.points());
521  }
522  return sum/s.size();
523 }
524 
525 
526 // Replace surface (localPoints, localTris) with single point. Returns
527 // point. Destructs arguments.
528 Foam::pointIndexHit Foam::isoSurface::collapseSurface
529 (
530  pointField& localPoints,
531  DynamicList<labelledTri, 64>& localTris
532 )
533 {
534  pointIndexHit info(false, vector::zero, localTris.size());
535 
536  if (localTris.size() == 1)
537  {
538  const labelledTri& tri = localTris[0];
539  info.setPoint(tri.centre(localPoints));
540  info.setHit();
541  }
542  else if (localTris.size() == 2)
543  {
544  // Check if the two triangles share an edge.
545  const labelledTri& tri0 = localTris[0];
546  const labelledTri& tri1 = localTris[0];
547 
548  labelPair shared = findCommonPoints(tri0, tri1);
549 
550  if (shared[0] != -1)
551  {
552  info.setPoint
553  (
554  0.5
555  * (
556  tri0.centre(localPoints)
557  + tri1.centre(localPoints)
558  )
559  );
560  info.setHit();
561  }
562  }
563  else if (localTris.size())
564  {
565  // Check if single region. Rare situation.
566  triSurface surf
567  (
568  localTris,
570  localPoints,
571  true
572  );
573  localTris.clearStorage();
574 
575  labelList faceZone;
576  label nZones = surf.markZones
577  (
578  boolList(surf.nEdges(), false),
579  faceZone
580  );
581 
582  if (nZones == 1)
583  {
584  info.setPoint(calcCentre(surf));
585  info.setHit();
586  }
587  }
588 
589  return info;
590 }
591 
592 
593 // Determine per cell centre whether all the intersections get collapsed
594 // to a single point
595 void Foam::isoSurface::calcSnappedCc
596 (
597  const labelList& boundaryRegion,
598  const volVectorField& meshC,
599  const volScalarField& cVals,
600  const scalarField& pVals,
601 
602  DynamicList<point>& snappedPoints,
603  labelList& snappedCc
604 ) const
605 {
606  const pointField& pts = mesh_.points();
607  const pointField& cc = mesh_.cellCentres();
608 
609  snappedCc.setSize(mesh_.nCells());
610  snappedCc = -1;
611 
612  // Work arrays
613  DynamicList<point, 64> localTriPoints(64);
614 
615  forAll(mesh_.cells(), cellI)
616  {
617  if (cellCutType_[cellI] == CUT)
618  {
619  scalar cVal = cVals[cellI];
620 
621  const cell& cFaces = mesh_.cells()[cellI];
622 
623  localTriPoints.clear();
624  label nOther = 0;
625  point otherPointSum = vector::zero;
626 
627  // Create points for all intersections close to cell centre
628  // (i.e. from pyramid edges)
629 
630  forAll(cFaces, cFaceI)
631  {
632  label faceI = cFaces[cFaceI];
633 
634  scalar nbrValue;
635  point nbrPoint;
636  getNeighbour
637  (
638  boundaryRegion,
639  meshC,
640  cVals,
641  cellI,
642  faceI,
643  nbrValue,
644  nbrPoint
645  );
646 
647  // Calculate intersection points of edges to cell centre
648  FixedList<scalar, 3> s;
649  FixedList<point, 3> pt;
650 
651  // From cc to neighbour cc.
652  s[2] = isoFraction(cVal, nbrValue);
653  pt[2] = (1.0-s[2])*cc[cellI] + s[2]*nbrPoint;
654 
655  const face& f = mesh_.faces()[cFaces[cFaceI]];
656 
657  forAll(f, fp)
658  {
659  // From cc to fp
660  label p0 = f[fp];
661  s[0] = isoFraction(cVal, pVals[p0]);
662  pt[0] = (1.0-s[0])*cc[cellI] + s[0]*pts[p0];
663 
664  // From cc to fp+1
665  label p1 = f[f.fcIndex(fp)];
666  s[1] = isoFraction(cVal, pVals[p1]);
667  pt[1] = (1.0-s[1])*cc[cellI] + s[1]*pts[p1];
668 
669  if
670  (
671  (s[0] >= 0.0 && s[0] <= 0.5)
672  && (s[1] >= 0.0 && s[1] <= 0.5)
673  && (s[2] >= 0.0 && s[2] <= 0.5)
674  )
675  {
676  localTriPoints.append(pt[0]);
677  localTriPoints.append(pt[1]);
678  localTriPoints.append(pt[2]);
679  }
680  else
681  {
682  // Get average of all other points
683  forAll(s, i)
684  {
685  if (s[i] >= 0.0 && s[i] <= 0.5)
686  {
687  otherPointSum += pt[i];
688  nOther++;
689  }
690  }
691  }
692  }
693  }
694 
695  if (localTriPoints.size() == 0)
696  {
697  // No complete triangles. Get average of all intersection
698  // points.
699  if (nOther > 0)
700  {
701  snappedCc[cellI] = snappedPoints.size();
702  snappedPoints.append(otherPointSum/nOther);
703 
704  //Pout<< " point:" << pointI
705  // << " replacing coord:" << mesh_.points()[pointI]
706  // << " by average:" << collapsedPoint[pointI] << endl;
707  }
708  }
709  else if (localTriPoints.size() == 3)
710  {
711  // Single triangle. No need for any analysis. Average points.
713  points.transfer(localTriPoints);
714  snappedCc[cellI] = snappedPoints.size();
715  snappedPoints.append(sum(points)/points.size());
716 
717  //Pout<< " point:" << pointI
718  // << " replacing coord:" << mesh_.points()[pointI]
719  // << " by average:" << collapsedPoint[pointI] << endl;
720  }
721  else
722  {
723  // Convert points into triSurface.
724 
725  // Merge points and compact out non-valid triangles
726  labelList triMap; // merged to unmerged triangle
727  labelList triPointReverseMap; // unmerged to merged point
728  triSurface surf
729  (
730  stitchTriPoints
731  (
732  false, // do not check for duplicate tris
733  localTriPoints,
734  triPointReverseMap,
735  triMap
736  )
737  );
738 
739  labelList faceZone;
740  label nZones = surf.markZones
741  (
742  boolList(surf.nEdges(), false),
743  faceZone
744  );
745 
746  if (nZones == 1)
747  {
748  snappedCc[cellI] = snappedPoints.size();
749  snappedPoints.append(calcCentre(surf));
750  //Pout<< " point:" << pointI << " nZones:" << nZones
751  // << " replacing coord:" << mesh_.points()[pointI]
752  // << " by average:" << collapsedPoint[pointI] << endl;
753  }
754  }
755  }
756  }
757 }
758 
759 
760 // Determine per meshpoint whether all the intersections get collapsed
761 // to a single point
762 void Foam::isoSurface::calcSnappedPoint
763 (
764  const PackedBoolList& isBoundaryPoint,
765  const labelList& boundaryRegion,
766  const volVectorField& meshC,
767  const volScalarField& cVals,
768  const scalarField& pVals,
769 
770  DynamicList<point>& snappedPoints,
771  labelList& snappedPoint
772 ) const
773 {
774  const pointField& pts = mesh_.points();
775  const pointField& cc = mesh_.cellCentres();
776 
777  pointField collapsedPoint(mesh_.nPoints(), point::max);
778 
779 
780  // Work arrays
781  DynamicList<point, 64> localTriPoints(100);
782 
783  forAll(mesh_.pointFaces(), pointI)
784  {
785  if (isBoundaryPoint.get(pointI) == 1)
786  {
787  continue;
788  }
789 
790  const labelList& pFaces = mesh_.pointFaces()[pointI];
791 
792  bool anyCut = false;
793 
794  forAll(pFaces, i)
795  {
796  label faceI = pFaces[i];
797 
798  if (faceCutType_[faceI] == CUT)
799  {
800  anyCut = true;
801  break;
802  }
803  }
804 
805  if (!anyCut)
806  {
807  continue;
808  }
809 
810 
811  localTriPoints.clear();
812  label nOther = 0;
813  point otherPointSum = vector::zero;
814 
815  forAll(pFaces, pFaceI)
816  {
817  // Create points for all intersections close to point
818  // (i.e. from pyramid edges)
819 
820  label faceI = pFaces[pFaceI];
821  const face& f = mesh_.faces()[faceI];
822  label own = mesh_.faceOwner()[faceI];
823 
824  // Get neighbour value
825  scalar nbrValue;
826  point nbrPoint;
827  getNeighbour
828  (
829  boundaryRegion,
830  meshC,
831  cVals,
832  own,
833  faceI,
834  nbrValue,
835  nbrPoint
836  );
837 
838  // Calculate intersection points of edges emanating from point
839  FixedList<scalar, 4> s;
840  FixedList<point, 4> pt;
841 
842  label fp = findIndex(f, pointI);
843  s[0] = isoFraction(pVals[pointI], cVals[own]);
844  pt[0] = (1.0-s[0])*pts[pointI] + s[0]*cc[own];
845 
846  s[1] = isoFraction(pVals[pointI], nbrValue);
847  pt[1] = (1.0-s[1])*pts[pointI] + s[1]*nbrPoint;
848 
849  label nextPointI = f[f.fcIndex(fp)];
850  s[2] = isoFraction(pVals[pointI], pVals[nextPointI]);
851  pt[2] = (1.0-s[2])*pts[pointI] + s[2]*pts[nextPointI];
852 
853  label prevPointI = f[f.rcIndex(fp)];
854  s[3] = isoFraction(pVals[pointI], pVals[prevPointI]);
855  pt[3] = (1.0-s[3])*pts[pointI] + s[3]*pts[prevPointI];
856 
857  if
858  (
859  (s[0] >= 0.0 && s[0] <= 0.5)
860  && (s[1] >= 0.0 && s[1] <= 0.5)
861  && (s[2] >= 0.0 && s[2] <= 0.5)
862  )
863  {
864  localTriPoints.append(pt[0]);
865  localTriPoints.append(pt[1]);
866  localTriPoints.append(pt[2]);
867  }
868  if
869  (
870  (s[0] >= 0.0 && s[0] <= 0.5)
871  && (s[1] >= 0.0 && s[1] <= 0.5)
872  && (s[3] >= 0.0 && s[3] <= 0.5)
873  )
874  {
875  localTriPoints.append(pt[3]);
876  localTriPoints.append(pt[0]);
877  localTriPoints.append(pt[1]);
878  }
879 
880  // Get average of points as fallback
881  forAll(s, i)
882  {
883  if (s[i] >= 0.0 && s[i] <= 0.5)
884  {
885  otherPointSum += pt[i];
886  nOther++;
887  }
888  }
889  }
890 
891  if (localTriPoints.size() == 0)
892  {
893  // No complete triangles. Get average of all intersection
894  // points.
895  if (nOther > 0)
896  {
897  collapsedPoint[pointI] = otherPointSum/nOther;
898  }
899  }
900  else if (localTriPoints.size() == 3)
901  {
902  // Single triangle. No need for any analysis. Average points.
904  points.transfer(localTriPoints);
905  collapsedPoint[pointI] = sum(points)/points.size();
906  }
907  else
908  {
909  // Convert points into triSurface.
910 
911  // Merge points and compact out non-valid triangles
912  labelList triMap; // merged to unmerged triangle
913  labelList triPointReverseMap; // unmerged to merged point
914  triSurface surf
915  (
916  stitchTriPoints
917  (
918  false, // do not check for duplicate tris
919  localTriPoints,
920  triPointReverseMap,
921  triMap
922  )
923  );
924 
925  labelList faceZone;
926  label nZones = surf.markZones
927  (
928  boolList(surf.nEdges(), false),
929  faceZone
930  );
931 
932  if (nZones == 1)
933  {
934  collapsedPoint[pointI] = calcCentre(surf);
935  }
936  }
937  }
938 
939 
940  // Synchronise snap location
941  syncUnseparatedPoints(collapsedPoint, point::max);
942 
943 
944  snappedPoint.setSize(mesh_.nPoints());
945  snappedPoint = -1;
946 
947  forAll(collapsedPoint, pointI)
948  {
949  if (collapsedPoint[pointI] != point::max)
950  {
951  snappedPoint[pointI] = snappedPoints.size();
952  snappedPoints.append(collapsedPoint[pointI]);
953  }
954  }
955 }
956 
957 
958 Foam::triSurface Foam::isoSurface::stitchTriPoints
959 (
960  const bool checkDuplicates,
961  const List<point>& triPoints,
962  labelList& triPointReverseMap, // unmerged to merged point
963  labelList& triMap // merged to unmerged triangle
964 ) const
965 {
966  label nTris = triPoints.size()/3;
967 
968  if ((triPoints.size() % 3) != 0)
969  {
970  FatalErrorIn("isoSurface::stitchTriPoints(..)")
971  << "Problem: number of points " << triPoints.size()
972  << " not a multiple of 3." << abort(FatalError);
973  }
974 
975  pointField newPoints;
977  (
978  triPoints,
979  mergeDistance_,
980  false,
981  triPointReverseMap,
982  newPoints
983  );
984 
985  // Check that enough merged.
986  if (debug)
987  {
988  pointField newNewPoints;
989  labelList oldToNew;
990  bool hasMerged = mergePoints
991  (
992  newPoints,
993  mergeDistance_,
994  true,
995  oldToNew,
996  newNewPoints
997  );
998 
999  if (hasMerged)
1000  {
1001  FatalErrorIn("isoSurface::stitchTriPoints(..)")
1002  << "Merged points contain duplicates"
1003  << " when merging with distance " << mergeDistance_ << endl
1004  << "merged:" << newPoints.size() << " re-merged:"
1005  << newNewPoints.size()
1006  << abort(FatalError);
1007  }
1008  }
1009 
1010 
1011  List<labelledTri> tris;
1012  {
1013  DynamicList<labelledTri> dynTris(nTris);
1014  label rawPointI = 0;
1015  DynamicList<label> newToOldTri(nTris);
1016 
1017  for (label oldTriI = 0; oldTriI < nTris; oldTriI++)
1018  {
1019  labelledTri tri
1020  (
1021  triPointReverseMap[rawPointI],
1022  triPointReverseMap[rawPointI+1],
1023  triPointReverseMap[rawPointI+2],
1024  0
1025  );
1026  rawPointI += 3;
1027 
1028  if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
1029  {
1030  newToOldTri.append(oldTriI);
1031  dynTris.append(tri);
1032  }
1033  }
1034 
1035  triMap.transfer(newToOldTri);
1036  tris.transfer(dynTris);
1037  }
1038 
1039 
1040 
1041  // Determine 'flat hole' situation (see RMT paper).
1042  // Two unconnected triangles get connected because (some of) the edges
1043  // separating them get collapsed. Below only checks for duplicate triangles,
1044  // not non-manifold edge connectivity.
1045  if (checkDuplicates)
1046  {
1047  labelListList pointFaces;
1048  invertManyToMany(newPoints.size(), tris, pointFaces);
1049 
1050  // Filter out duplicates.
1051  DynamicList<label> newToOldTri(tris.size());
1052 
1053  forAll(tris, triI)
1054  {
1055  const labelledTri& tri = tris[triI];
1056  const labelList& pFaces = pointFaces[tri[0]];
1057 
1058  // Find the maximum of any duplicates. Maximum since the tris
1059  // below triI
1060  // get overwritten so we cannot use them in a comparison.
1061  label dupTriI = -1;
1062  forAll(pFaces, i)
1063  {
1064  label nbrTriI = pFaces[i];
1065 
1066  if (nbrTriI > triI && (tris[nbrTriI] == tri))
1067  {
1068  //Pout<< "Duplicate : " << triI << " verts:" << tri
1069  // << " to " << nbrTriI << " verts:" << tris[nbrTriI]
1070  // << endl;
1071  dupTriI = nbrTriI;
1072  break;
1073  }
1074  }
1075 
1076  if (dupTriI == -1)
1077  {
1078  // There is no (higher numbered) duplicate triangle
1079  label newTriI = newToOldTri.size();
1080  newToOldTri.append(triMap[triI]);
1081  tris[newTriI] = tris[triI];
1082  }
1083  }
1084 
1085  triMap.transfer(newToOldTri);
1086  tris.setSize(triMap.size());
1087 
1088  if (debug)
1089  {
1090  Pout<< "isoSurface : merged from " << nTris
1091  << " down to " << tris.size() << " unique triangles." << endl;
1092  }
1093 
1094 
1095  if (debug)
1096  {
1097  triSurface surf(tris, geometricSurfacePatchList(0), newPoints);
1098 
1099  forAll(surf, faceI)
1100  {
1101  const labelledTri& f = surf[faceI];
1102  const labelList& fFaces = surf.faceFaces()[faceI];
1103 
1104  forAll(fFaces, i)
1105  {
1106  label nbrFaceI = fFaces[i];
1107 
1108  if (nbrFaceI <= faceI)
1109  {
1110  // lower numbered faces already checked
1111  continue;
1112  }
1113 
1114  const labelledTri& nbrF = surf[nbrFaceI];
1115 
1116  if (f == nbrF)
1117  {
1118  FatalErrorIn("validTri(const triSurface&, const label)")
1119  << "Check : "
1120  << " triangle " << faceI << " vertices " << f
1121  << " fc:" << f.centre(surf.points())
1122  << " has the same vertices as triangle " << nbrFaceI
1123  << " vertices " << nbrF
1124  << " fc:" << nbrF.centre(surf.points())
1125  << abort(FatalError);
1126  }
1127  }
1128  }
1129  }
1130  }
1131 
1132  return triSurface(tris, geometricSurfacePatchList(0), newPoints, true);
1133 }
1134 
1135 
1136 // Does face use valid vertices?
1137 bool Foam::isoSurface::validTri(const triSurface& surf, const label faceI)
1138 {
1139  // Simple check on indices ok.
1140 
1141  const labelledTri& f = surf[faceI];
1142 
1143  if
1144  (
1145  (f[0] < 0) || (f[0] >= surf.points().size())
1146  || (f[1] < 0) || (f[1] >= surf.points().size())
1147  || (f[2] < 0) || (f[2] >= surf.points().size())
1148  )
1149  {
1150  WarningIn("validTri(const triSurface&, const label)")
1151  << "triangle " << faceI << " vertices " << f
1152  << " uses point indices outside point range 0.."
1153  << surf.points().size()-1 << endl;
1154 
1155  return false;
1156  }
1157 
1158  if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
1159  {
1160  WarningIn("validTri(const triSurface&, const label)")
1161  << "triangle " << faceI
1162  << " uses non-unique vertices " << f
1163  << endl;
1164  return false;
1165  }
1166 
1167  // duplicate triangle check
1168 
1169  const labelList& fFaces = surf.faceFaces()[faceI];
1170 
1171  // Check if faceNeighbours use same points as this face.
1172  // Note: discards normal information - sides of baffle are merged.
1173  forAll(fFaces, i)
1174  {
1175  label nbrFaceI = fFaces[i];
1176 
1177  if (nbrFaceI <= faceI)
1178  {
1179  // lower numbered faces already checked
1180  continue;
1181  }
1182 
1183  const labelledTri& nbrF = surf[nbrFaceI];
1184 
1185  if
1186  (
1187  ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2]))
1188  && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2]))
1189  && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
1190  )
1191  {
1192  WarningIn("validTri(const triSurface&, const label)")
1193  << "triangle " << faceI << " vertices " << f
1194  << " fc:" << f.centre(surf.points())
1195  << " has the same vertices as triangle " << nbrFaceI
1196  << " vertices " << nbrF
1197  << " fc:" << nbrF.centre(surf.points())
1198  << endl;
1199 
1200  return false;
1201  }
1202  }
1203  return true;
1204 }
1205 
1206 
1207 void Foam::isoSurface::calcAddressing
1208 (
1209  const triSurface& surf,
1210  List<FixedList<label, 3> >& faceEdges,
1211  labelList& edgeFace0,
1212  labelList& edgeFace1,
1213  Map<labelList>& edgeFacesRest
1214 ) const
1215 {
1216  const pointField& points = surf.points();
1217 
1218  pointField edgeCentres(3*surf.size());
1219  label edgeI = 0;
1220  forAll(surf, triI)
1221  {
1222  const labelledTri& tri = surf[triI];
1223  edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
1224  edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
1225  edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
1226  }
1227 
1228  pointField mergedCentres;
1229  labelList oldToMerged;
1230  bool hasMerged = mergePoints
1231  (
1232  edgeCentres,
1233  mergeDistance_,
1234  false,
1235  oldToMerged,
1236  mergedCentres
1237  );
1238 
1239  if (debug)
1240  {
1241  Pout<< "isoSurface : detected "
1242  << mergedCentres.size()
1243  << " geometric edges on " << surf.size() << " triangles." << endl;
1244  }
1245 
1246  if (!hasMerged)
1247  {
1248  return;
1249  }
1250 
1251 
1252  // Determine faceEdges
1253  faceEdges.setSize(surf.size());
1254  edgeI = 0;
1255  forAll(surf, triI)
1256  {
1257  faceEdges[triI][0] = oldToMerged[edgeI++];
1258  faceEdges[triI][1] = oldToMerged[edgeI++];
1259  faceEdges[triI][2] = oldToMerged[edgeI++];
1260  }
1261 
1262 
1263  // Determine edgeFaces
1264  edgeFace0.setSize(mergedCentres.size());
1265  edgeFace0 = -1;
1266  edgeFace1.setSize(mergedCentres.size());
1267  edgeFace1 = -1;
1268  edgeFacesRest.clear();
1269 
1270  // Overflow edge faces for geometric shared edges that turned
1271  // out to be different anyway.
1272  EdgeMap<labelList> extraEdgeFaces(mergedCentres.size()/100);
1273 
1274  forAll(oldToMerged, oldEdgeI)
1275  {
1276  label triI = oldEdgeI / 3;
1277  label edgeI = oldToMerged[oldEdgeI];
1278 
1279  if (edgeFace0[edgeI] == -1)
1280  {
1281  // First triangle for edge
1282  edgeFace0[edgeI] = triI;
1283  }
1284  else
1285  {
1286  //- Check that the two triangles actually topologically
1287  // share an edge
1288  const labelledTri& prevTri = surf[edgeFace0[edgeI]];
1289  const labelledTri& tri = surf[triI];
1290 
1291  label fp = oldEdgeI % 3;
1292 
1293  edge e(tri[fp], tri[tri.fcIndex(fp)]);
1294 
1295  label prevTriIndex = -1;
1296 
1297  forAll(prevTri, i)
1298  {
1299  if (edge(prevTri[i], prevTri[prevTri.fcIndex(i)]) == e)
1300  {
1301  prevTriIndex = i;
1302  break;
1303  }
1304  }
1305 
1306  if (prevTriIndex == -1)
1307  {
1308  // Different edge. Store for later.
1309  EdgeMap<labelList>::iterator iter = extraEdgeFaces.find(e);
1310 
1311  if (iter != extraEdgeFaces.end())
1312  {
1313  labelList& eFaces = iter();
1314  label sz = eFaces.size();
1315  eFaces.setSize(sz+1);
1316  eFaces[sz] = triI;
1317  }
1318  else
1319  {
1320  extraEdgeFaces.insert(e, labelList(1, triI));
1321  }
1322  }
1323  else if (edgeFace1[edgeI] == -1)
1324  {
1325  edgeFace1[edgeI] = triI;
1326  }
1327  else
1328  {
1329  //WarningIn("orientSurface(triSurface&)")
1330  // << "Edge " << edgeI << " with centre "
1331  // << mergedCentres[edgeI]
1332  // << " used by more than two triangles: "
1333  // << edgeFace0[edgeI] << ", "
1334  // << edgeFace1[edgeI] << " and " << triI << endl;
1335  Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);
1336 
1337  if (iter != edgeFacesRest.end())
1338  {
1339  labelList& eFaces = iter();
1340  label sz = eFaces.size();
1341  eFaces.setSize(sz+1);
1342  eFaces[sz] = triI;
1343  }
1344  else
1345  {
1346  edgeFacesRest.insert(edgeI, labelList(1, triI));
1347  }
1348  }
1349  }
1350  }
1351 
1352  // Add extraEdgeFaces
1353  edgeI = edgeFace0.size();
1354 
1355  edgeFace0.setSize(edgeI + extraEdgeFaces.size());
1356  edgeFace1.setSize(edgeI + extraEdgeFaces.size(), -1);
1357 
1358  forAllConstIter(EdgeMap<labelList>, extraEdgeFaces, iter)
1359  {
1360  const labelList& eFaces = iter();
1361 
1362  // The current edge will become edgeI. Replace all occurrences in
1363  // faceEdges
1364  forAll(eFaces, i)
1365  {
1366  label triI = eFaces[i];
1367  const labelledTri& tri = surf[triI];
1368 
1369  FixedList<label, 3>& fEdges = faceEdges[triI];
1370  forAll(tri, fp)
1371  {
1372  edge e(tri[fp], tri[tri.fcIndex(fp)]);
1373 
1374  if (e == iter.key())
1375  {
1376  fEdges[fp] = edgeI;
1377  break;
1378  }
1379  }
1380  }
1381 
1382 
1383  // Add face to edgeFaces
1384 
1385  edgeFace0[edgeI] = eFaces[0];
1386 
1387  if (eFaces.size() >= 2)
1388  {
1389  edgeFace1[edgeI] = eFaces[1];
1390 
1391  if (eFaces.size() > 2)
1392  {
1393  edgeFacesRest.insert
1394  (
1395  edgeI,
1396  SubList<label>(eFaces, eFaces.size()-2, 2)
1397  );
1398  }
1399  }
1400 
1401  edgeI++;
1402  }
1403 }
1404 
1405 
1406 void Foam::isoSurface::walkOrientation
1407 (
1408  const triSurface& surf,
1409  const List<FixedList<label, 3> >& faceEdges,
1410  const labelList& edgeFace0,
1411  const labelList& edgeFace1,
1412  const label seedTriI,
1413  labelList& flipState
1414 )
1415 {
1416  // Do walk for consistent orientation.
1417  DynamicList<label> changedFaces(surf.size());
1418 
1419  changedFaces.append(seedTriI);
1420 
1421  while (changedFaces.size())
1422  {
1423  DynamicList<label> newChangedFaces(changedFaces.size());
1424 
1425  forAll(changedFaces, i)
1426  {
1427  label triI = changedFaces[i];
1428  const labelledTri& tri = surf[triI];
1429  const FixedList<label, 3>& fEdges = faceEdges[triI];
1430 
1431  forAll(fEdges, fp)
1432  {
1433  label edgeI = fEdges[fp];
1434 
1435  // my points:
1436  label p0 = tri[fp];
1437  label p1 = tri[tri.fcIndex(fp)];
1438 
1439  label nbrI =
1440  (
1441  edgeFace0[edgeI] != triI
1442  ? edgeFace0[edgeI]
1443  : edgeFace1[edgeI]
1444  );
1445 
1446  if (nbrI != -1 && flipState[nbrI] == -1)
1447  {
1448  const labelledTri& nbrTri = surf[nbrI];
1449 
1450  // nbr points
1451  label nbrFp = findIndex(nbrTri, p0);
1452 
1453  if (nbrFp == -1)
1454  {
1455  FatalErrorIn("isoSurface::walkOrientation(..)")
1456  << "triI:" << triI
1457  << " tri:" << tri
1458  << " p0:" << p0
1459  << " p1:" << p1
1460  << " fEdges:" << fEdges
1461  << " edgeI:" << edgeI
1462  << " edgeFace0:" << edgeFace0[edgeI]
1463  << " edgeFace1:" << edgeFace1[edgeI]
1464  << " nbrI:" << nbrI
1465  << " nbrTri:" << nbrTri
1466  << abort(FatalError);
1467  }
1468 
1469 
1470  label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
1471 
1472  bool sameOrientation = (p1 == nbrP1);
1473 
1474  if (flipState[triI] == 0)
1475  {
1476  flipState[nbrI] = (sameOrientation ? 0 : 1);
1477  }
1478  else
1479  {
1480  flipState[nbrI] = (sameOrientation ? 1 : 0);
1481  }
1482  newChangedFaces.append(nbrI);
1483  }
1484  }
1485  }
1486 
1487  changedFaces.transfer(newChangedFaces);
1488  }
1489 }
1490 
1491 
1492 void Foam::isoSurface::orientSurface
1493 (
1494  triSurface& surf,
1495  const List<FixedList<label, 3> >& faceEdges,
1496  const labelList& edgeFace0,
1497  const labelList& edgeFace1,
1498  const Map<labelList>& edgeFacesRest
1499 )
1500 {
1501  // -1 : unvisited
1502  // 0 : leave as is
1503  // 1 : flip
1504  labelList flipState(surf.size(), -1);
1505 
1506  label seedTriI = 0;
1507 
1508  while (true)
1509  {
1510  // Find first unvisited triangle
1511  for
1512  (
1513  ;
1514  seedTriI < surf.size() && flipState[seedTriI] != -1;
1515  seedTriI++
1516  )
1517  {}
1518 
1519  if (seedTriI == surf.size())
1520  {
1521  break;
1522  }
1523 
1524  // Note: Determine orientation of seedTriI?
1525  // for now assume it is ok
1526  flipState[seedTriI] = 0;
1527 
1528  walkOrientation
1529  (
1530  surf,
1531  faceEdges,
1532  edgeFace0,
1533  edgeFace1,
1534  seedTriI,
1535  flipState
1536  );
1537  }
1538 
1539  // Do actual flipping
1540  surf.clearOut();
1541  forAll(surf, triI)
1542  {
1543  if (flipState[triI] == 1)
1544  {
1545  labelledTri tri(surf[triI]);
1546 
1547  surf[triI][0] = tri[0];
1548  surf[triI][1] = tri[2];
1549  surf[triI][2] = tri[1];
1550  }
1551  else if (flipState[triI] == -1)
1552  {
1553  FatalErrorIn
1554  (
1555  "isoSurface::orientSurface(triSurface&, const label)"
1556  ) << "problem" << abort(FatalError);
1557  }
1558  }
1559 }
1560 
1561 
1562 // Checks if triangle is connected through edgeI only.
1563 bool Foam::isoSurface::danglingTriangle
1564 (
1565  const FixedList<label, 3>& fEdges,
1566  const labelList& edgeFace1
1567 )
1568 {
1569  label nOpen = 0;
1570  forAll(fEdges, i)
1571  {
1572  if (edgeFace1[fEdges[i]] == -1)
1573  {
1574  nOpen++;
1575  }
1576  }
1577 
1578  if (nOpen == 1 || nOpen == 2 || nOpen == 3)
1579  {
1580  return true;
1581  }
1582  else
1583  {
1584  return false;
1585  }
1586 }
1587 
1588 
1589 // Mark triangles to keep. Returns number of dangling triangles.
1590 Foam::label Foam::isoSurface::markDanglingTriangles
1591 (
1592  const List<FixedList<label, 3> >& faceEdges,
1593  const labelList& edgeFace0,
1594  const labelList& edgeFace1,
1595  const Map<labelList>& edgeFacesRest,
1596  boolList& keepTriangles
1597 )
1598 {
1599  keepTriangles.setSize(faceEdges.size());
1600  keepTriangles = true;
1601 
1602  label nDangling = 0;
1603 
1604  // Remove any dangling triangles
1605  forAllConstIter(Map<labelList>, edgeFacesRest, iter)
1606  {
1607  // These are all the non-manifold edges. Filter out all triangles
1608  // with only one connected edge (= this edge)
1609 
1610  label edgeI = iter.key();
1611  const labelList& otherEdgeFaces = iter();
1612 
1613  // Remove all dangling triangles
1614  if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
1615  {
1616  keepTriangles[edgeFace0[edgeI]] = false;
1617  nDangling++;
1618  }
1619  if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
1620  {
1621  keepTriangles[edgeFace1[edgeI]] = false;
1622  nDangling++;
1623  }
1624  forAll(otherEdgeFaces, i)
1625  {
1626  label triI = otherEdgeFaces[i];
1627  if (danglingTriangle(faceEdges[triI], edgeFace1))
1628  {
1629  keepTriangles[triI] = false;
1630  nDangling++;
1631  }
1632  }
1633  }
1634  return nDangling;
1635 }
1636 
1637 
1638 Foam::triSurface Foam::isoSurface::subsetMesh
1639 (
1640  const triSurface& s,
1641  const labelList& newToOldFaces,
1642  labelList& oldToNewPoints,
1643  labelList& newToOldPoints
1644 )
1645 {
1646  const boolList include
1647  (
1648  createWithValues<boolList>
1649  (
1650  s.size(),
1651  false,
1652  newToOldFaces,
1653  true
1654  )
1655  );
1656 
1657  newToOldPoints.setSize(s.points().size());
1658  oldToNewPoints.setSize(s.points().size());
1659  oldToNewPoints = -1;
1660  {
1661  label pointI = 0;
1662 
1663  forAll(include, oldFacei)
1664  {
1665  if (include[oldFacei])
1666  {
1667  // Renumber labels for triangle
1668  const labelledTri& tri = s[oldFacei];
1669 
1670  forAll(tri, fp)
1671  {
1672  label oldPointI = tri[fp];
1673 
1674  if (oldToNewPoints[oldPointI] == -1)
1675  {
1676  oldToNewPoints[oldPointI] = pointI;
1677  newToOldPoints[pointI++] = oldPointI;
1678  }
1679  }
1680  }
1681  }
1682  newToOldPoints.setSize(pointI);
1683  }
1684 
1685  // Extract points
1686  pointField newPoints(newToOldPoints.size());
1687  forAll(newToOldPoints, i)
1688  {
1689  newPoints[i] = s.points()[newToOldPoints[i]];
1690  }
1691  // Extract faces
1692  List<labelledTri> newTriangles(newToOldFaces.size());
1693 
1694  forAll(newToOldFaces, i)
1695  {
1696  // Get old vertex labels
1697  const labelledTri& tri = s[newToOldFaces[i]];
1698 
1699  newTriangles[i][0] = oldToNewPoints[tri[0]];
1700  newTriangles[i][1] = oldToNewPoints[tri[1]];
1701  newTriangles[i][2] = oldToNewPoints[tri[2]];
1702  newTriangles[i].region() = tri.region();
1703  }
1704 
1705  // Reuse storage.
1706  return triSurface(newTriangles, s.patches(), newPoints, true);
1707 }
1708 
1709 
1710 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1711 
1714  const volScalarField& cVals,
1715  const scalarField& pVals,
1716  const scalar iso,
1717  const bool regularise,
1718  const scalar mergeTol
1719 )
1720 :
1721  mesh_(cVals.mesh()),
1722  pVals_(pVals),
1723  iso_(iso),
1724  regularise_(regularise),
1725  mergeDistance_(mergeTol*mesh_.bounds().mag())
1726 {
1727  if (debug)
1728  {
1729  Pout<< "isoSurface:" << nl
1730  << " isoField : " << cVals.name() << nl
1731  << " cell min/max : "
1732  << min(cVals.internalField()) << " / "
1733  << max(cVals.internalField()) << nl
1734  << " point min/max : "
1735  << min(pVals_) << " / "
1736  << max(pVals_) << nl
1737  << " isoValue : " << iso << nl
1738  << " regularise : " << regularise_ << nl
1739  << " mergeTol : " << mergeTol << nl
1740  << endl;
1741  }
1742 
1743  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1744 
1745 
1746  // Rewrite input field
1747  // ~~~~~~~~~~~~~~~~~~~
1748  // Rewrite input volScalarField to have interpolated values
1749  // on separated patches.
1750 
1751  cValsPtr_.reset(adaptPatchFields(cVals).ptr());
1752 
1753 
1754  // Construct cell centres field consistent with cVals
1755  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1756  // Generate field to interpolate. This is identical to the mesh.C()
1757  // except on separated coupled patches and on empty patches.
1758 
1759  slicedVolVectorField meshC
1760  (
1761  IOobject
1762  (
1763  "C",
1764  mesh_.pointsInstance(),
1765  mesh_.meshSubDir,
1766  mesh_,
1769  false
1770  ),
1771  mesh_,
1772  dimLength,
1773  mesh_.cellCentres(),
1774  mesh_.faceCentres()
1775  );
1776  forAll(patches, patchI)
1777  {
1778  const polyPatch& pp = patches[patchI];
1779 
1780  // Adapt separated coupled (proc and cyclic) patches
1781  if (pp.coupled())
1782  {
1783  fvPatchVectorField& pfld = const_cast<fvPatchVectorField&>
1784  (
1785  meshC.boundaryField()[patchI]
1786  );
1787 
1788  PackedBoolList isCollocated
1789  (
1790  collocatedFaces(refCast<const coupledPolyPatch>(pp))
1791  );
1792 
1793  forAll(isCollocated, i)
1794  {
1795  if (!isCollocated[i])
1796  {
1797  pfld[i] = mesh_.faceCentres()[pp.start()+i];
1798  }
1799  }
1800  }
1801  else if (isA<emptyPolyPatch>(pp))
1802  {
1804 
1805  bType& bfld = const_cast<bType&>(meshC.boundaryField());
1806 
1807  // Clear old value. Cannot resize it since is a slice.
1808  bfld.set(patchI, NULL);
1809 
1810  // Set new value we can change
1811  bfld.set
1812  (
1813  patchI,
1815  (
1816  mesh_.boundary()[patchI],
1817  meshC
1818  )
1819  );
1820 
1821  // Change to face centres
1822  bfld[patchI] = pp.patchSlice(mesh_.faceCentres());
1823  }
1824  }
1825 
1826 
1827 
1828  // Pre-calculate patch-per-face to avoid whichPatch call.
1829  labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces());
1830 
1831  forAll(patches, patchI)
1832  {
1833  const polyPatch& pp = patches[patchI];
1834 
1835  label faceI = pp.start();
1836 
1837  forAll(pp, i)
1838  {
1839  boundaryRegion[faceI-mesh_.nInternalFaces()] = patchI;
1840  faceI++;
1841  }
1842  }
1843 
1844 
1845 
1846  // Determine if any cut through face/cell
1847  calcCutTypes(boundaryRegion, meshC, cValsPtr_(), pVals_);
1848 
1849 
1850  DynamicList<point> snappedPoints(nCutCells_);
1851 
1852  // Per cc -1 or a point inside snappedPoints.
1853  labelList snappedCc;
1854  if (regularise_)
1855  {
1856  calcSnappedCc
1857  (
1858  boundaryRegion,
1859  meshC,
1860  cValsPtr_(),
1861  pVals_,
1862 
1863  snappedPoints,
1864  snappedCc
1865  );
1866  }
1867  else
1868  {
1869  snappedCc.setSize(mesh_.nCells());
1870  snappedCc = -1;
1871  }
1872 
1873 
1874 
1875  if (debug)
1876  {
1877  Pout<< "isoSurface : shifted " << snappedPoints.size()
1878  << " cell centres to intersection." << endl;
1879  }
1880 
1881  label nCellSnaps = snappedPoints.size();
1882 
1883 
1884  // Per point -1 or a point inside snappedPoints.
1885  labelList snappedPoint;
1886  if (regularise_)
1887  {
1888  // Determine if point is on boundary.
1889  PackedBoolList isBoundaryPoint(mesh_.nPoints());
1890 
1891  forAll(patches, patchI)
1892  {
1893  // Mark all boundary points that are not physically coupled
1894  // (so anything but collocated coupled patches)
1895 
1896  if (patches[patchI].coupled())
1897  {
1898  const coupledPolyPatch& cpp =
1899  refCast<const coupledPolyPatch>
1900  (
1901  patches[patchI]
1902  );
1903 
1904  PackedBoolList isCollocated(collocatedFaces(cpp));
1905 
1906  forAll(isCollocated, i)
1907  {
1908  if (!isCollocated[i])
1909  {
1910  const face& f = mesh_.faces()[cpp.start()+i];
1911 
1912  forAll(f, fp)
1913  {
1914  isBoundaryPoint.set(f[fp], 1);
1915  }
1916  }
1917  }
1918  }
1919  else
1920  {
1921  const polyPatch& pp = patches[patchI];
1922 
1923  forAll(pp, i)
1924  {
1925  const face& f = mesh_.faces()[pp.start()+i];
1926 
1927  forAll(f, fp)
1928  {
1929  isBoundaryPoint.set(f[fp], 1);
1930  }
1931  }
1932  }
1933  }
1934 
1935  calcSnappedPoint
1936  (
1937  isBoundaryPoint,
1938  boundaryRegion,
1939  meshC,
1940  cValsPtr_(),
1941  pVals_,
1942 
1943  snappedPoints,
1944  snappedPoint
1945  );
1946  }
1947  else
1948  {
1949  snappedPoint.setSize(mesh_.nPoints());
1950  snappedPoint = -1;
1951  }
1952 
1953  if (debug)
1954  {
1955  Pout<< "isoSurface : shifted " << snappedPoints.size()-nCellSnaps
1956  << " vertices to intersection." << endl;
1957  }
1958 
1959 
1960 
1961  DynamicList<point> triPoints(nCutCells_);
1962  DynamicList<label> triMeshCells(nCutCells_);
1963 
1964  generateTriPoints
1965  (
1966  cValsPtr_(),
1967  pVals_,
1968 
1969  meshC,
1970  mesh_.points(),
1971 
1972  snappedPoints,
1973  snappedCc,
1974  snappedPoint,
1975 
1976  triPoints,
1977  triMeshCells
1978  );
1979 
1980  if (debug)
1981  {
1982  Pout<< "isoSurface : generated " << triMeshCells.size()
1983  << " unmerged triangles from " << triPoints.size()
1984  << " unmerged points." << endl;
1985  }
1986 
1987 
1988  // Merge points and compact out non-valid triangles
1989  labelList triMap; // merged to unmerged triangle
1991  (
1992  stitchTriPoints
1993  (
1994  true, // check for duplicate tris
1995  triPoints,
1996  triPointMergeMap_, // unmerged to merged point
1997  triMap
1998  )
1999  );
2000 
2001  if (debug)
2002  {
2003  Pout<< "isoSurface : generated " << triMap.size()
2004  << " merged triangles." << endl;
2005  }
2006 
2007  meshCells_.setSize(triMap.size());
2008  forAll(triMap, i)
2009  {
2010  meshCells_[i] = triMeshCells[triMap[i]];
2011  }
2012 
2013  if (debug)
2014  {
2015  Pout<< "isoSurface : checking " << size()
2016  << " triangles for validity." << endl;
2017 
2018  forAll(*this, triI)
2019  {
2020  // Copied from surfaceCheck
2021  validTri(*this, triI);
2022  }
2023  }
2024 
2025 
2026  if (false)
2027  {
2028  List<FixedList<label, 3> > faceEdges;
2029  labelList edgeFace0, edgeFace1;
2030  Map<labelList> edgeFacesRest;
2031 
2032 
2033  while (true)
2034  {
2035  // Calculate addressing
2036  calcAddressing
2037  (
2038  *this,
2039  faceEdges,
2040  edgeFace0,
2041  edgeFace1,
2042  edgeFacesRest
2043  );
2044 
2045  // See if any dangling triangles
2046  boolList keepTriangles;
2047  label nDangling = markDanglingTriangles
2048  (
2049  faceEdges,
2050  edgeFace0,
2051  edgeFace1,
2052  edgeFacesRest,
2053  keepTriangles
2054  );
2055 
2056  if (debug)
2057  {
2058  Pout<< "isoSurface : detected " << nDangling
2059  << " dangling triangles." << endl;
2060  }
2061 
2062  if (nDangling == 0)
2063  {
2064  break;
2065  }
2066 
2067  // Create face map (new to old)
2068  labelList subsetTriMap(findIndices(keepTriangles, true));
2069 
2070  labelList subsetPointMap;
2071  labelList reversePointMap;
2073  (
2074  subsetMesh
2075  (
2076  *this,
2077  subsetTriMap,
2078  reversePointMap,
2079  subsetPointMap
2080  )
2081  );
2082  meshCells_ = labelField(meshCells_, subsetTriMap);
2083  inplaceRenumber(reversePointMap, triPointMergeMap_);
2084  }
2085 
2086  orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
2087  }
2088 
2089 
2090  if (debug)
2091  {
2092  fileName stlFile = mesh_.time().path() + ".stl";
2093  Pout<< "Dumping surface to " << stlFile << endl;
2094  triSurface::write(stlFile);
2095  }
2096 }
2097 
2098 
2099 // ************************************************************************* //
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:380
label size() const
Return the number of elements in the VectorSpace = nCmpt.
Definition: VectorSpaceI.H:67
const pointField & points
This boundary condition is not designed to be evaluated; it is assmued that the value is assigned via...
Foam::surfaceFields.
void invertManyToMany(const label len, const UList< InList > &, List< OutList > &)
Invert many-to-many.
bool set(const label) const
Is element set.
Definition: PtrListI.H:107
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:310
vector point
Point is a vector.
Definition: point.H:41
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject( name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE ))
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:390
dimensioned< scalar > mag(const dimensioned< Type > &)
Merge points. See below.
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
Specialization of GeometricField which holds slices of given complete fields in a form that they act ...
isoSurface(const volScalarField &cellIsoVals, const scalarField &pointIsoVals, const scalar iso, const bool regularise, const scalar mergeTol=1e-6)
Construct from cell values and point values. Uses boundaryField.
Definition: isoSurface.C:1713
A HashTable to objects of type <T> with a label key.
Definition: Map.H:49
Triangulated surface description with patch information.
Definition: triSurface.H:57
Field< label > labelField
Specialisation of Field<T> for label.
Definition: labelField.H:49
A bit-packed bool list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
InternalField & internalField()
Return internal field.
static const Vector max
Definition: Vector.H:82
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
friend Ostream & operator(Ostream &, const UList< T > &)
List< geometricSurfacePatch > geometricSurfacePatchList
const Mesh & mesh() const
Return mesh.
fileName path() const
Return directory path name (part before last /)
Definition: fileName.C:287
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Namespace for OpenFOAM.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
static const char nl
Definition: Ostream.H:260
const double e
Elementary charge.
Definition: doubleFloat.H:78
void setSize(const label)
Reset size of List.
Definition: List.C:318
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
#define WarningIn(functionName)
Report a warning using Foam::Warning.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
Foam::polyBoundaryMesh.
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
#define forAll(list, i)
Definition: UList.H:421
Macros for easy insertion into run-time selection tables.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:65
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const word & name() const
Return name.
Definition: IOobject.H:260
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
List< edge > edgeList
Definition: edgeList.H:38
const List< T >::subList patchSlice(const UList< T > &l) const
Slice list to patch.
Definition: polyPatch.H:333
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
label mergePoints(const UList< Type > &points, const scalar mergeTol, const bool verbose, labelList &pointMap, const Type &origin=Type::zero)
Sorts and merges points. All points closer than/equal mergeTol get merged.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
error FatalError
List< label > labelList
A List of labels.
Definition: labelList.H:56
A class for handling file names.
Definition: fileName.H:69
static const Vector zero
Definition: Vector.H:80
HashTable< T, label, Hash< label > >::iterator iterator
Definition: Map.H:56
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:139
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
label nPoints
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
friend class iterator
Declare friendship with the iterator.
Definition: HashTable.H:189
labelList findIndices(const ListType &, typename ListType::const_reference, const label start=0)
Find all occurences of given element. Linear search.
void inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
List< bool > boolList
Bool container classes.
Definition: boolList.H:50