triSurfaceTools.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2021 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 "triSurfaceTools.H"
27 
28 #include "triSurface.H"
29 #include "OFstream.H"
30 #include "mergePoints.H"
31 #include "polyMesh.H"
32 #include "plane.H"
33 #include "geompack.H"
34 #include "polygonTriangulate.H"
35 
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 
41 
42 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
43 
44 /*
45  Refine by splitting all three edges of triangle ('red' refinement).
46  Neighbouring triangles (which are not marked for refinement get split
47  in half ('green') refinement. (R. Verfuerth, "A review of a posteriori
48  error estimation and adaptive mesh refinement techniques",
49  Wiley-Teubner, 1996)
50 */
51 
52 // Facei gets refined ('red'). Propagate edge refinement.
53 void Foam::triSurfaceTools::calcRefineStatus
54 (
55  const triSurface& surf,
56  const label facei,
57  List<refineType>& refine
58 )
59 {
60  if (refine[facei] == RED)
61  {
62  // Already marked for refinement. Do nothing.
63  }
64  else
65  {
66  // Not marked or marked for 'green' refinement. Refine.
67  refine[facei] = RED;
68 
69  const labelList& myNeighbours = surf.faceFaces()[facei];
70 
71  forAll(myNeighbours, myNeighbourI)
72  {
73  label neighbourFacei = myNeighbours[myNeighbourI];
74 
75  if (refine[neighbourFacei] == GREEN)
76  {
77  // Change to red refinement and propagate
78  calcRefineStatus(surf, neighbourFacei, refine);
79  }
80  else if (refine[neighbourFacei] == NONE)
81  {
82  refine[neighbourFacei] = GREEN;
83  }
84  }
85  }
86 }
87 
88 
89 // Split facei along edgeI at position newPointi
90 void Foam::triSurfaceTools::greenRefine
91 (
92  const triSurface& surf,
93  const label facei,
94  const label edgeI,
95  const label newPointi,
96  DynamicList<labelledTri>& newFaces
97 )
98 {
99  const labelledTri& f = surf.localFaces()[facei];
100  const edge& e = surf.edges()[edgeI];
101 
102  // Find index of edge in face.
103 
104  label fp0 = findIndex(f, e[0]);
105  label fp1 = f.fcIndex(fp0);
106  label fp2 = f.fcIndex(fp1);
107 
108  if (f[fp1] == e[1])
109  {
110  // Edge oriented like face
111  newFaces.append
112  (
113  labelledTri
114  (
115  f[fp0],
116  newPointi,
117  f[fp2],
118  f.region()
119  )
120  );
121  newFaces.append
122  (
123  labelledTri
124  (
125  newPointi,
126  f[fp1],
127  f[fp2],
128  f.region()
129  )
130  );
131  }
132  else
133  {
134  newFaces.append
135  (
136  labelledTri
137  (
138  f[fp2],
139  newPointi,
140  f[fp1],
141  f.region()
142  )
143  );
144  newFaces.append
145  (
146  labelledTri
147  (
148  newPointi,
149  f[fp0],
150  f[fp1],
151  f.region()
152  )
153  );
154  }
155 }
156 
157 
158 // Refine all triangles marked for refinement.
159 Foam::triSurface Foam::triSurfaceTools::doRefine
160 (
161  const triSurface& surf,
162  const List<refineType>& refineStatus
163 )
164 {
165  // Storage for new points. (start after old points)
166  DynamicList<point> newPoints(surf.nPoints());
167  forAll(surf.localPoints(), pointi)
168  {
169  newPoints.append(surf.localPoints()[pointi]);
170  }
171  label newVertI = surf.nPoints();
172 
173  // Storage for new faces
174  DynamicList<labelledTri> newFaces(surf.size());
175 
176 
177  // Point index for midpoint on edge
178  labelList edgeMid(surf.nEdges(), -1);
179 
180  forAll(refineStatus, facei)
181  {
182  if (refineStatus[facei] == RED)
183  {
184  // Create new vertices on all edges to be refined.
185  const labelList& fEdges = surf.faceEdges()[facei];
186 
187  forAll(fEdges, i)
188  {
189  label edgeI = fEdges[i];
190 
191  if (edgeMid[edgeI] == -1)
192  {
193  const edge& e = surf.edges()[edgeI];
194 
195  // Create new point on mid of edge
196  newPoints.append
197  (
198  0.5
199  * (
200  surf.localPoints()[e.start()]
201  + surf.localPoints()[e.end()]
202  )
203  );
204  edgeMid[edgeI] = newVertI++;
205  }
206  }
207 
208  // Now we have new mid edge vertices for all edges on face
209  // so create triangles for RED rerfinement.
210 
211  const edgeList& edges = surf.edges();
212 
213  // Corner triangles
214  newFaces.append
215  (
216  labelledTri
217  (
218  edgeMid[fEdges[0]],
219  edges[fEdges[0]].commonVertex(edges[fEdges[1]]),
220  edgeMid[fEdges[1]],
221  surf[facei].region()
222  )
223  );
224 
225  newFaces.append
226  (
227  labelledTri
228  (
229  edgeMid[fEdges[1]],
230  edges[fEdges[1]].commonVertex(edges[fEdges[2]]),
231  edgeMid[fEdges[2]],
232  surf[facei].region()
233  )
234  );
235 
236  newFaces.append
237  (
238  labelledTri
239  (
240  edgeMid[fEdges[2]],
241  edges[fEdges[2]].commonVertex(edges[fEdges[0]]),
242  edgeMid[fEdges[0]],
243  surf[facei].region()
244  )
245  );
246 
247  // Inner triangle
248  newFaces.append
249  (
250  labelledTri
251  (
252  edgeMid[fEdges[0]],
253  edgeMid[fEdges[1]],
254  edgeMid[fEdges[2]],
255  surf[facei].region()
256  )
257  );
258 
259 
260  // Create triangles for GREEN refinement.
261  forAll(fEdges, i)
262  {
263  const label edgeI = fEdges[i];
264 
265  label otherFacei = otherFace(surf, facei, edgeI);
266 
267  if ((otherFacei != -1) && (refineStatus[otherFacei] == GREEN))
268  {
269  greenRefine
270  (
271  surf,
272  otherFacei,
273  edgeI,
274  edgeMid[edgeI],
275  newFaces
276  );
277  }
278  }
279  }
280  }
281 
282  // Copy unmarked triangles since keep original vertex numbering.
283  forAll(refineStatus, facei)
284  {
285  if (refineStatus[facei] == NONE)
286  {
287  newFaces.append(surf.localFaces()[facei]);
288  }
289  }
290 
291  newFaces.shrink();
292  newPoints.shrink();
293 
294 
295  // Transfer DynamicLists to straight ones.
297  allPoints.transfer(newPoints);
298  newPoints.clear();
299 
300  return triSurface(newFaces, surf.patches(), allPoints, true);
301 }
302 
303 
304 // Angle between two neighbouring triangles,
305 // angle between shared-edge end points and left and right face end points
306 Foam::scalar Foam::triSurfaceTools::faceCosAngle
307 (
308  const point& pStart,
309  const point& pEnd,
310  const point& pLeft,
311  const point& pRight
312 )
313 {
314  const vector common(pEnd - pStart);
315  const vector base0(pLeft - pStart);
316  const vector base1(pRight - pStart);
317 
318  vector n0(common ^ base0);
319  n0 /= Foam::mag(n0);
320 
321  vector n1(base1 ^ common);
322  n1 /= Foam::mag(n1);
323 
324  return n0 & n1;
325 }
326 
327 
328 // Protect edges around vertex from collapsing.
329 // Moving vertI to a different position
330 // - affects obviously the cost of the faces using it
331 // - but also their neighbours since the edge between the faces changes
332 void Foam::triSurfaceTools::protectNeighbours
333 (
334  const triSurface& surf,
335  const label vertI,
336  labelList& faceStatus
337 )
338 {
339 // const labelList& myFaces = surf.pointFaces()[vertI];
340 // forAll(myFaces, i)
341 // {
342 // label facei = myFaces[i];
343 //
344 // if ((faceStatus[facei] == ANYEDGE) || (faceStatus[facei] >= 0))
345 // {
346 // faceStatus[facei] = NOEDGE;
347 // }
348 // }
349 
350  const labelList& myEdges = surf.pointEdges()[vertI];
351  forAll(myEdges, i)
352  {
353  const labelList& myFaces = surf.edgeFaces()[myEdges[i]];
354 
355  forAll(myFaces, myFacei)
356  {
357  label facei = myFaces[myFacei];
358 
359  if ((faceStatus[facei] == ANYEDGE) || (faceStatus[facei] >= 0))
360  {
361  faceStatus[facei] = NOEDGE;
362  }
363  }
364  }
365 }
366 
367 
368 //
369 // Edge collapse helper functions
370 //
371 
372 // Get all faces that will get collapsed if edgeI collapses.
373 Foam::labelHashSet Foam::triSurfaceTools::getCollapsedFaces
374 (
375  const triSurface& surf,
376  label edgeI
377 )
378 {
379  const edge& e = surf.edges()[edgeI];
380  label v1 = e.start();
381  label v2 = e.end();
382 
383  // Faces using edge will certainly get collapsed.
384  const labelList& myFaces = surf.edgeFaces()[edgeI];
385 
386  labelHashSet facesToBeCollapsed(2*myFaces.size());
387 
388  forAll(myFaces, myFacei)
389  {
390  facesToBeCollapsed.insert(myFaces[myFacei]);
391  }
392 
393  // From faces using v1 check if they share an edge with faces
394  // using v2.
395  // - share edge: are part of 'splay' tree and will collapse if edge
396  // collapses
397  const labelList& v1Faces = surf.pointFaces()[v1];
398 
399  forAll(v1Faces, v1Facei)
400  {
401  label face1I = v1Faces[v1Facei];
402 
403  label otherEdgeI = oppositeEdge(surf, face1I, v1);
404 
405  // Step across edge to other face
406  label face2I = otherFace(surf, face1I, otherEdgeI);
407 
408  if (face2I != -1)
409  {
410  // found face on other side of edge. Now check if uses v2.
411  if (oppositeVertex(surf, face2I, otherEdgeI) == v2)
412  {
413  // triangles face1I and face2I form splay tree and will
414  // collapse.
415  facesToBeCollapsed.insert(face1I);
416  facesToBeCollapsed.insert(face2I);
417  }
418  }
419  }
420 
421  return facesToBeCollapsed;
422 }
423 
424 
425 // Return value of faceUsed for faces using vertI
426 Foam::label Foam::triSurfaceTools::vertexUsesFace
427 (
428  const triSurface& surf,
429  const labelHashSet& faceUsed,
430  const label vertI
431 )
432 {
433  const labelList& myFaces = surf.pointFaces()[vertI];
434 
435  forAll(myFaces, myFacei)
436  {
437  label face1I = myFaces[myFacei];
438 
439  if (faceUsed.found(face1I))
440  {
441  return face1I;
442  }
443  }
444  return -1;
445 }
446 
447 
448 // Calculate new edge-face addressing resulting from edge collapse
449 void Foam::triSurfaceTools::getMergedEdges
450 (
451  const triSurface& surf,
452  const label edgeI,
453  const labelHashSet& collapsedFaces,
454  HashTable<label, label, Hash<label>>& edgeToEdge,
455  HashTable<label, label, Hash<label>>& edgeToFace
456 )
457 {
458  const edge& e = surf.edges()[edgeI];
459  label v1 = e.start();
460  label v2 = e.end();
461 
462  const labelList& v1Faces = surf.pointFaces()[v1];
463  const labelList& v2Faces = surf.pointFaces()[v2];
464 
465  // Mark all (non collapsed) faces using v2
466  labelHashSet v2FacesHash(v2Faces.size());
467 
468  forAll(v2Faces, v2Facei)
469  {
470  if (!collapsedFaces.found(v2Faces[v2Facei]))
471  {
472  v2FacesHash.insert(v2Faces[v2Facei]);
473  }
474  }
475 
476 
477  forAll(v1Faces, v1Facei)
478  {
479  label face1I = v1Faces[v1Facei];
480 
481  if (collapsedFaces.found(face1I))
482  {
483  continue;
484  }
485 
486  //
487  // Check if face1I uses a vertex connected to a face using v2
488  //
489 
490  label vert1I = -1;
491  label vert2I = -1;
493  (
494  surf,
495  face1I,
496  v1,
497  vert1I,
498  vert2I
499  );
500  // Pout<< "Face:" << surf.localFaces()[face1I] << " other vertices:"
501  // << vert1I << ' ' << vert2I << endl;
502 
503  // Check vert1, vert2 for usage by v2Face.
504 
505  label commonVert = vert1I;
506  label face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
507  if (face2I == -1)
508  {
509  commonVert = vert2I;
510  face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
511  }
512 
513  if (face2I != -1)
514  {
515  // Found one: commonVert is used by both face1I and face2I
516  label edge1I = getEdge(surf, v1, commonVert);
517  label edge2I = getEdge(surf, v2, commonVert);
518 
519  edgeToEdge.insert(edge1I, edge2I);
520  edgeToEdge.insert(edge2I, edge1I);
521 
522  edgeToFace.insert(edge1I, face2I);
523  edgeToFace.insert(edge2I, face1I);
524  }
525  }
526 }
527 
528 
529 // Calculates (cos of) angle across edgeI of facei,
530 // taking into account updated addressing (resulting from edge collapse)
531 Foam::scalar Foam::triSurfaceTools::edgeCosAngle
532 (
533  const triSurface& surf,
534  const label v1,
535  const point& pt,
536  const labelHashSet& collapsedFaces,
537  const HashTable<label, label, Hash<label>>& edgeToEdge,
538  const HashTable<label, label, Hash<label>>& edgeToFace,
539  const label facei,
540  const label edgeI
541 )
542 {
543  const pointField& localPoints = surf.localPoints();
544 
545  label A = surf.edges()[edgeI].start();
546  label B = surf.edges()[edgeI].end();
547  label C = oppositeVertex(surf, facei, edgeI);
548 
549  label D = -1;
550 
551  label face2I = -1;
552 
553  if (edgeToEdge.found(edgeI))
554  {
555  // Use merged addressing
556  label edge2I = edgeToEdge[edgeI];
557  face2I = edgeToFace[edgeI];
558 
559  D = oppositeVertex(surf, face2I, edge2I);
560  }
561  else
562  {
563  // Use normal edge-face addressing
564  face2I = otherFace(surf, facei, edgeI);
565 
566  if ((face2I != -1) && !collapsedFaces.found(face2I))
567  {
568  D = oppositeVertex(surf, face2I, edgeI);
569  }
570  }
571 
572  scalar cosAngle = 1;
573 
574  if (D != -1)
575  {
576  if (A == v1)
577  {
578  cosAngle = faceCosAngle
579  (
580  pt,
581  localPoints[B],
582  localPoints[C],
583  localPoints[D]
584  );
585  }
586  else if (B == v1)
587  {
588  cosAngle = faceCosAngle
589  (
590  localPoints[A],
591  pt,
592  localPoints[C],
593  localPoints[D]
594  );
595  }
596  else if (C == v1)
597  {
598  cosAngle = faceCosAngle
599  (
600  localPoints[A],
601  localPoints[B],
602  pt,
603  localPoints[D]
604  );
605  }
606  else if (D == v1)
607  {
608  cosAngle = faceCosAngle
609  (
610  localPoints[A],
611  localPoints[B],
612  localPoints[C],
613  pt
614  );
615  }
616  else
617  {
619  << "face " << facei << " does not use vertex "
620  << v1 << " of collapsed edge" << abort(FatalError);
621  }
622  }
623  return cosAngle;
624 }
625 
626 
627 Foam::scalar Foam::triSurfaceTools::collapseMinCosAngle
628 (
629  const triSurface& surf,
630  const label v1,
631  const point& pt,
632  const labelHashSet& collapsedFaces,
633  const HashTable<label, label, Hash<label>>& edgeToEdge,
634  const HashTable<label, label, Hash<label>>& edgeToFace
635 )
636 {
637  const labelList& v1Faces = surf.pointFaces()[v1];
638 
639  scalar minCos = 1;
640 
641  forAll(v1Faces, v1Facei)
642  {
643  label facei = v1Faces[v1Facei];
644 
645  if (collapsedFaces.found(facei))
646  {
647  continue;
648  }
649 
650  const labelList& myEdges = surf.faceEdges()[facei];
651 
652  forAll(myEdges, myEdgeI)
653  {
654  label edgeI = myEdges[myEdgeI];
655 
656  minCos =
657  min
658  (
659  minCos,
660  edgeCosAngle
661  (
662  surf,
663  v1,
664  pt,
665  collapsedFaces,
666  edgeToEdge,
667  edgeToFace,
668  facei,
669  edgeI
670  )
671  );
672  }
673  }
674 
675  return minCos;
676 }
677 
678 
679 // Calculate max dihedral angle after collapsing edge to v1 (at pt).
680 // Note that all edges of all faces using v1 are affected.
681 bool Foam::triSurfaceTools::collapseCreatesFold
682 (
683  const triSurface& surf,
684  const label v1,
685  const point& pt,
686  const labelHashSet& collapsedFaces,
687  const HashTable<label, label, Hash<label>>& edgeToEdge,
688  const HashTable<label, label, Hash<label>>& edgeToFace,
689  const scalar minCos
690 )
691 {
692  const labelList& v1Faces = surf.pointFaces()[v1];
693 
694  forAll(v1Faces, v1Facei)
695  {
696  label facei = v1Faces[v1Facei];
697 
698  if (collapsedFaces.found(facei))
699  {
700  continue;
701  }
702 
703  const labelList& myEdges = surf.faceEdges()[facei];
704 
705  forAll(myEdges, myEdgeI)
706  {
707  label edgeI = myEdges[myEdgeI];
708 
709  if
710  (
711  edgeCosAngle
712  (
713  surf,
714  v1,
715  pt,
716  collapsedFaces,
717  edgeToEdge,
718  edgeToFace,
719  facei,
720  edgeI
721  )
722  < minCos
723  )
724  {
725  return true;
726  }
727  }
728  }
729 
730  return false;
731 }
732 
733 
736 // a vertex.
737 //bool Foam::triSurfaceTools::collapseCreatesDuplicates
738 //(
739 // const triSurface& surf,
740 // const label edgeI,
741 // const labelHashSet& collapsedFaces
742 //)
743 //{
744 //Pout<< "duplicate : edgeI:" << edgeI << surf.edges()[edgeI]
745 // << " collapsedFaces:" << collapsedFaces.toc() << endl;
746 //
747 // // Mark neighbours of faces to be collapsed, i.e. get the first layer
748 // // of triangles outside collapsedFaces.
749 // // neighbours actually contains the
750 // // edge with which triangle connects to collapsedFaces.
751 //
752 // HashTable<label, label, Hash<label>> neighbours;
753 //
754 // labelList collapsed = collapsedFaces.toc();
755 //
756 // forAll(collapsed, collapseI)
757 // {
758 // const label facei = collapsed[collapseI];
759 //
760 // const labelList& myEdges = surf.faceEdges()[facei];
761 //
762 // Pout<< "collapsing facei:" << facei << " uses edges:" << myEdges
763 // << endl;
764 //
765 // forAll(myEdges, myEdgeI)
766 // {
767 // const labelList& myFaces = surf.edgeFaces()[myEdges[myEdgeI]];
768 //
769 // Pout<< "Edge " << myEdges[myEdgeI] << " is used by faces "
770 // << myFaces << endl;
771 //
772 // if ((myEdges[myEdgeI] != edgeI) && (myFaces.size() == 2))
773 // {
774 // // Get the other face
775 // label neighbourFacei = myFaces[0];
776 //
777 // if (neighbourFacei == facei)
778 // {
779 // neighbourFacei = myFaces[1];
780 // }
781 //
782 // // Is 'outside' face if not in collapsedFaces itself
783 // if (!collapsedFaces.found(neighbourFacei))
784 // {
785 // neighbours.insert(neighbourFacei, myEdges[myEdgeI]);
786 // }
787 // }
788 // }
789 // }
790 //
791 // // Now neighbours contains first layer of triangles outside of
792 // // collapseFaces
793 // // There should be
794 // // -two if edgeI is a boundary edge
795 // // since the outside 'edge' of collapseFaces should
796 // // form a triangle and the face connected to edgeI is not inserted.
797 // // -four if edgeI is not a boundary edge since then the outside edge forms
798 // // a diamond.
799 //
800 // // Check if any of the faces in neighbours share an edge. (n^2)
801 //
802 // labelList neighbourList = neighbours.toc();
803 //
804 //Pout<< "edgeI:" << edgeI << " neighbourList:" << neighbourList << endl;
805 //
806 //
807 // forAll(neighbourList, i)
808 // {
809 // const labelList& faceIEdges = surf.faceEdges()[neighbourList[i]];
810 //
811 // for (label j = i+1; j < neighbourList.size(); i++)
812 // {
813 // const labelList& faceJEdges = surf.faceEdges()[neighbourList[j]];
814 //
815 // // Check if facei and faceJ share an edge
816 // forAll(faceIEdges, fI)
817 // {
818 // forAll(faceJEdges, fJ)
819 // {
820 // Pout<< " comparing " << faceIEdges[fI] << " to "
821 // << faceJEdges[fJ] << endl;
822 //
823 // if (faceIEdges[fI] == faceJEdges[fJ])
824 // {
825 // return true;
826 // }
827 // }
828 // }
829 // }
830 // }
831 // Pout<< "Found no match. Returning false" << endl;
832 // return false;
833 //}
834 
835 
836 // Finds the triangle edge cut by the plane between a point inside / on edge
837 // of a triangle and a point outside. Returns:
838 // - cut through edge/point
839 // - miss()
840 Foam::surfaceLocation Foam::triSurfaceTools::cutEdge
841 (
842  const triSurface& s,
843  const label triI,
844  const label excludeEdgeI,
845  const label excludePointi,
846 
847  const point& triPoint,
848  const plane& cutPlane,
849  const point& toPoint
850 )
851 {
852  const pointField& points = s.points();
853  const labelledTri& f = s[triI];
854  const labelList& fEdges = s.faceEdges()[triI];
855 
856  // Get normal distance to planeN
857  FixedList<scalar, 3> d;
858 
859  scalar norm = 0;
860  forAll(d, fp)
861  {
862  d[fp] = (points[f[fp]]-cutPlane.refPoint()) & cutPlane.normal();
863  norm += mag(d[fp]);
864  }
865 
866  // Normalise and truncate
867  forAll(d, i)
868  {
869  d[i] /= norm;
870 
871  if (mag(d[i]) < 1e-6)
872  {
873  d[i] = 0.0;
874  }
875  }
876 
877  // Return information
878  surfaceLocation cut;
879 
880  if (excludePointi != -1)
881  {
882  // Excluded point. Test only opposite edge.
883 
884  label fp0 = findIndex(s.localFaces()[triI], excludePointi);
885 
886  if (fp0 == -1)
887  {
889  << " localF:" << s.localFaces()[triI] << abort(FatalError);
890  }
891 
892  label fp1 = f.fcIndex(fp0);
893  label fp2 = f.fcIndex(fp1);
894 
895 
896  if (d[fp1] == 0.0)
897  {
898  cut.setHit();
899  cut.setPoint(points[f[fp1]]);
900  cut.elementType() = triPointRef::POINT;
901  cut.setIndex(s.localFaces()[triI][fp1]);
902  }
903  else if (d[fp2] == 0.0)
904  {
905  cut.setHit();
906  cut.setPoint(points[f[fp2]]);
907  cut.elementType() = triPointRef::POINT;
908  cut.setIndex(s.localFaces()[triI][fp2]);
909  }
910  else if
911  (
912  (d[fp1] < 0 && d[fp2] < 0)
913  || (d[fp1] > 0 && d[fp2] > 0)
914  )
915  {
916  // Both same sign. Not crossing edge at all.
917  // cut already set to miss().
918  }
919  else
920  {
921  cut.setHit();
922  cut.setPoint
923  (
924  (d[fp2]*points[f[fp1]] - d[fp1]*points[f[fp2]])
925  / (d[fp2] - d[fp1])
926  );
927  cut.elementType() = triPointRef::EDGE;
928  cut.setIndex(fEdges[fp1]);
929  }
930  }
931  else
932  {
933  // Find the two intersections
934  FixedList<surfaceLocation, 2> inters;
935  label interI = 0;
936 
937  forAll(f, fp0)
938  {
939  label fp1 = f.fcIndex(fp0);
940 
941  if (d[fp0] == 0)
942  {
943  if (interI >= 2)
944  {
946  << "problem : triangle has three intersections." << nl
947  << "triangle:" << f.tri(points)
948  << " d:" << d << abort(FatalError);
949  }
950  inters[interI].setHit();
951  inters[interI].setPoint(points[f[fp0]]);
952  inters[interI].elementType() = triPointRef::POINT;
953  inters[interI].setIndex(s.localFaces()[triI][fp0]);
954  interI++;
955  }
956  else if
957  (
958  (d[fp0] < 0 && d[fp1] > 0)
959  || (d[fp0] > 0 && d[fp1] < 0)
960  )
961  {
962  if (interI >= 2)
963  {
965  << "problem : triangle has three intersections." << nl
966  << "triangle:" << f.tri(points)
967  << " d:" << d << abort(FatalError);
968  }
969  inters[interI].setHit();
970  inters[interI].setPoint
971  (
972  (d[fp0]*points[f[fp1]] - d[fp1]*points[f[fp0]])
973  / (d[fp0] - d[fp1])
974  );
975  inters[interI].elementType() = triPointRef::EDGE;
976  inters[interI].setIndex(fEdges[fp0]);
977  interI++;
978  }
979  }
980 
981 
982  if (interI == 0)
983  {
984  // Return miss
985  }
986  else if (interI == 1)
987  {
988  // Only one intersection. Should not happen!
989  cut = inters[0];
990  }
991  else if (interI == 2)
992  {
993  // Handle excludeEdgeI
994  if
995  (
996  inters[0].elementType() == triPointRef::EDGE
997  && inters[0].index() == excludeEdgeI
998  )
999  {
1000  cut = inters[1];
1001  }
1002  else if
1003  (
1004  inters[1].elementType() == triPointRef::EDGE
1005  && inters[1].index() == excludeEdgeI
1006  )
1007  {
1008  cut = inters[0];
1009  }
1010  else
1011  {
1012  // Two cuts. Find nearest.
1013  if
1014  (
1015  magSqr(inters[0].rawPoint() - toPoint)
1016  < magSqr(inters[1].rawPoint() - toPoint)
1017  )
1018  {
1019  cut = inters[0];
1020  }
1021  else
1022  {
1023  cut = inters[1];
1024  }
1025  }
1026  }
1027  }
1028  return cut;
1029 }
1030 
1031 
1032 // 'Snap' point on to endPoint.
1033 void Foam::triSurfaceTools::snapToEnd
1034 (
1035  const triSurface& s,
1036  const surfaceLocation& end,
1037  surfaceLocation& current
1038 )
1039 {
1040  if (end.elementType() == triPointRef::NONE)
1041  {
1042  if (current.elementType() == triPointRef::NONE)
1043  {
1044  // endpoint on triangle; current on triangle
1045  if (current.index() == end.index())
1046  {
1047  // if (debug)
1048  //{
1049  // Pout<< "snapToEnd : snapping:" << current << " onto:"
1050  // << end << endl;
1051  //}
1052  current = end;
1053  current.setHit();
1054  }
1055  }
1056  // No need to handle current on edge/point since tracking handles this.
1057  }
1058  else if (end.elementType() == triPointRef::EDGE)
1059  {
1060  if (current.elementType() == triPointRef::NONE)
1061  {
1062  // endpoint on edge; current on triangle
1063  const labelList& fEdges = s.faceEdges()[current.index()];
1064 
1065  if (findIndex(fEdges, end.index()) != -1)
1066  {
1067  // if (debug)
1068  //{
1069  // Pout<< "snapToEnd : snapping:" << current << " onto:"
1070  // << end << endl;
1071  //}
1072  current = end;
1073  current.setHit();
1074  }
1075  }
1076  else if (current.elementType() == triPointRef::EDGE)
1077  {
1078  // endpoint on edge; current on edge
1079  if (current.index() == end.index())
1080  {
1081  // if (debug)
1082  //{
1083  // Pout<< "snapToEnd : snapping:" << current << " onto:"
1084  // << end << endl;
1085  //}
1086  current = end;
1087  current.setHit();
1088  }
1089  }
1090  else
1091  {
1092  // endpoint on edge; current on point
1093  const edge& e = s.edges()[end.index()];
1094 
1095  if (current.index() == e[0] || current.index() == e[1])
1096  {
1097  // if (debug)
1098  //{
1099  // Pout<< "snapToEnd : snapping:" << current << " onto:"
1100  // << end << endl;
1101  //}
1102  current = end;
1103  current.setHit();
1104  }
1105  }
1106  }
1107  else // end.elementType() == POINT
1108  {
1109  if (current.elementType() == triPointRef::NONE)
1110  {
1111  // endpoint on point; current on triangle
1112  const triSurface::FaceType& f = s.localFaces()[current.index()];
1113 
1114  if (findIndex(f, end.index()) != -1)
1115  {
1116  // if (debug)
1117  //{
1118  // Pout<< "snapToEnd : snapping:" << current << " onto:"
1119  // << end << endl;
1120  //}
1121  current = end;
1122  current.setHit();
1123  }
1124  }
1125  else if (current.elementType() == triPointRef::EDGE)
1126  {
1127  // endpoint on point; current on edge
1128  const edge& e = s.edges()[current.index()];
1129 
1130  if (end.index() == e[0] || end.index() == e[1])
1131  {
1132  // if (debug)
1133  //{
1134  // Pout<< "snapToEnd : snapping:" << current << " onto:"
1135  // << end << endl;
1136  //}
1137  current = end;
1138  current.setHit();
1139  }
1140  }
1141  else
1142  {
1143  // endpoint on point; current on point
1144  if (current.index() == end.index())
1145  {
1146  // if (debug)
1147  //{
1148  // Pout<< "snapToEnd : snapping:" << current << " onto:"
1149  // << end << endl;
1150  //}
1151  current = end;
1152  current.setHit();
1153  }
1154  }
1155  }
1156 }
1157 
1158 
1159 // Start:
1160 // - location
1161 // - element type (triangle/edge/point) and index
1162 // - triangle to exclude
1163 Foam::surfaceLocation Foam::triSurfaceTools::visitFaces
1164 (
1165  const triSurface& s,
1166  const labelList& eFaces,
1167  const surfaceLocation& start,
1168  const label excludeEdgeI,
1169  const label excludePointi,
1170  const surfaceLocation& end,
1171  const plane& cutPlane
1172 )
1173 {
1174  surfaceLocation nearest;
1175 
1176  scalar minDistSqr = Foam::sqr(great);
1177 
1178  forAll(eFaces, i)
1179  {
1180  label triI = eFaces[i];
1181 
1182  // Make sure we don't revisit previous face
1183  if (triI != start.triangle())
1184  {
1185  if (end.elementType() == triPointRef::NONE && end.index() == triI)
1186  {
1187  // Endpoint is in this triangle. Jump there.
1188  nearest = end;
1189  nearest.setHit();
1190  nearest.triangle() = triI;
1191  break;
1192  }
1193  else
1194  {
1195  // Which edge is cut.
1196 
1197  surfaceLocation cutInfo = cutEdge
1198  (
1199  s,
1200  triI,
1201  excludeEdgeI, // excludeEdgeI
1202  excludePointi, // excludePointi
1203  start.rawPoint(),
1204  cutPlane,
1205  end.rawPoint()
1206  );
1207 
1208  // If crossing an edge we expect next edge to be cut.
1209  if (excludeEdgeI != -1 && !cutInfo.hit())
1210  {
1212  << "Triangle:" << triI
1213  << " excludeEdge:" << excludeEdgeI
1214  << " point:" << start.rawPoint()
1215  << " plane:" << cutPlane
1216  << " . No intersection!" << abort(FatalError);
1217  }
1218 
1219  if (cutInfo.hit())
1220  {
1221  scalar distSqr = magSqr(cutInfo.rawPoint()-end.rawPoint());
1222 
1223  if (distSqr < minDistSqr)
1224  {
1225  minDistSqr = distSqr;
1226  nearest = cutInfo;
1227  nearest.triangle() = triI;
1228  nearest.setMiss();
1229  }
1230  }
1231  }
1232  }
1233  }
1234 
1235  if (nearest.triangle() == -1)
1236  {
1237  // Did not move from edge. Give warning? Return something special?
1238  // For now responsibility of caller to make sure that nothing has
1239  // moved.
1240  }
1241 
1242  return nearest;
1243 }
1244 
1245 
1246 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1247 
1248 
1249 // Write pointField to file
1252  const fileName& fName,
1253  const pointField& pts
1254 )
1255 {
1256  OFstream outFile(fName);
1257 
1258  forAll(pts, pointi)
1259  {
1260  const point& pt = pts[pointi];
1261 
1262  outFile<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
1263  }
1264  Pout<< "Written " << pts.size() << " vertices to file " << fName << endl;
1265 }
1266 
1267 
1268 // Write vertex subset to OBJ format file
1271  const triSurface& surf,
1272  const fileName& fName,
1273  const boolList& markedVerts
1274 )
1275 {
1276  OFstream outFile(fName);
1277 
1278  label nVerts = 0;
1279  forAll(markedVerts, vertI)
1280  {
1281  if (markedVerts[vertI])
1282  {
1283  const point& pt = surf.localPoints()[vertI];
1284 
1285  outFile<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
1286 
1287  nVerts++;
1288  }
1289  }
1290  Pout<< "Written " << nVerts << " vertices to file " << fName << endl;
1291 }
1292 
1293 
1294 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1295 // Addressing helper functions:
1296 
1297 
1298 // Get all triangles using vertices of edge
1301  const triSurface& surf,
1302  const label edgeI,
1303  labelList& edgeTris
1304 )
1305 {
1306  const edge& e = surf.edges()[edgeI];
1307  const labelList& myFaces = surf.edgeFaces()[edgeI];
1308 
1309  label face1I = myFaces[0];
1310  label face2I = -1;
1311  if (myFaces.size() == 2)
1312  {
1313  face2I = myFaces[1];
1314  }
1315 
1316  const labelList& startFaces = surf.pointFaces()[e.start()];
1317  const labelList& endFaces = surf.pointFaces()[e.end()];
1318 
1319  // Number of triangles is sum of pointfaces - common faces
1320  // (= faces using edge)
1321  edgeTris.setSize(startFaces.size() + endFaces.size() - myFaces.size());
1322 
1323  label nTris = 0;
1324  forAll(startFaces, startFacei)
1325  {
1326  edgeTris[nTris++] = startFaces[startFacei];
1327  }
1328 
1329  forAll(endFaces, endFacei)
1330  {
1331  label facei = endFaces[endFacei];
1332 
1333  if ((facei != face1I) && (facei != face2I))
1334  {
1335  edgeTris[nTris++] = facei;
1336  }
1337  }
1338 }
1339 
1340 
1341 // Get all vertices connected to vertices of edge
1344  const triSurface& surf,
1345  const edge& e
1346 )
1347 {
1348  const edgeList& edges = surf.edges();
1349 
1350  label v1 = e.start();
1351  label v2 = e.end();
1352 
1353  // Get all vertices connected to v1 or v2 through an edge
1354  labelHashSet vertexNeighbours;
1355 
1356  const labelList& v1Edges = surf.pointEdges()[v1];
1357 
1358  forAll(v1Edges, v1EdgeI)
1359  {
1360  const edge& e = edges[v1Edges[v1EdgeI]];
1361  vertexNeighbours.insert(e.otherVertex(v1));
1362  }
1363 
1364  const labelList& v2Edges = surf.pointEdges()[v2];
1365 
1366  forAll(v2Edges, v2EdgeI)
1367  {
1368  const edge& e = edges[v2Edges[v2EdgeI]];
1369 
1370  label vertI = e.otherVertex(v2);
1371 
1372  vertexNeighbours.insert(vertI);
1373  }
1374  return vertexNeighbours.toc();
1375 }
1376 
1377 
1379 //void Foam::triSurfaceTools::orderVertices
1380 //(
1381 // const labelledTri& f,
1382 // const label v1,
1383 // const label v2,
1384 // label& vA,
1385 // label& vB
1386 //)
1387 //{
1388 // // Order v1, v2 in anticlockwise order.
1389 // bool reverse = false;
1390 //
1391 // if (f[0] == v1)
1392 // {
1393 // if (f[1] != v2)
1394 // {
1395 // reverse = true;
1396 // }
1397 // }
1398 // else if (f[1] == v1)
1399 // {
1400 // if (f[2] != v2)
1401 // {
1402 // reverse = true;
1403 // }
1404 // }
1405 // else
1406 // {
1407 // if (f[0] != v2)
1408 // {
1409 // reverse = true;
1410 // }
1411 // }
1412 //
1413 // if (reverse)
1414 // {
1415 // vA = v2;
1416 // vB = v1;
1417 // }
1418 // else
1419 // {
1420 // vA = v1;
1421 // vB = v2;
1422 // }
1423 //}
1424 
1425 
1426 // Get the other face using edgeI
1429  const triSurface& surf,
1430  const label facei,
1431  const label edgeI
1432 )
1433 {
1434  const labelList& myFaces = surf.edgeFaces()[edgeI];
1435 
1436  if (myFaces.size() != 2)
1437  {
1438  return -1;
1439  }
1440  else
1441  {
1442  if (facei == myFaces[0])
1443  {
1444  return myFaces[1];
1445  }
1446  else
1447  {
1448  return myFaces[0];
1449  }
1450  }
1451 }
1452 
1453 
1454 // Get the two edges on facei counterclockwise after edgeI
1457  const triSurface& surf,
1458  const label facei,
1459  const label edgeI,
1460  label& e1,
1461  label& e2
1462 )
1463 {
1464  const labelList& eFaces = surf.faceEdges()[facei];
1465 
1466  label i0 = findIndex(eFaces, edgeI);
1467 
1468  if (i0 == -1)
1469  {
1471  << "Edge " << surf.edges()[edgeI] << " not in face "
1472  << surf.localFaces()[facei] << abort(FatalError);
1473  }
1474 
1475  label i1 = eFaces.fcIndex(i0);
1476  label i2 = eFaces.fcIndex(i1);
1477 
1478  e1 = eFaces[i1];
1479  e2 = eFaces[i2];
1480 }
1481 
1482 
1483 // Get the two vertices on facei counterclockwise vertI
1486  const triSurface& surf,
1487  const label facei,
1488  const label vertI,
1489  label& vert1I,
1490  label& vert2I
1491 )
1492 {
1493  const labelledTri& f = surf.localFaces()[facei];
1494 
1495  if (vertI == f[0])
1496  {
1497  vert1I = f[1];
1498  vert2I = f[2];
1499  }
1500  else if (vertI == f[1])
1501  {
1502  vert1I = f[2];
1503  vert2I = f[0];
1504  }
1505  else if (vertI == f[2])
1506  {
1507  vert1I = f[0];
1508  vert2I = f[1];
1509  }
1510  else
1511  {
1513  << "Vertex " << vertI << " not in face " << f << abort(FatalError);
1514  }
1515 }
1516 
1517 
1518 // Get edge opposite vertex
1521  const triSurface& surf,
1522  const label facei,
1523  const label vertI
1524 )
1525 {
1526  const labelList& myEdges = surf.faceEdges()[facei];
1527 
1528  forAll(myEdges, myEdgeI)
1529  {
1530  label edgeI = myEdges[myEdgeI];
1531 
1532  const edge& e = surf.edges()[edgeI];
1533 
1534  if ((e.start() != vertI) && (e.end() != vertI))
1535  {
1536  return edgeI;
1537  }
1538  }
1539 
1541  << "Cannot find vertex " << vertI << " in edges of face " << facei
1542  << abort(FatalError);
1543 
1544  return -1;
1545 }
1546 
1547 
1548 // Get vertex opposite edge
1551  const triSurface& surf,
1552  const label facei,
1553  const label edgeI
1554 )
1555 {
1556  const triSurface::FaceType& f = surf.localFaces()[facei];
1557  const edge& e = surf.edges()[edgeI];
1558 
1559  forAll(f, fp)
1560  {
1561  label vertI = f[fp];
1562 
1563  if (vertI != e.start() && vertI != e.end())
1564  {
1565  return vertI;
1566  }
1567  }
1568 
1570  << "Cannot find vertex opposite edge " << edgeI << " vertices " << e
1571  << " in face " << facei << " vertices " << f << abort(FatalError);
1572 
1573  return -1;
1574 }
1575 
1576 
1577 // Returns edge label connecting v1, v2
1580  const triSurface& surf,
1581  const label v1,
1582  const label v2
1583 )
1584 {
1585  const labelList& v1Edges = surf.pointEdges()[v1];
1586 
1587  forAll(v1Edges, v1EdgeI)
1588  {
1589  label edgeI = v1Edges[v1EdgeI];
1590  const edge& e = surf.edges()[edgeI];
1591 
1592  if ((e.start() == v2) || (e.end() == v2))
1593  {
1594  return edgeI;
1595  }
1596  }
1597  return -1;
1598 }
1599 
1600 
1601 // Return index of triangle (or -1) using all three edges
1604  const triSurface& surf,
1605  const label e0I,
1606  const label e1I,
1607  const label e2I
1608 )
1609 {
1610  if ((e0I == e1I) || (e0I == e2I) || (e1I == e2I))
1611  {
1613  << "Duplicate edge labels : e0:" << e0I << " e1:" << e1I
1614  << " e2:" << e2I
1615  << abort(FatalError);
1616  }
1617 
1618  const labelList& eFaces = surf.edgeFaces()[e0I];
1619 
1620  forAll(eFaces, eFacei)
1621  {
1622  label facei = eFaces[eFacei];
1623 
1624  const labelList& myEdges = surf.faceEdges()[facei];
1625 
1626  if
1627  (
1628  (myEdges[0] == e1I)
1629  || (myEdges[1] == e1I)
1630  || (myEdges[2] == e1I)
1631  )
1632  {
1633  if
1634  (
1635  (myEdges[0] == e2I)
1636  || (myEdges[1] == e2I)
1637  || (myEdges[2] == e2I)
1638  )
1639  {
1640  return facei;
1641  }
1642  }
1643  }
1644  return -1;
1645 }
1646 
1647 
1648 // Collapse indicated edges. Return new tri surface.
1651  const triSurface& surf,
1652  const labelList& collapsableEdges
1653 )
1654 {
1655  pointField edgeMids(surf.nEdges());
1656 
1657  forAll(edgeMids, edgeI)
1658  {
1659  const edge& e = surf.edges()[edgeI];
1660 
1661  edgeMids[edgeI] =
1662  0.5
1663  * (
1664  surf.localPoints()[e.start()]
1665  + surf.localPoints()[e.end()]
1666  );
1667  }
1668 
1669 
1670  labelList faceStatus(surf.size(), ANYEDGE);
1671 
1673  // forAll(edges, edgeI)
1674  //{
1675  // const labelList& neighbours = edgeFaces[edgeI];
1676  //
1677  // if ((neighbours.size() != 2) && (neighbours.size() != 1))
1678  // {
1679  // FatalErrorInFunction
1680  // << abort(FatalError);
1681  // }
1682  //
1683  // if (neighbours.size() == 2)
1684  // {
1685  // if (surf[neighbours[0]].region() != surf[neighbours[1]].region())
1686  // {
1687  // // Neighbours on different regions. For now don't allow
1688  // // any collapse.
1689  // // Pout<< "protecting face " << neighbours[0]
1690  // // << ' ' << neighbours[1] << endl;
1691  // faceStatus[neighbours[0]] = NOEDGE;
1692  // faceStatus[neighbours[1]] = NOEDGE;
1693  // }
1694  // }
1695  //}
1696 
1697  return collapseEdges(surf, collapsableEdges, edgeMids, faceStatus);
1698 }
1699 
1700 
1701 // Collapse indicated edges. Return new tri surface.
1704  const triSurface& surf,
1705  const labelList& collapseEdgeLabels,
1706  const pointField& edgeMids,
1707  labelList& faceStatus
1708 )
1709 {
1710  const labelListList& edgeFaces = surf.edgeFaces();
1711  const pointField& localPoints = surf.localPoints();
1712  const edgeList& edges = surf.edges();
1713 
1714  // Storage for new points
1715  pointField newPoints(localPoints);
1716 
1717  // Map for old to new points
1718  labelList pointMap(localPoints.size());
1719  forAll(localPoints, pointi)
1720  {
1721  pointMap[pointi] = pointi;
1722  }
1723 
1724 
1725  // Do actual 'collapsing' of edges
1726 
1727  forAll(collapseEdgeLabels, collapseEdgeI)
1728  {
1729  const label edgeI = collapseEdgeLabels[collapseEdgeI];
1730 
1731  if ((edgeI < 0) || (edgeI >= surf.nEdges()))
1732  {
1734  << "Edge label outside valid range." << endl
1735  << "edge label:" << edgeI << endl
1736  << "total number of edges:" << surf.nEdges() << endl
1737  << abort(FatalError);
1738  }
1739 
1740  const labelList& neighbours = edgeFaces[edgeI];
1741 
1742  if (neighbours.size() == 2)
1743  {
1744  const label stat0 = faceStatus[neighbours[0]];
1745  const label stat1 = faceStatus[neighbours[1]];
1746 
1747  // Check faceStatus to make sure this one can be collapsed
1748  if
1749  (
1750  ((stat0 == ANYEDGE) || (stat0 == edgeI))
1751  && ((stat1 == ANYEDGE) || (stat1 == edgeI))
1752  )
1753  {
1754  const edge& e = edges[edgeI];
1755 
1756  // Set up mapping to 'collapse' points of edge
1757  if
1758  (
1759  (pointMap[e.start()] != e.start())
1760  || (pointMap[e.end()] != e.end())
1761  )
1762  {
1764  << "points already mapped. Double collapse." << endl
1765  << "edgeI:" << edgeI
1766  << " start:" << e.start()
1767  << " end:" << e.end()
1768  << " pointMap[start]:" << pointMap[e.start()]
1769  << " pointMap[end]:" << pointMap[e.end()]
1770  << abort(FatalError);
1771  }
1772 
1773  const label minVert = min(e.start(), e.end());
1774  pointMap[e.start()] = minVert;
1775  pointMap[e.end()] = minVert;
1776 
1777  // Move shared vertex to mid of edge
1778  newPoints[minVert] = edgeMids[edgeI];
1779 
1780  // Protect neighbouring faces
1781  protectNeighbours(surf, e.start(), faceStatus);
1782  protectNeighbours(surf, e.end(), faceStatus);
1783  protectNeighbours
1784  (
1785  surf,
1786  oppositeVertex(surf, neighbours[0], edgeI),
1787  faceStatus
1788  );
1789  protectNeighbours
1790  (
1791  surf,
1792  oppositeVertex(surf, neighbours[1], edgeI),
1793  faceStatus
1794  );
1795 
1796  // Mark all collapsing faces
1797  labelList collapseFaces =
1798  getCollapsedFaces
1799  (
1800  surf,
1801  edgeI
1802  ).toc();
1803 
1804  forAll(collapseFaces, collapseI)
1805  {
1806  faceStatus[collapseFaces[collapseI]] = COLLAPSED;
1807  }
1808  }
1809  }
1810  }
1811 
1812 
1813  // Storage for new triangles
1814  List<labelledTri> newTris(surf.size());
1815  label newTriI = 0;
1816 
1817  const List<labelledTri>& localFaces = surf.localFaces();
1818 
1819 
1820  // Get only non-collapsed triangles and renumber vertex labels.
1821  forAll(localFaces, facei)
1822  {
1823  const labelledTri& f = localFaces[facei];
1824 
1825  const label a = pointMap[f[0]];
1826  const label b = pointMap[f[1]];
1827  const label c = pointMap[f[2]];
1828 
1829  if
1830  (
1831  (a != b) && (a != c) && (b != c)
1832  && (faceStatus[facei] != COLLAPSED)
1833  )
1834  {
1835  // uncollapsed triangle
1836  newTris[newTriI++] = labelledTri(a, b, c, f.region());
1837  }
1838  else
1839  {
1840  // Pout<< "Collapsed triangle " << facei
1841  // << " vertices:" << f << endl;
1842  }
1843  }
1844  newTris.setSize(newTriI);
1845 
1846 
1847 
1848  // Pack faces
1849 
1850  triSurface tempSurf(newTris, surf.patches(), newPoints);
1851 
1852  return
1853  triSurface
1854  (
1855  tempSurf.localFaces(),
1856  tempSurf.patches(),
1857  tempSurf.localPoints()
1858  );
1859 }
1860 
1861 
1862 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1863 
1866  const triSurface& surf,
1867  const labelList& refineFaces
1868 )
1869 {
1870  List<refineType> refineStatus(surf.size(), NONE);
1871 
1872  // Mark & propagate refinement
1873  forAll(refineFaces, refineFacei)
1874  {
1875  calcRefineStatus(surf, refineFaces[refineFacei], refineStatus);
1876  }
1877 
1878  // Do actual refinement
1879  return doRefine(surf, refineStatus);
1880 }
1881 
1882 
1883 Foam::triSurface Foam::triSurfaceTools::greenRefine
1885  const triSurface& surf,
1886  const labelList& refineEdges
1887 )
1888 {
1889  // Storage for marking faces
1890  List<refineType> refineStatus(surf.size(), NONE);
1891 
1892  // Storage for new faces
1893  DynamicList<labelledTri> newFaces(0);
1894 
1895  pointField newPoints(surf.localPoints());
1896  newPoints.setSize(surf.nPoints() + surf.nEdges());
1897  label newPointi = surf.nPoints();
1898 
1899 
1900  // Refine edges
1901  forAll(refineEdges, refineEdgeI)
1902  {
1903  label edgeI = refineEdges[refineEdgeI];
1904 
1905  const labelList& myFaces = surf.edgeFaces()[edgeI];
1906 
1907  bool neighbourIsRefined= false;
1908 
1909  forAll(myFaces, myFacei)
1910  {
1911  if (refineStatus[myFaces[myFacei]] != NONE)
1912  {
1913  neighbourIsRefined = true;
1914  }
1915  }
1916 
1917  // Only refine if none of the faces is refined
1918  if (!neighbourIsRefined)
1919  {
1920  // Refine edge
1921  const edge& e = surf.edges()[edgeI];
1922 
1923  point mid =
1924  0.5
1925  * (
1926  surf.localPoints()[e.start()]
1927  + surf.localPoints()[e.end()]
1928  );
1929 
1930  newPoints[newPointi] = mid;
1931 
1932  // Refine faces using edge
1933  forAll(myFaces, myFacei)
1934  {
1935  // Add faces to newFaces
1936  greenRefine
1937  (
1938  surf,
1939  myFaces[myFacei],
1940  edgeI,
1941  newPointi,
1942  newFaces
1943  );
1944 
1945  // Mark as refined
1946  refineStatus[myFaces[myFacei]] = GREEN;
1947  }
1948 
1949  newPointi++;
1950  }
1951  }
1952 
1953  // Add unrefined faces
1954  forAll(surf.localFaces(), facei)
1955  {
1956  if (refineStatus[facei] == NONE)
1957  {
1958  newFaces.append(surf.localFaces()[facei]);
1959  }
1960  }
1961 
1962  newFaces.shrink();
1963  newPoints.setSize(newPointi);
1964 
1965  return triSurface(newFaces, surf.patches(), newPoints, true);
1966 }
1967 
1968 
1969 
1970 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1971 // Geometric helper functions:
1972 
1973 
1974 // Returns element in edgeIndices with minimum length
1977  const triSurface& surf,
1978  const labelList& edgeIndices
1979 )
1980 {
1981  scalar minLength = great;
1982  label minIndex = -1;
1983  forAll(edgeIndices, i)
1984  {
1985  const edge& e = surf.edges()[edgeIndices[i]];
1986 
1987  scalar length =
1988  mag
1989  (
1990  surf.localPoints()[e.end()]
1991  - surf.localPoints()[e.start()]
1992  );
1993 
1994  if (length < minLength)
1995  {
1996  minLength = length;
1997  minIndex = i;
1998  }
1999  }
2000  return edgeIndices[minIndex];
2001 }
2002 
2003 
2004 // Returns element in edgeIndices with maximum length
2007  const triSurface& surf,
2008  const labelList& edgeIndices
2009 )
2010 {
2011  scalar maxLength = -great;
2012  label maxIndex = -1;
2013  forAll(edgeIndices, i)
2014  {
2015  const edge& e = surf.edges()[edgeIndices[i]];
2016 
2017  scalar length =
2018  mag
2019  (
2020  surf.localPoints()[e.end()]
2021  - surf.localPoints()[e.start()]
2022  );
2023 
2024  if (length > maxLength)
2025  {
2026  maxLength = length;
2027  maxIndex = i;
2028  }
2029  }
2030  return edgeIndices[maxIndex];
2031 }
2032 
2033 
2034 // Merge points and reconstruct surface
2037  const triSurface& surf,
2038  const scalar mergeTol
2039 )
2040 {
2041  pointField newPoints(surf.nPoints());
2042 
2043  labelList pointMap(surf.nPoints());
2044 
2045  bool hasMerged = Foam::mergePoints
2046  (
2047  surf.localPoints(),
2048  mergeTol,
2049  false,
2050  pointMap,
2051  newPoints
2052  );
2053 
2054  if (hasMerged)
2055  {
2056  // Pack the triangles
2057 
2058  // Storage for new triangles
2059  List<labelledTri> newTriangles(surf.size());
2060  label newTriangleI = 0;
2061 
2062  forAll(surf, facei)
2063  {
2064  const labelledTri& f = surf.localFaces()[facei];
2065 
2066  label newA = pointMap[f[0]];
2067  label newB = pointMap[f[1]];
2068  label newC = pointMap[f[2]];
2069 
2070  if ((newA != newB) && (newA != newC) && (newB != newC))
2071  {
2072  newTriangles[newTriangleI++] =
2073  labelledTri(newA, newB, newC, f.region());
2074  }
2075  }
2076  newTriangles.setSize(newTriangleI);
2077 
2078  return triSurface
2079  (
2080  newTriangles,
2081  surf.patches(),
2082  newPoints,
2083  true // reuse storage
2084  );
2085  }
2086  else
2087  {
2088  return surf;
2089  }
2090 }
2091 
2092 
2093 // Calculate normal on triangle
2096  const triSurface& surf,
2097  const label nearestFacei,
2098  const point& nearestPt
2099 )
2100 {
2101  const triSurface::FaceType& f = surf[nearestFacei];
2102  const pointField& points = surf.points();
2103 
2104  label nearType, nearLabel;
2105 
2106  f.nearestPointClassify(nearestPt, points, nearType, nearLabel);
2107 
2108  if (nearType == triPointRef::NONE)
2109  {
2110  // Nearest to face
2111  return surf.faceNormals()[nearestFacei];
2112  }
2113  else if (nearType == triPointRef::EDGE)
2114  {
2115  // Nearest to edge. Assume order of faceEdges same as face vertices.
2116  label edgeI = surf.faceEdges()[nearestFacei][nearLabel];
2117 
2118  // Calculate edge normal by averaging face normals
2119  const labelList& eFaces = surf.edgeFaces()[edgeI];
2120 
2121  vector edgeNormal(Zero);
2122 
2123  forAll(eFaces, i)
2124  {
2125  edgeNormal += surf.faceNormals()[eFaces[i]];
2126  }
2127  return edgeNormal/(mag(edgeNormal) + vSmall);
2128  }
2129  else
2130  {
2131  // Nearest to point
2132  const triSurface::FaceType& localF = surf.localFaces()[nearestFacei];
2133  return surf.pointNormals()[localF[nearLabel]];
2134  }
2135 }
2136 
2137 
2140  const triSurface& surf,
2141  const point& sample,
2142  const point& nearestPoint,
2143  const label edgeI
2144 )
2145 {
2146  const labelList& eFaces = surf.edgeFaces()[edgeI];
2147 
2148  if (eFaces.size() != 2)
2149  {
2150  // Surface not closed.
2151  return UNKNOWN;
2152  }
2153  else
2154  {
2155  const vectorField& faceNormals = surf.faceNormals();
2156 
2157  // Compare to bisector. This is actually correct since edge is
2158  // nearest so there is a knife-edge.
2159 
2160  vector n = 0.5*(faceNormals[eFaces[0]] + faceNormals[eFaces[1]]);
2161 
2162  if (((sample - nearestPoint) & n) > 0)
2163  {
2164  return OUTSIDE;
2165  }
2166  else
2167  {
2168  return INSIDE;
2169  }
2170  }
2171 }
2172 
2173 
2174 // Calculate normal on triangle
2177  const triSurface& surf,
2178  const point& sample,
2179  const label nearestFacei
2180 )
2181 {
2182  const triSurface::FaceType& f = surf[nearestFacei];
2183  const pointField& points = surf.points();
2184 
2185  // Find where point is on face
2186  label nearType, nearLabel;
2187 
2188  pointHit pHit = f.nearestPointClassify(sample, points, nearType, nearLabel);
2189 
2190  const point& nearestPoint(pHit.rawPoint());
2191 
2192  if (nearType == triPointRef::NONE)
2193  {
2194  vector sampleNearestVec = (sample - nearestPoint);
2195 
2196  // Nearest to face interior. Use faceNormal to determine side
2197  scalar c = sampleNearestVec & surf.faceNormals()[nearestFacei];
2198 
2199  // // If the sample is essentially on the face, do not check for
2200  // // it being perpendicular.
2201 
2202  // scalar magSampleNearestVec = mag(sampleNearestVec);
2203 
2204  // if (magSampleNearestVec > small)
2205  // {
2206  // c /= magSampleNearestVec*mag(surf.faceNormals()[nearestFacei]);
2207 
2208  // if (mag(c) < 0.99)
2209  // {
2210  // FatalErrorInFunction
2211  // << "nearestPoint identified as being on triangle face "
2212  // << "but vector from nearestPoint to sample is not "
2213  // << "perpendicular to the normal." << nl
2214  // << "sample: " << sample << nl
2215  // << "nearestPoint: " << nearestPoint << nl
2216  // << "sample - nearestPoint: "
2217  // << sample - nearestPoint << nl
2218  // << "normal: " << surf.faceNormals()[nearestFacei] << nl
2219  // << "mag(sample - nearestPoint): "
2220  // << mag(sample - nearestPoint) << nl
2221  // << "normalised dot product: " << c << nl
2222  // << "triangle vertices: " << nl
2223  // << " " << points[f[0]] << nl
2224  // << " " << points[f[1]] << nl
2225  // << " " << points[f[2]] << nl
2226  // << abort(FatalError);
2227  // }
2228  // }
2229 
2230  if (c > 0)
2231  {
2232  return OUTSIDE;
2233  }
2234  else
2235  {
2236  return INSIDE;
2237  }
2238  }
2239  else if (nearType == triPointRef::EDGE)
2240  {
2241  // Nearest to edge nearLabel. Note that this can only be a knife-edge
2242  // situation since otherwise the nearest point could never be the edge.
2243 
2244  // Get the edge. Assume order of faceEdges same as face vertices.
2245  label edgeI = surf.faceEdges()[nearestFacei][nearLabel];
2246 
2247  // if (debug)
2248  // {
2249  // // Check order of faceEdges same as face vertices.
2250  // const edge& e = surf.edges()[edgeI];
2251  // const labelList& meshPoints = surf.meshPoints();
2252  // const edge meshEdge(meshPoints[e[0]], meshPoints[e[1]]);
2253 
2254  // if
2255  // (
2256  // meshEdge
2257  // != edge(f[nearLabel], f[f.fcIndex(nearLabel)])
2258  // )
2259  // {
2260  // FatalErrorInFunction
2261  // << "Edge:" << edgeI << " local vertices:" << e
2262  // << " mesh vertices:" << meshEdge
2263  // << " not at position " << nearLabel
2264  // << " in face " << f
2265  // << abort(FatalError);
2266  // }
2267  // }
2268 
2269  return edgeSide(surf, sample, nearestPoint, edgeI);
2270  }
2271  else
2272  {
2273  // Nearest to point. Could use pointNormal here but is not correct.
2274  // Instead determine which edge using point is nearest and use test
2275  // above (nearType == triPointRef::EDGE).
2276 
2277 
2278  const triSurface::FaceType& localF = surf.localFaces()[nearestFacei];
2279  label nearPointi = localF[nearLabel];
2280 
2281  const edgeList& edges = surf.edges();
2282  const pointField& localPoints = surf.localPoints();
2283  const point& base = localPoints[nearPointi];
2284 
2285  const labelList& pEdges = surf.pointEdges()[nearPointi];
2286 
2287  scalar minDistSqr = Foam::sqr(great);
2288  label minEdgeI = -1;
2289 
2290  forAll(pEdges, i)
2291  {
2292  label edgeI = pEdges[i];
2293 
2294  const edge& e = edges[edgeI];
2295 
2296  label otherPointi = e.otherVertex(nearPointi);
2297 
2298  // Get edge normal.
2299  vector eVec(localPoints[otherPointi] - base);
2300  scalar magEVec = mag(eVec);
2301 
2302  if (magEVec > vSmall)
2303  {
2304  eVec /= magEVec;
2305 
2306  // Get point along vector and determine closest.
2307  const point perturbPoint = base + eVec;
2308 
2309  scalar distSqr = Foam::magSqr(sample - perturbPoint);
2310 
2311  if (distSqr < minDistSqr)
2312  {
2313  minDistSqr = distSqr;
2314  minEdgeI = edgeI;
2315  }
2316  }
2317  }
2318 
2319  if (minEdgeI == -1)
2320  {
2322  << "Problem: did not find edge closer than " << minDistSqr
2323  << abort(FatalError);
2324  }
2325 
2326  return edgeSide(surf, sample, nearestPoint, minEdgeI);
2327  }
2328 }
2329 
2330 
2333  const polyBoundaryMesh& bMesh,
2334  const labelHashSet& includePatches,
2335  const bool verbose
2336 )
2337 {
2338  const polyMesh& mesh = bMesh.mesh();
2339 
2340  // Storage for surfaceMesh. Size estimate.
2341  DynamicList<labelledTri> triangles
2342  (
2343  mesh.nFaces() - mesh.nInternalFaces()
2344  );
2345 
2346  polygonTriangulate triEngine;
2347 
2348  label newPatchi = 0;
2349 
2350  forAllConstIter(labelHashSet, includePatches, iter)
2351  {
2352  const label patchi = iter.key();
2353  const polyPatch& patch = bMesh[patchi];
2354  const pointField& points = patch.points();
2355 
2356  label nTriTotal = 0;
2357 
2358  forAll(patch, patchFacei)
2359  {
2360  const face& f = patch[patchFacei];
2361 
2362  triEngine.triangulate(UIndirectList<point>(points, f));
2363 
2364  forAll(triEngine.triPoints(), triFacei)
2365  {
2366  triangles.append
2367  (
2368  labelledTri(triEngine.triPoints(triFacei, f), newPatchi)
2369  );
2370 
2371  nTriTotal++;
2372  }
2373  }
2374 
2375  if (verbose)
2376  {
2377  Pout<< patch.name() << " : generated " << nTriTotal
2378  << " triangles from " << patch.size() << " faces with"
2379  << " new patchid " << newPatchi << endl;
2380  }
2381 
2382  newPatchi++;
2383  }
2384  triangles.shrink();
2385 
2386  // Create globally numbered tri surface
2387  triSurface rawSurface(triangles, mesh.points());
2388 
2389  // Create locally numbered tri surface
2390  triSurface surface
2391  (
2392  rawSurface.localFaces(),
2393  rawSurface.localPoints()
2394  );
2395 
2396  // Add patch names to surface
2397  surface.patches().setSize(newPatchi);
2398 
2399  newPatchi = 0;
2400 
2401  forAllConstIter(labelHashSet, includePatches, iter)
2402  {
2403  const label patchi = iter.key();
2404  const polyPatch& patch = bMesh[patchi];
2405 
2406  surface.patches()[newPatchi].name() = patch.name();
2407  surface.patches()[newPatchi].geometricType() = patch.type();
2408 
2409  newPatchi++;
2410  }
2411 
2412  return surface;
2413 }
2414 
2415 
2418  const polyBoundaryMesh& bMesh,
2419  const labelHashSet& includePatches,
2420  const boundBox& bBox,
2421  const bool verbose
2422 )
2423 {
2424  const polyMesh& mesh = bMesh.mesh();
2425 
2426  // Storage for surfaceMesh. Size estimate.
2427  DynamicList<labelledTri> triangles
2428  (
2429  mesh.nFaces() - mesh.nInternalFaces()
2430  );
2431 
2432  polygonTriangulate triEngine;
2433 
2434  label newPatchi = 0;
2435 
2436  forAllConstIter(labelHashSet, includePatches, iter)
2437  {
2438  const label patchi = iter.key();
2439  const polyPatch& patch = bMesh[patchi];
2440  const pointField& points = patch.points();
2441 
2442  label nTriTotal = 0;
2443 
2444  forAll(patch, patchFacei)
2445  {
2446  const face& f = patch[patchFacei];
2447 
2448  if (bBox.containsAny(points, f))
2449  {
2450  triEngine.triangulate(UIndirectList<point>(points, f));
2451 
2452  forAll(triEngine.triPoints(), triFacei)
2453  {
2454  triangles.append
2455  (
2456  labelledTri(triEngine.triPoints(triFacei, f), newPatchi)
2457  );
2458 
2459  nTriTotal++;
2460  }
2461  }
2462  }
2463 
2464  if (verbose)
2465  {
2466  Pout<< patch.name() << " : generated " << nTriTotal
2467  << " triangles from " << patch.size() << " faces with"
2468  << " new patchid " << newPatchi << endl;
2469  }
2470 
2471  newPatchi++;
2472  }
2473  triangles.shrink();
2474 
2475  // Create globally numbered tri surface
2476  triSurface rawSurface(triangles, mesh.points());
2477 
2478  // Create locally numbered tri surface
2479  triSurface surface
2480  (
2481  rawSurface.localFaces(),
2482  rawSurface.localPoints()
2483  );
2484 
2485  // Add patch names to surface
2486  surface.patches().setSize(newPatchi);
2487 
2488  newPatchi = 0;
2489 
2490  forAllConstIter(labelHashSet, includePatches, iter)
2491  {
2492  const label patchi = iter.key();
2493  const polyPatch& patch = bMesh[patchi];
2494 
2495  surface.patches()[newPatchi].name() = patch.name();
2496  surface.patches()[newPatchi].geometricType() = patch.type();
2497 
2498  newPatchi++;
2499  }
2500 
2501  return surface;
2502 }
2503 
2504 
2505 // triangulation of boundaryMesh
2508  const polyBoundaryMesh& bMesh,
2509  const labelHashSet& includePatches,
2510  const bool verbose
2511 )
2512 {
2513  const polyMesh& mesh = bMesh.mesh();
2514 
2515  // Storage for new points = meshpoints + face centres.
2516  const pointField& points = mesh.points();
2517  const pointField& faceCentres = mesh.faceCentres();
2518 
2519  pointField newPoints(points.size() + faceCentres.size());
2520 
2521  label newPointi = 0;
2522 
2523  forAll(points, pointi)
2524  {
2525  newPoints[newPointi++] = points[pointi];
2526  }
2527  forAll(faceCentres, facei)
2528  {
2529  newPoints[newPointi++] = faceCentres[facei];
2530  }
2531 
2532 
2533  // Count number of faces.
2534  DynamicList<labelledTri> triangles
2535  (
2536  mesh.nFaces() - mesh.nInternalFaces()
2537  );
2538 
2539  label newPatchi = 0;
2540 
2541  forAllConstIter(labelHashSet, includePatches, iter)
2542  {
2543  const label patchi = iter.key();
2544  const polyPatch& patch = bMesh[patchi];
2545 
2546  label nTriTotal = 0;
2547 
2548  forAll(patch, patchFacei)
2549  {
2550  // Face in global coords.
2551  const face& f = patch[patchFacei];
2552 
2553  // Index in newPointi of face centre.
2554  label fc = points.size() + patchFacei + patch.start();
2555 
2556  forAll(f, fp)
2557  {
2558  label fp1 = f.fcIndex(fp);
2559 
2560  triangles.append(labelledTri(f[fp], f[fp1], fc, newPatchi));
2561 
2562  nTriTotal++;
2563  }
2564  }
2565 
2566  if (verbose)
2567  {
2568  Pout<< patch.name() << " : generated " << nTriTotal
2569  << " triangles from " << patch.size() << " faces with"
2570  << " new patchid " << newPatchi << endl;
2571  }
2572 
2573  newPatchi++;
2574  }
2575  triangles.shrink();
2576 
2577 
2578  // Create globally numbered tri surface
2579  triSurface rawSurface(triangles, newPoints);
2580 
2581  // Create locally numbered tri surface
2582  triSurface surface
2583  (
2584  rawSurface.localFaces(),
2585  rawSurface.localPoints()
2586  );
2587 
2588  // Add patch names to surface
2589  surface.patches().setSize(newPatchi);
2590 
2591  newPatchi = 0;
2592 
2593  forAllConstIter(labelHashSet, includePatches, iter)
2594  {
2595  const label patchi = iter.key();
2596  const polyPatch& patch = bMesh[patchi];
2597 
2598  surface.patches()[newPatchi].name() = patch.name();
2599  surface.patches()[newPatchi].geometricType() = patch.type();
2600 
2601  newPatchi++;
2602  }
2603 
2604  return surface;
2605 }
2606 
2607 
2609 {
2610  // Vertices in geompack notation. Note that could probably just use
2611  // pts.begin() if double precision.
2612  List<doubleScalar> geompackVertices(2*pts.size());
2613  label doubleI = 0;
2614  forAll(pts, i)
2615  {
2616  geompackVertices[doubleI++] = pts[i][0];
2617  geompackVertices[doubleI++] = pts[i][1];
2618  }
2619 
2620  // Storage for triangles
2621  int m2 = 3;
2622  List<int> triangle_node(m2*3*pts.size());
2623  List<int> triangle_neighbor(m2*3*pts.size());
2624 
2625  // Triangulate
2626  int nTris = 0;
2627  int err = dtris2
2628  (
2629  pts.size(),
2630  geompackVertices.begin(),
2631  &nTris,
2632  triangle_node.begin(),
2633  triangle_neighbor.begin()
2634  );
2635 
2636  if (err != 0)
2637  {
2639  << "Failed dtris2 with vertices:" << pts.size()
2640  << abort(FatalError);
2641  }
2642 
2643  // Trim
2644  triangle_node.setSize(3*nTris);
2645  triangle_neighbor.setSize(3*nTris);
2646 
2647  // Convert to triSurface.
2648  List<labelledTri> faces(nTris);
2649 
2650  forAll(faces, i)
2651  {
2652  faces[i] = labelledTri
2653  (
2654  triangle_node[3*i]-1,
2655  triangle_node[3*i+1]-1,
2656  triangle_node[3*i+2]-1,
2657  0
2658  );
2659  }
2660 
2661  pointField points(pts.size());
2662  forAll(pts, i)
2663  {
2664  points[i][0] = pts[i][0];
2665  points[i][1] = pts[i][1];
2666  points[i][2] = 0.0;
2667  }
2668 
2669  return triSurface(faces, points);
2670 }
2671 
2672 
2675  const triPointRef& tri,
2676  const point& p,
2677  FixedList<scalar, 3>& weights
2678 )
2679 {
2680  // calculate triangle edge vectors and triangle face normal
2681  // the 'i':th edge is opposite node i
2683  edge[0] = tri.c()-tri.b();
2684  edge[1] = tri.a()-tri.c();
2685  edge[2] = tri.b()-tri.a();
2686 
2687  vector triangleFaceNormal = edge[1] ^ edge[2];
2688 
2689  // calculate edge normal (pointing inwards)
2690  FixedList<vector, 3> normal;
2691  for (label i=0; i<3; i++)
2692  {
2693  normal[i] = triangleFaceNormal ^ edge[i];
2694  normal[i] /= mag(normal[i]) + vSmall;
2695  }
2696 
2697  weights[0] = ((p-tri.b()) & normal[0]) / max(vSmall, normal[0] & edge[1]);
2698  weights[1] = ((p-tri.c()) & normal[1]) / max(vSmall, normal[1] & edge[2]);
2699  weights[2] = ((p-tri.a()) & normal[2]) / max(vSmall, normal[2] & edge[0]);
2700 }
2701 
2702 
2703 // Calculate weighting factors from samplePts to triangle it is in.
2704 // Uses linear search.
2707  const triSurface& s,
2708  const pointField& samplePts,
2709  List<FixedList<label, 3>>& allVerts,
2710  List<FixedList<scalar, 3>>& allWeights
2711 )
2712 {
2713  allVerts.setSize(samplePts.size());
2714  allWeights.setSize(samplePts.size());
2715 
2716  const pointField& points = s.points();
2717 
2718  forAll(samplePts, i)
2719  {
2720  const point& samplePt = samplePts[i];
2721 
2722 
2723  FixedList<label, 3>& verts = allVerts[i];
2724  FixedList<scalar, 3>& weights = allWeights[i];
2725 
2726  scalar minDistance = great;
2727 
2728  forAll(s, facei)
2729  {
2730  const labelledTri& f = s[facei];
2731 
2732  triPointRef tri(f.tri(points));
2733 
2734  label nearType, nearLabel;
2735 
2736  pointHit nearest = tri.nearestPointClassify
2737  (
2738  samplePt,
2739  nearType,
2740  nearLabel
2741  );
2742 
2743  if (nearest.hit())
2744  {
2745  // samplePt inside triangle
2746  verts[0] = f[0];
2747  verts[1] = f[1];
2748  verts[2] = f[2];
2749 
2750  calcInterpolationWeights(tri, nearest.rawPoint(), weights);
2751 
2752  // Pout<< "calcScalingFactors : samplePt:" << samplePt
2753  // << " inside triangle:" << facei
2754  // << " verts:" << verts
2755  // << " weights:" << weights
2756  // << endl;
2757 
2758  break;
2759  }
2760  else if (nearest.distance() < minDistance)
2761  {
2762  minDistance = nearest.distance();
2763 
2764  // Outside triangle. Store nearest.
2765 
2766  if (nearType == triPointRef::POINT)
2767  {
2768  verts[0] = f[nearLabel];
2769  weights[0] = 1;
2770  verts[1] = -1;
2771  weights[1] = -great;
2772  verts[2] = -1;
2773  weights[2] = -great;
2774 
2775  // Pout<< "calcScalingFactors : samplePt:" << samplePt
2776  // << " distance:" << nearest.distance()
2777  // << " from point:" << points[f[nearLabel]]
2778  // << endl;
2779  }
2780  else if (nearType == triPointRef::EDGE)
2781  {
2782  verts[0] = f[nearLabel];
2783  verts[1] = f[f.fcIndex(nearLabel)];
2784  verts[2] = -1;
2785 
2786  const point& p0 = points[verts[0]];
2787  const point& p1 = points[verts[1]];
2788 
2789  scalar s = min
2790  (
2791  1,
2792  max
2793  (
2794  0,
2795  mag(nearest.rawPoint() - p0)/mag(p1 - p0)
2796  )
2797  );
2798 
2799  // Interpolate
2800  weights[0] = 1 - s;
2801  weights[1] = s;
2802  weights[2] = -great;
2803 
2804  // Pout<< "calcScalingFactors : samplePt:" << samplePt
2805  // << " distance:" << nearest.distance()
2806  // << " from edge:" << p0 << p1 << " s:" << s
2807  // << endl;
2808  }
2809  else
2810  {
2811  // triangle. Can only happen because of truncation errors.
2812  verts[0] = f[0];
2813  verts[1] = f[1];
2814  verts[2] = f[2];
2815 
2816  calcInterpolationWeights(tri, nearest.rawPoint(), weights);
2817 
2818  // Pout<< "calcScalingFactors : samplePt:" << samplePt
2819  // << " distance:" << nearest.distance()
2820  // << " to verts:" << verts
2821  // << " weights:" << weights
2822  // << endl;
2823 
2824  break;
2825  }
2826  }
2827  }
2828  }
2829 }
2830 
2831 
2832 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
2833 // Tracking:
2834 
2835 
2836 // Test point on surface to see if is on face,edge or point.
2839  const triSurface& s,
2840  const label triI,
2841  const point& trianglePoint
2842 )
2843 {
2844  surfaceLocation nearest;
2845 
2846  // Nearest point could be on point or edge. Retest.
2847  label index, elemType;
2848  // bool inside =
2849  triPointRef(s[triI].tri(s.points())).classify
2850  (
2851  trianglePoint,
2852  elemType,
2853  index
2854  );
2855 
2856  nearest.setPoint(trianglePoint);
2857 
2858  if (elemType == triPointRef::NONE)
2859  {
2860  nearest.setHit();
2861  nearest.setIndex(triI);
2862  nearest.elementType() = triPointRef::NONE;
2863  }
2864  else if (elemType == triPointRef::EDGE)
2865  {
2866  nearest.setMiss();
2867  nearest.setIndex(s.faceEdges()[triI][index]);
2868  nearest.elementType() = triPointRef::EDGE;
2869  }
2870  else // if (elemType == triPointRef::POINT)
2871  {
2872  nearest.setMiss();
2873  nearest.setIndex(s.localFaces()[triI][index]);
2874  nearest.elementType() = triPointRef::POINT;
2875  }
2876 
2877  return nearest;
2878 }
2879 
2880 
2883  const triSurface& s,
2884  const surfaceLocation& start,
2885  const surfaceLocation& end,
2886  const plane& cutPlane
2887 )
2888 {
2889  // Start off from starting point
2890  surfaceLocation nearest = start;
2891  nearest.setMiss();
2892 
2893  // See if in same triangle as endpoint. If so snap.
2894  snapToEnd(s, end, nearest);
2895 
2896  if (!nearest.hit())
2897  {
2898  // Not yet at end point
2899 
2900  if (start.elementType() == triPointRef::NONE)
2901  {
2902  // Start point is inside triangle. Trivial cases already handled
2903  // above.
2904 
2905  // end point is on edge or point so cross current triangle to
2906  // see which edge is cut.
2907 
2908  nearest = cutEdge
2909  (
2910  s,
2911  start.index(), // triangle
2912  -1, // excludeEdge
2913  -1, // excludePoint
2914  start.rawPoint(),
2915  cutPlane,
2916  end.rawPoint()
2917  );
2918  nearest.elementType() = triPointRef::EDGE;
2919  nearest.triangle() = start.index();
2920  nearest.setMiss();
2921  }
2922  else if (start.elementType() == triPointRef::EDGE)
2923  {
2924  // Pick connected triangle that is most in direction.
2925  const labelList& eFaces = s.edgeFaces()[start.index()];
2926 
2927  nearest = visitFaces
2928  (
2929  s,
2930  eFaces,
2931  start,
2932  start.index(), // excludeEdgeI
2933  -1, // excludePointi
2934  end,
2935  cutPlane
2936  );
2937  }
2938  else // start.elementType() == triPointRef::POINT
2939  {
2940  const labelList& pFaces = s.pointFaces()[start.index()];
2941 
2942  nearest = visitFaces
2943  (
2944  s,
2945  pFaces,
2946  start,
2947  -1, // excludeEdgeI
2948  start.index(), // excludePointi
2949  end,
2950  cutPlane
2951  );
2952  }
2953  snapToEnd(s, end, nearest);
2954  }
2955  return nearest;
2956 }
2957 
2958 
2961  const triSurface& s,
2962  const surfaceLocation& endInfo,
2963  const plane& cutPlane,
2964  surfaceLocation& hitInfo
2965 )
2966 {
2967  // OFstream str("track.obj");
2968  // label vertI = 0;
2969  // meshTools::writeOBJ(str, hitInfo.rawPoint());
2970  // vertI++;
2971 
2972  // Track across surface.
2973  while (true)
2974  {
2975  // Pout<< "Tracking from:" << nl
2976  // << " " << hitInfo.info()
2977  // << endl;
2978 
2979  hitInfo = trackToEdge
2980  (
2981  s,
2982  hitInfo,
2983  endInfo,
2984  cutPlane
2985  );
2986 
2987  // meshTools::writeOBJ(str, hitInfo.rawPoint());
2988  // vertI++;
2989  // str<< "l " << vertI-1 << ' ' << vertI << nl;
2990 
2991  // Pout<< "Tracked to:" << nl
2992  // << " " << hitInfo.info() << endl;
2993 
2994  if (hitInfo.hit() || hitInfo.triangle() == -1)
2995  {
2996  break;
2997  }
2998  }
2999 }
3000 
3001 
3002 // ************************************************************************* //
const labelListList & pointEdges() const
Return point-edge addressing.
static labelList getVertexVertices(const triSurface &surf, const edge &e)
Get all vertices (local numbering) connected to vertices of edge.
label nPoints() const
Return number of points supporting patch faces.
A HashTable with keys but without contents.
Definition: HashSet.H:59
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A triangle primitive used to calculate face areas and swept volumes.
Definition: triangle.H:57
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
const word & name() const
Return name.
const word & name() const
Return name.
Definition: IOobject.H:303
A class for handling file names.
Definition: fileName.H:79
const Point & c() const
Return third vertex.
Definition: triangleI.H:82
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
label nInternalFaces() const
label nFaces() const
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Output to file stream.
Definition: OFstream.H:82
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
void setIndex(const label index)
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
pointHit nearestPointClassify(const point &p, label &nearType, label &nearLabel) const
Find the nearest point to p on the triangle and classify it:
Definition: triangleI.H:518
static void track(const triSurface &, const surfaceLocation &endInfo, const plane &cutPlane, surfaceLocation &hitInfo)
Track from edge to edge across surface. Uses trackToEdge.
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:58
const Cmpt & z() const
Definition: VectorI.H:87
static triSurface redGreenRefine(const triSurface &surf, const labelList &refineFaces)
Refine face by splitting all edges. Neighbouring face is.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
label otherVertex(const label a) const
Given one vertex, return the other.
Definition: edgeI.H:140
static label maxEdge(const triSurface &surf, const labelList &edgeIndices)
Returns element in edgeIndices with minimum length.
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
Definition: plane.H:60
const Field< PointType > & localPoints() const
Return pointField of points in patch.
static vector surfaceNormal(const triSurface &surf, const label nearestFacei, const point &nearestPt)
Triangle (unit) normal. If nearest point to triangle on edge use.
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:111
void setPoint(const Point &p)
const dimensionedScalar c
Speed of light in a vacuum.
static const label NOEDGE
static triSurface delaunay2D(const List< vector2D > &)
Do unconstrained Delaunay of points. Returns triSurface with 3D.
static triSurface collapseEdges(const triSurface &surf, const labelList &collapsableEdges)
Create new triSurface by collapsing edges to edge mids.
static label otherFace(const triSurface &surf, const label facei, const label edgeI)
Get face connected to edge not facei.
const Cmpt & y() const
Definition: VectorI.H:81
triSurface triangulateFaceCentre(const polyBoundaryMesh &mBesh, const labelHashSet &includePatches, const bool verbose=false)
Face-centre triangulation of (selected patches of) boundaryMesh.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1131
static label oppositeVertex(const triSurface &surf, const label facei, const label edgeI)
Get vertex (local numbering) opposite edge.
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:211
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))
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
static triSurface mergePoints(const triSurface &surf, const scalar mergeTol)
Merge points within distance.
dynamicFvMesh & mesh
triPointRef::proxType & elementType()
const pointField & points
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
List< edge > edgeList
Definition: edgeList.H:38
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
static sideType edgeSide(const triSurface &surf, const point &sample, const point &nearestPoint, const label edgeI)
If nearest point is on edgeI, determine on which side of surface.
static label minEdge(const triSurface &surf, const labelList &edgeIndices)
Returns element in edgeIndices with minimum length.
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:178
iterator begin()
Return an iterator to begin traversing the UList.
Definition: UListI.H:216
bool hit() const
Is there a hit.
const Field< PointType > & points() const
Return reference to global points.
static void otherVertices(const triSurface &surf, const label facei, const label vertI, label &vert1I, label &vert2I)
Get the two vertices (local numbering) on facei counterclockwise.
const polyMesh & mesh() const
Return the mesh reference.
const labelListList & edgeFaces() const
Return edge-face addressing.
bool hit() const
Is there a hit.
Definition: PointHit.H:120
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
static const label COLLAPSED
Triangle with additional region number.
Definition: labelledTri.H:57
triPointRef tri(const pointField &) const
Return the triangle.
Definition: triFaceI.H:152
const Field< PointType > & faceNormals() const
Return face normals for patch.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
errorManip< error > abort(error &err)
Definition: errorManip.H:131
static surfaceLocation trackToEdge(const triSurface &, const surfaceLocation &start, const surfaceLocation &end, const plane &cutPlane)
Track on surface to get closer to point.
const Cmpt & x() const
Definition: VectorI.H:75
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: FixedListI.H:125
const Point & rawPoint() const
Return point with no checking.
Foam::polyBoundaryMesh.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< small) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &mergedCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
const labelListList & pointFaces() const
Return point-face addressing.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
const Point & rawPoint() const
Return point with no checking.
Definition: PointHit.H:158
label nEdges() const
Return number of edges in patch.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:252
sideType
On which side of surface.
static const char nl
Definition: Ostream.H:260
Merge points. See below.
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
static void calcInterpolationWeights(const triPointRef &, const point &, FixedList< scalar, 3 > &weights)
Calculate linear interpolation weights for point (guaranteed to be.
static label getTriangle(const triSurface &surf, const label e0I, const label e1I, const label e2I)
Return index of triangle (or -1) using all three edges.
const Point & a() const
Return first vertex.
Definition: triangleI.H:70
const Field< PointType > & pointNormals() const
Return point normals for patch.
const vectorField & faceCentres() const
void setSize(const label)
Reset size of List.
Definition: List.C:281
const geometricSurfacePatchList & patches() const
Definition: triSurface.H:322
int dtris2(int point_num, double point_xy[], int *tri_num, int tri_vert[], int tri_nabe[])
Definition: geompack.C:949
const Point & b() const
Return second vertex.
Definition: triangleI.H:76
label patchi
std::remove_reference< ::Foam::List< labelledTri > >::type::value_type FaceType
static sideType surfaceSide(const triSurface &surf, const point &sample, const label nearestFacei)
Given nearest point (to sample) on surface determines which side.
vector point
Point is a vector.
Definition: point.H:41
Triangulation of three-dimensional polygons.
label newPointi
Definition: readKivaGrid.H:501
label end() const
Return end vertex label.
Definition: edgeI.H:92
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
static void otherEdges(const triSurface &surf, const label facei, const label edgeI, label &e1, label &e2)
Get the two edges on facei counterclockwise after edgeI.
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
A List with indirect addressing.
Definition: fvMatrix.H:106
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:309
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.
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
const labelListList & faceEdges() const
Return face-edge addressing.
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
tmp< pointField > allPoints(const Triangulation &t)
Extract all points in vertex-index order.
void setSize(const label)
Dummy setSize function.
Definition: FixedListI.H:183
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
static triSurface triangulate(const polyBoundaryMesh &mBesh, const labelHashSet &includePatches, const bool verbose=false)
Simple triangulation of (selected patches of) boundaryMesh. Needs.
scalar distance() const
Return distance to hit.
Definition: PointHit.H:139
Contains information about location on a triSurface.
List< Key > toc() const
Return the table of contents.
Definition: HashTable.C:202
static void writeOBJ(const fileName &fName, const pointField &pts)
Write pointField to OBJ format file.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
Triangulated surface description with patch information.
Definition: triSurface.H:66
static const label ANYEDGE
Face collapse status.
static label getEdge(const triSurface &surf, const label vert1I, const label vert2I)
Returns edge label connecting v1, v2 (local numbering)
static void getVertexTriangles(const triSurface &surf, const label edgeI, labelList &edgeTris)
Get all triangles using edge endpoint.
static surfaceLocation classify(const triSurface &, const label triI, const point &trianglePoint)
Test point on plane of triangle to see if on edge or point or inside.
label index() const
Return index.
static label oppositeEdge(const triSurface &surf, const label facei, const label vertI)
Get edge opposite vertex (local numbering)
label start() const
Return start vertex label.
Definition: edgeI.H:81
bool containsAny(const UList< point > &) const
Contains any of the points? (inside or on edge)
Definition: boundBox.C:261