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