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