triSurfaceTools.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2016 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 
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 responsability 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 dont 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 
2331 // triangulation of boundaryMesh
2334  const polyBoundaryMesh& bMesh,
2335  const labelHashSet& includePatches,
2336  const bool verbose
2337 )
2338 {
2339  const polyMesh& mesh = bMesh.mesh();
2340 
2341  // Storage for surfaceMesh. Size estimate.
2342  DynamicList<labelledTri> triangles
2343  (
2344  mesh.nFaces() - mesh.nInternalFaces()
2345  );
2346 
2347  label newPatchi = 0;
2348 
2349  forAllConstIter(labelHashSet, includePatches, iter)
2350  {
2351  const label patchi = iter.key();
2352  const polyPatch& patch = bMesh[patchi];
2353  const pointField& points = patch.points();
2354 
2355  label nTriTotal = 0;
2356 
2357  forAll(patch, patchFacei)
2358  {
2359  const face& f = patch[patchFacei];
2360 
2361  faceList triFaces(f.nTriangles(points));
2362 
2363  label nTri = 0;
2364 
2365  f.triangles(points, nTri, triFaces);
2366 
2367  forAll(triFaces, triFacei)
2368  {
2369  const face& f = triFaces[triFacei];
2370 
2371  triangles.append(labelledTri(f[0], f[1], f[2], newPatchi));
2372 
2373  nTriTotal++;
2374  }
2375  }
2376 
2377  if (verbose)
2378  {
2379  Pout<< patch.name() << " : generated " << nTriTotal
2380  << " triangles from " << patch.size() << " faces with"
2381  << " new patchid " << newPatchi << endl;
2382  }
2383 
2384  newPatchi++;
2385  }
2386  triangles.shrink();
2387 
2388  // Create globally numbered tri surface
2389  triSurface rawSurface(triangles, mesh.points());
2390 
2391  // Create locally numbered tri surface
2392  triSurface surface
2393  (
2394  rawSurface.localFaces(),
2395  rawSurface.localPoints()
2396  );
2397 
2398  // Add patch names to surface
2399  surface.patches().setSize(newPatchi);
2400 
2401  newPatchi = 0;
2402 
2403  forAllConstIter(labelHashSet, includePatches, iter)
2404  {
2405  const label patchi = iter.key();
2406  const polyPatch& patch = bMesh[patchi];
2407 
2408  surface.patches()[newPatchi].name() = patch.name();
2409  surface.patches()[newPatchi].geometricType() = patch.type();
2410 
2411  newPatchi++;
2412  }
2413 
2414  return surface;
2415 }
2416 
2417 
2418 // triangulation of boundaryMesh
2421  const polyBoundaryMesh& bMesh,
2422  const labelHashSet& includePatches,
2423  const bool verbose
2424 )
2425 {
2426  const polyMesh& mesh = bMesh.mesh();
2427 
2428  // Storage for new points = meshpoints + face centres.
2429  const pointField& points = mesh.points();
2430  const pointField& faceCentres = mesh.faceCentres();
2431 
2432  pointField newPoints(points.size() + faceCentres.size());
2433 
2434  label newPointi = 0;
2435 
2436  forAll(points, pointi)
2437  {
2438  newPoints[newPointi++] = points[pointi];
2439  }
2440  forAll(faceCentres, facei)
2441  {
2442  newPoints[newPointi++] = faceCentres[facei];
2443  }
2444 
2445 
2446  // Count number of faces.
2447  DynamicList<labelledTri> triangles
2448  (
2449  mesh.nFaces() - mesh.nInternalFaces()
2450  );
2451 
2452  label newPatchi = 0;
2453 
2454  forAllConstIter(labelHashSet, includePatches, iter)
2455  {
2456  const label patchi = iter.key();
2457  const polyPatch& patch = bMesh[patchi];
2458 
2459  label nTriTotal = 0;
2460 
2461  forAll(patch, patchFacei)
2462  {
2463  // Face in global coords.
2464  const face& f = patch[patchFacei];
2465 
2466  // Index in newPointi of face centre.
2467  label fc = points.size() + patchFacei + patch.start();
2468 
2469  forAll(f, fp)
2470  {
2471  label fp1 = f.fcIndex(fp);
2472 
2473  triangles.append(labelledTri(f[fp], f[fp1], fc, newPatchi));
2474 
2475  nTriTotal++;
2476  }
2477  }
2478 
2479  if (verbose)
2480  {
2481  Pout<< patch.name() << " : generated " << nTriTotal
2482  << " triangles from " << patch.size() << " faces with"
2483  << " new patchid " << newPatchi << endl;
2484  }
2485 
2486  newPatchi++;
2487  }
2488  triangles.shrink();
2489 
2490 
2491  // Create globally numbered tri surface
2492  triSurface rawSurface(triangles, newPoints);
2493 
2494  // Create locally numbered tri surface
2495  triSurface surface
2496  (
2497  rawSurface.localFaces(),
2498  rawSurface.localPoints()
2499  );
2500 
2501  // Add patch names to surface
2502  surface.patches().setSize(newPatchi);
2503 
2504  newPatchi = 0;
2505 
2506  forAllConstIter(labelHashSet, includePatches, iter)
2507  {
2508  const label patchi = iter.key();
2509  const polyPatch& patch = bMesh[patchi];
2510 
2511  surface.patches()[newPatchi].name() = patch.name();
2512  surface.patches()[newPatchi].geometricType() = patch.type();
2513 
2514  newPatchi++;
2515  }
2516 
2517  return surface;
2518 }
2519 
2520 
2522 {
2523  // Vertices in geompack notation. Note that could probably just use
2524  // pts.begin() if double precision.
2525  List<doubleScalar> geompackVertices(2*pts.size());
2526  label doubleI = 0;
2527  forAll(pts, i)
2528  {
2529  geompackVertices[doubleI++] = pts[i][0];
2530  geompackVertices[doubleI++] = pts[i][1];
2531  }
2532 
2533  // Storage for triangles
2534  int m2 = 3;
2535  List<int> triangle_node(m2*3*pts.size());
2536  List<int> triangle_neighbor(m2*3*pts.size());
2537 
2538  // Triangulate
2539  int nTris = 0;
2540  int err = dtris2
2541  (
2542  pts.size(),
2543  geompackVertices.begin(),
2544  &nTris,
2545  triangle_node.begin(),
2546  triangle_neighbor.begin()
2547  );
2548 
2549  if (err != 0)
2550  {
2552  << "Failed dtris2 with vertices:" << pts.size()
2553  << abort(FatalError);
2554  }
2555 
2556  // Trim
2557  triangle_node.setSize(3*nTris);
2558  triangle_neighbor.setSize(3*nTris);
2559 
2560  // Convert to triSurface.
2561  List<labelledTri> faces(nTris);
2562 
2563  forAll(faces, i)
2564  {
2565  faces[i] = labelledTri
2566  (
2567  triangle_node[3*i]-1,
2568  triangle_node[3*i+1]-1,
2569  triangle_node[3*i+2]-1,
2570  0
2571  );
2572  }
2573 
2574  pointField points(pts.size());
2575  forAll(pts, i)
2576  {
2577  points[i][0] = pts[i][0];
2578  points[i][1] = pts[i][1];
2579  points[i][2] = 0.0;
2580  }
2581 
2582  return triSurface(faces, points);
2583 }
2584 
2585 
2588  const triPointRef& tri,
2589  const point& p,
2590  FixedList<scalar, 3>& weights
2591 )
2592 {
2593  // calculate triangle edge vectors and triangle face normal
2594  // the 'i':th edge is opposite node i
2596  edge[0] = tri.c()-tri.b();
2597  edge[1] = tri.a()-tri.c();
2598  edge[2] = tri.b()-tri.a();
2599 
2600  vector triangleFaceNormal = edge[1] ^ edge[2];
2601 
2602  // calculate edge normal (pointing inwards)
2604  for (label i=0; i<3; i++)
2605  {
2606  normal[i] = triangleFaceNormal ^ edge[i];
2607  normal[i] /= mag(normal[i]) + VSMALL;
2608  }
2609 
2610  weights[0] = ((p-tri.b()) & normal[0]) / max(VSMALL, normal[0] & edge[1]);
2611  weights[1] = ((p-tri.c()) & normal[1]) / max(VSMALL, normal[1] & edge[2]);
2612  weights[2] = ((p-tri.a()) & normal[2]) / max(VSMALL, normal[2] & edge[0]);
2613 }
2614 
2615 
2616 // Calculate weighting factors from samplePts to triangle it is in.
2617 // Uses linear search.
2620  const triSurface& s,
2621  const pointField& samplePts,
2622  List<FixedList<label, 3>>& allVerts,
2623  List<FixedList<scalar, 3>>& allWeights
2624 )
2625 {
2626  allVerts.setSize(samplePts.size());
2627  allWeights.setSize(samplePts.size());
2628 
2629  const pointField& points = s.points();
2630 
2631  forAll(samplePts, i)
2632  {
2633  const point& samplePt = samplePts[i];
2634 
2635 
2636  FixedList<label, 3>& verts = allVerts[i];
2637  FixedList<scalar, 3>& weights = allWeights[i];
2638 
2639  scalar minDistance = GREAT;
2640 
2641  forAll(s, facei)
2642  {
2643  const labelledTri& f = s[facei];
2644 
2645  triPointRef tri(f.tri(points));
2646 
2647  label nearType, nearLabel;
2648 
2649  pointHit nearest = tri.nearestPointClassify
2650  (
2651  samplePt,
2652  nearType,
2653  nearLabel
2654  );
2655 
2656  if (nearest.hit())
2657  {
2658  // samplePt inside triangle
2659  verts[0] = f[0];
2660  verts[1] = f[1];
2661  verts[2] = f[2];
2662 
2663  calcInterpolationWeights(tri, nearest.rawPoint(), weights);
2664 
2665  //Pout<< "calcScalingFactors : samplePt:" << samplePt
2666  // << " inside triangle:" << facei
2667  // << " verts:" << verts
2668  // << " weights:" << weights
2669  // << endl;
2670 
2671  break;
2672  }
2673  else if (nearest.distance() < minDistance)
2674  {
2675  minDistance = nearest.distance();
2676 
2677  // Outside triangle. Store nearest.
2678 
2679  if (nearType == triPointRef::POINT)
2680  {
2681  verts[0] = f[nearLabel];
2682  weights[0] = 1;
2683  verts[1] = -1;
2684  weights[1] = -GREAT;
2685  verts[2] = -1;
2686  weights[2] = -GREAT;
2687 
2688  //Pout<< "calcScalingFactors : samplePt:" << samplePt
2689  // << " distance:" << nearest.distance()
2690  // << " from point:" << points[f[nearLabel]]
2691  // << endl;
2692  }
2693  else if (nearType == triPointRef::EDGE)
2694  {
2695  verts[0] = f[nearLabel];
2696  verts[1] = f[f.fcIndex(nearLabel)];
2697  verts[2] = -1;
2698 
2699  const point& p0 = points[verts[0]];
2700  const point& p1 = points[verts[1]];
2701 
2702  scalar s = min
2703  (
2704  1,
2705  max
2706  (
2707  0,
2708  mag(nearest.rawPoint() - p0)/mag(p1 - p0)
2709  )
2710  );
2711 
2712  // Interpolate
2713  weights[0] = 1 - s;
2714  weights[1] = s;
2715  weights[2] = -GREAT;
2716 
2717  //Pout<< "calcScalingFactors : samplePt:" << samplePt
2718  // << " distance:" << nearest.distance()
2719  // << " from edge:" << p0 << p1 << " s:" << s
2720  // << endl;
2721  }
2722  else
2723  {
2724  // triangle. Can only happen because of truncation errors.
2725  verts[0] = f[0];
2726  verts[1] = f[1];
2727  verts[2] = f[2];
2728 
2729  calcInterpolationWeights(tri, nearest.rawPoint(), weights);
2730 
2731  //Pout<< "calcScalingFactors : samplePt:" << samplePt
2732  // << " distance:" << nearest.distance()
2733  // << " to verts:" << verts
2734  // << " weights:" << weights
2735  // << endl;
2736 
2737  break;
2738  }
2739  }
2740  }
2741  }
2742 }
2743 
2744 
2745 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
2746 // Tracking:
2747 
2748 
2749 // Test point on surface to see if is on face,edge or point.
2752  const triSurface& s,
2753  const label triI,
2754  const point& trianglePoint
2755 )
2756 {
2757  surfaceLocation nearest;
2758 
2759  // Nearest point could be on point or edge. Retest.
2760  label index, elemType;
2761  //bool inside =
2762  triPointRef(s[triI].tri(s.points())).classify
2763  (
2764  trianglePoint,
2765  elemType,
2766  index
2767  );
2768 
2769  nearest.setPoint(trianglePoint);
2770 
2771  if (elemType == triPointRef::NONE)
2772  {
2773  nearest.setHit();
2774  nearest.setIndex(triI);
2775  nearest.elementType() = triPointRef::NONE;
2776  }
2777  else if (elemType == triPointRef::EDGE)
2778  {
2779  nearest.setMiss();
2780  nearest.setIndex(s.faceEdges()[triI][index]);
2781  nearest.elementType() = triPointRef::EDGE;
2782  }
2783  else //if (elemType == triPointRef::POINT)
2784  {
2785  nearest.setMiss();
2786  nearest.setIndex(s.localFaces()[triI][index]);
2787  nearest.elementType() = triPointRef::POINT;
2788  }
2789 
2790  return nearest;
2791 }
2792 
2793 
2796  const triSurface& s,
2797  const surfaceLocation& start,
2798  const surfaceLocation& end,
2799  const plane& cutPlane
2800 )
2801 {
2802  // Start off from starting point
2803  surfaceLocation nearest = start;
2804  nearest.setMiss();
2805 
2806  // See if in same triangle as endpoint. If so snap.
2807  snapToEnd(s, end, nearest);
2808 
2809  if (!nearest.hit())
2810  {
2811  // Not yet at end point
2812 
2813  if (start.elementType() == triPointRef::NONE)
2814  {
2815  // Start point is inside triangle. Trivial cases already handled
2816  // above.
2817 
2818  // end point is on edge or point so cross currrent triangle to
2819  // see which edge is cut.
2820 
2821  nearest = cutEdge
2822  (
2823  s,
2824  start.index(), // triangle
2825  -1, // excludeEdge
2826  -1, // excludePoint
2827  start.rawPoint(),
2828  cutPlane,
2829  end.rawPoint()
2830  );
2831  nearest.elementType() = triPointRef::EDGE;
2832  nearest.triangle() = start.index();
2833  nearest.setMiss();
2834  }
2835  else if (start.elementType() == triPointRef::EDGE)
2836  {
2837  // Pick connected triangle that is most in direction.
2838  const labelList& eFaces = s.edgeFaces()[start.index()];
2839 
2840  nearest = visitFaces
2841  (
2842  s,
2843  eFaces,
2844  start,
2845  start.index(), // excludeEdgeI
2846  -1, // excludePointi
2847  end,
2848  cutPlane
2849  );
2850  }
2851  else // start.elementType() == triPointRef::POINT
2852  {
2853  const labelList& pFaces = s.pointFaces()[start.index()];
2854 
2855  nearest = visitFaces
2856  (
2857  s,
2858  pFaces,
2859  start,
2860  -1, // excludeEdgeI
2861  start.index(), // excludePointi
2862  end,
2863  cutPlane
2864  );
2865  }
2866  snapToEnd(s, end, nearest);
2867  }
2868  return nearest;
2869 }
2870 
2871 
2874  const triSurface& s,
2875  const surfaceLocation& endInfo,
2876  const plane& cutPlane,
2877  surfaceLocation& hitInfo
2878 )
2879 {
2880  //OFstream str("track.obj");
2881  //label vertI = 0;
2882  //meshTools::writeOBJ(str, hitInfo.rawPoint());
2883  //vertI++;
2884 
2885  // Track across surface.
2886  while (true)
2887  {
2888  //Pout<< "Tracking from:" << nl
2889  // << " " << hitInfo.info()
2890  // << endl;
2891 
2892  hitInfo = trackToEdge
2893  (
2894  s,
2895  hitInfo,
2896  endInfo,
2897  cutPlane
2898  );
2899 
2900  //meshTools::writeOBJ(str, hitInfo.rawPoint());
2901  //vertI++;
2902  //str<< "l " << vertI-1 << ' ' << vertI << nl;
2903 
2904  //Pout<< "Tracked to:" << nl
2905  // << " " << hitInfo.info() << endl;
2906 
2907  if (hitInfo.hit() || hitInfo.triangle() == -1)
2908  {
2909  break;
2910  }
2911  }
2912 }
2913 
2914 
2915 // ************************************************************************* //
pointHit nearestPointClassify(const point &p, const pointField &points, label &nearType, label &nearLabel) const
Return nearest point to face and classify it:
Definition: triFaceI.H:301
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.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:58
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:291
A class for handling file names.
Definition: fileName.H:69
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:319
const Field< PointType > & pointNormals() const
Return point normals for patch.
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:163
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
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:510
static void track(const triSurface &, const surfaceLocation &endInfo, const plane &cutPlane, surfaceLocation &hitInfo)
Track from edge to edge across surface. Uses trackToEdge.
const Cmpt & z() const
Definition: VectorI.H:87
const Field< PointType > & faceNormals() const
Return face normals for patch.
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
edgeList edges() const
Return edges in face point ordering,.
Definition: triFaceI.H:318
label otherVertex(const label a) const
Given one vertex, return the other.
Definition: edgeI.H:103
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
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:116
label triangles(const pointField &points, label &triI, faceList &triFaces) const
Split into triangles using existing points.
Definition: face.C:829
void setPoint(const Point &p)
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:1011
const List< Face > & localFaces() const
Return patch faces addressing into local point list.
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:210
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.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
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:184
iterator begin()
Return an iterator to begin traversing the UList.
Definition: UListI.H:216
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, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
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:91
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
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:135
const Point & rawPoint() const
Return point with no checking.
Foam::polyBoundaryMesh.
const labelListList & pointFaces() const
Return point-face addressing.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
const Field< PointType > & localPoints() const
Return pointField of points in patch.
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:240
sideType
On which side of surface.
static const char nl
Definition: Ostream.H:262
Merge points. See below.
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence 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
A normal distribution model.
const vectorField & faceCentres() const
void setSize(const label)
Reset size of List.
Definition: List.C:281
const geometricSurfacePatchList & patches() const
Definition: triSurface.H:314
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
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
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 dimensionedScalar c
Speed of light in a vacuum.
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
label mergePoints(const UList< Type > &points, const scalar mergeTol, const bool verbose, labelList &pointMap, const Type &origin=Type::zero)
Sorts and merges points. All points closer than/equal mergeTol get merged.
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:193
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.
label nTriangles() const
Number of triangles after splitting.
Definition: faceI.H:130
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:65
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