polygonTriangulate.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2021 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "polygonTriangulate.H"
27 #include "tensor2D.H"
28 
29 // * * * * * * * * * * * Private Static Member Functions * * * * * * * * * * //
30 
31 template<class PointField>
32 Foam::scalar Foam::polygonTriangulate::area
33 (
34  const triFace& triPoints,
35  const PointField& points,
36  const vector& normal
37 )
38 {
39  const point& o = points[triPoints[0]];
40  const vector ao = points[triPoints[1]] - o;
41  const vector ab = points[triPoints[2]] - o;
42 
43  return (ao ^ ab) & normal;
44 }
45 
46 
47 template<class PointField>
48 Foam::scalar Foam::polygonTriangulate::quality
49 (
50  const triFace& triPoints,
51  const PointField& points,
52  const vector& normal
53 )
54 {
55  const scalar An = area(triPoints, points, normal);
56 
57  scalar PSqr = 0;
58  forAll(triPoints, triEdgei)
59  {
60  const point& p0 = points[triPoints[triEdgei]];
61  const point& p1 = points[triPoints[(triEdgei + 1) % 3]];
62  PSqr += magSqr(p0 - p1);
63  }
64 
65  static const scalar equilateralAnByPSqr = sqrt(scalar(3))/12;
66 
67  return PSqr != 0 ? An/PSqr/equilateralAnByPSqr : 0;
68 }
69 
70 
71 template<class PointField>
72 bool Foam::polygonTriangulate::intersection
73 (
74  const edge& edgePointsA,
75  const edge& edgePointsB,
76  const PointField& points,
77  const vector& normal
78 )
79 {
80  const vector tau0 = perpendicular(normal);
81  const vector tau1 = normal ^ tau0;
82 
83  const point& pointA0 = points[edgePointsA[0]];
84  const point& pointA1 = points[edgePointsA[1]];
85  const point& pointB0 = points[edgePointsB[0]];
86  const point& pointB1 = points[edgePointsB[1]];
87 
88  const point2D point2DA0(tau0 & pointA0, tau1 & pointA0);
89  const point2D point2DA1(tau0 & pointA1, tau1 & pointA1);
90  const point2D point2DB0(tau0 & pointB0, tau1 & pointB0);
91  const point2D point2DB1(tau0 & pointB1, tau1 & pointB1);
92 
93  const tensor2D M =
94  tensor2D(point2DA0 - point2DA1, point2DB1 - point2DB0).T();
95 
96  const scalar detM = det(M);
97  const tensor2D detMInvM = cof(M);
98 
99  const vector2D magDetMTu = sign(detM)*detMInvM & (point2DA0 - point2DB0);
100 
101  return 0 <= cmptMin(magDetMTu) && cmptMax(magDetMTu) <= mag(detM);
102 }
103 
104 
105 template<class PointField>
106 Foam::label Foam::polygonTriangulate::nIntersections
107 (
108  const edge& edgePoints,
109  const PointField& points,
110  const vector& normal
111 )
112 {
113  label n = 0;
114 
115  forAll(points, pi)
116  {
117  const edge otherEdgePoints(pi, points.fcIndex(pi));
118 
119  if (edgePoints.commonVertex(otherEdgePoints) == -1)
120  {
121  n += intersection(edgePoints, otherEdgePoints, points, normal);
122  }
123  }
124 
125  return n;
126 }
127 
128 
129 template<class PointField>
130 Foam::scalar Foam::polygonTriangulate::angle
131 (
132  const label pointi,
133  const PointField& points,
134  const vector& normal
135 )
136 {
137  const point& o = points[pointi];
138  const vector oa = points[points.rcIndex(pointi)] - o;
139  const vector ob = points[points.fcIndex(pointi)] - o;
140 
141  const vector oaNegNNOa = oa - normal*(normal & oa);
142  const vector obNegNNOb = ob - normal*(normal & ob);
143 
144  return
145  atan2(normal & (oaNegNNOa ^ obNegNNOb), - oaNegNNOa & obNegNNOb)
147 }
148 
149 
150 template<class PointField>
151 bool Foam::polygonTriangulate::ear
152 (
153  const label pointi,
154  const PointField& points,
155  const vector& normal
156 )
157 {
158  const point& o = points[pointi];
159  const vector oa = points[points.rcIndex(pointi)] - o;
160  const vector ob = points[points.fcIndex(pointi)] - o;
161 
162  const tensor A = tensor(ob, oa, normal).T();
163  const tensor T(A.y() ^ A.z(), A.z() ^ A.x(), A.x() ^ A.y());
164  const scalar detA = det(A);
165 
166  for
167  (
168  label pointj = points.fcIndex(points.fcIndex(pointi));
169  pointj != points.rcIndex(pointi);
170  pointj = points.fcIndex(pointj)
171  )
172  {
173  const vector detAY = (points[pointj] - o) & T;
174 
175  if (detAY.x() > 0 && detAY.y() > 0 && detAY.x() + detAY.y() < detA)
176  {
177  return false;
178  }
179  }
180 
181  return true;
182 }
183 
184 
185 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
186 
188 (
189  Random& rndGen,
190  const label n,
191  const scalar error
192 )
193 {
194  // Get random points on a unit disk
195  List<point> points(n);
196  forAll(points, pointi)
197  {
198  const scalar theta =
199  2*constant::mathematical::pi*rndGen.sample01<scalar>();
200  const scalar r = sqrt(rndGen.sample01<scalar>());
201  points[pointi] = point(r*cos(theta), r*sin(theta), 0);
202  }
203 
204  // Initialise intersected polygon
205  labelList pointis(identity(n));
206 
207  // Reorder until non-intersected
208  bool incomplete = true;
209  while (incomplete)
210  {
211  incomplete = false;
212 
213  for (label piA0 = 0; piA0 < n; ++ piA0)
214  {
215  const label piA1 = points.fcIndex(piA0);
216 
217  for (label piB0 = piA0 + 2; piB0 < (piA0 == 0 ? n - 1 : n); ++ piB0)
218  {
219  const label piB1 = points.fcIndex(piB0);
220 
221  if
222  (
224  (
225  edge(pointis[piA0], pointis[piA1]),
226  edge(pointis[piB0], pointis[piB1]),
227  points,
228  vector(0, 0, 1)
229  )
230  )
231  {
232  SubList<label> subOrder(pointis, piB0 + 1 - piA1, piA1);
233  inplaceReverseList(subOrder);
234  incomplete = true;
235  }
236 
237  if (incomplete) break;
238  }
239 
240  if (incomplete) break;
241  }
242  }
243 
244  // Add noise and potential self-intersection
245  forAll(points, pointi)
246  {
247  const scalar theta =
248  2*constant::mathematical::pi*rndGen.sample01<scalar>();
249  const scalar r = sqrt(rndGen.sample01<scalar>());
250  points[pointi] =
251  (1 - error)*points[pointi]
252  + error*point(r*cos(theta), r*sin(theta), rndGen.sample01<scalar>());
253  }
254 
255  // Reorder and return
256  return List<point>(points, pointis);
257 }
258 
259 
260 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
261 
262 template<class PointField>
263 void Foam::polygonTriangulate::optimiseTriangulation
264 (
265  const label trii,
266  const PointField& points,
267  const vector& normal,
268  UList<triFace>& triPoints,
270  UList<labelPair>& edgeTris
271 )
272 {
273  const triFace& t = triPoints[trii];
274  const scalar q = quality(t, points, normal);
275 
276  forAll(triEdges[trii], triEdgei)
277  {
278  const label edgei = triEdges[trii][triEdgei];
279 
280  const label otherTrii = edgeTris[edgei][edgeTris[edgei][0] == trii];
281 
282  if (otherTrii == -1) continue;
283 
284  const label otherTriEdgei = findIndex(triEdges[otherTrii], edgei);
285 
286  const triFace& otherT = triPoints[otherTrii];
287  const scalar otherQ = quality(otherT, points, normal);
288 
289  const label piA = triPoints[trii][(triEdgei + 1) % 3];
290  const label piB = triPoints[trii][(triEdgei + 2) % 3];
291  const label piC = triPoints[otherTrii][(otherTriEdgei + 1) % 3];
292  const label piD = triPoints[otherTrii][(otherTriEdgei + 2) % 3];
293 
294  const triFace flipT0(piA, piB, piD);
295  const triFace flipT1(piC, piD, piB);
296 
297  if
298  (
299  triPointsSet_.found(flipT0)
300  || triPointsSet_.found(flipT1)
301  ) continue;
302 
303  const scalar flipQ0 = quality(flipT0, points, normal);
304  const scalar flipQ1 = quality(flipT1, points, normal);
305 
306  if (min(q, otherQ) < min(flipQ0, flipQ1))
307  {
308  triPointsSet_.set(t);
309  triPointsSet_.set(otherT);
310 
311  const label eiA = triEdges[trii][(triEdgei + 1) % 3];
312  const label eiB = triEdges[trii][(triEdgei + 2) % 3];
313  const label eiC = triEdges[otherTrii][(otherTriEdgei + 1) % 3];
314  const label eiD = triEdges[otherTrii][(otherTriEdgei + 2) % 3];
315 
316  triPoints[trii] = {piA, piB, piD};
317  triEdges[trii] = {eiA, edgei, eiD};
318  edgeTris[eiD][edgeTris[eiD][0] != otherTrii] = trii;
319 
320  triPoints[otherTrii] = {piC, piD, piB};
321  triEdges[otherTrii] = {eiC, edgei, eiB};
322  edgeTris[eiB][edgeTris[eiB][0] != trii] = otherTrii;
323 
324  optimiseTriangulation
325  (
326  trii,
327  points,
328  normal,
329  triPoints,
330  triEdges,
331  edgeTris
332  );
333  optimiseTriangulation
334  (
335  otherTrii,
336  points,
337  normal,
338  triPoints,
339  triEdges,
340  edgeTris
341  );
342 
343  break;
344  }
345  }
346 }
347 
348 
349 void Foam::polygonTriangulate::simpleTriangulate
350 (
351  const UList<point>& allPoints,
352  const vector& normal,
353  UList<triFace>& triPoints,
355  UList<labelPair>& edgeTris,
356  const bool optimal
357 )
358 {
359  // Get the polygon size
360  const label n = allPoints.size();
361 
362  // Clear the triangulation
363  triPoints = triFace(-1, -1, -1);
364  triEdges = FixedList<label, 3>({-1, -1, -1});
365  edgeTris = labelPair(-1, -1);
366 
367  // Initialise workspace
368  pointis_.setSize(n);
369  edges_.setSize(n);
370  angle_.setSize(n);
371  ear_.setSize(n);
372  forAll(pointis_, i)
373  {
374  pointis_[i] = i;
375  edges_[i] = i;
376  angle_[i] = angle(i, allPoints, normal);
377  ear_[i] = ear(i, allPoints, normal);
378  }
379 
380  // Create the indirect point list for the current un-triangulated polygon
381  UIndirectList<point> points(allPoints, pointis_);
382 
383  // Generate triangles
384  forAll(triPoints, trii)
385  {
386  // Find the ear with the smallest angle
387  scalar minEarAngle = vGreat;
388  label minEarAnglei = -1;
389  for (label i = 0; i < pointis_.size(); ++ i)
390  {
391  if (ear_[i] && minEarAngle > angle_[i])
392  {
393  minEarAngle = angle_[i];
394  minEarAnglei = i;
395  }
396  }
397 
398  // If nothing is an ear, then this face is degenerate. Default to
399  // the smallest angle, ignoring the ear test.
400  if (minEarAnglei == -1)
401  {
402  minEarAnglei = findMin(angle_);
403  }
404 
405  // Get the adjacent points
406  label minEarAnglei0 = pointis_.rcIndex(minEarAnglei);
407  label minEarAnglei1 = pointis_.fcIndex(minEarAnglei);
408 
409  // Add a triangle
410  triPoints[trii] = triFace
411  (
412  pointis_[minEarAnglei0],
413  pointis_[minEarAnglei],
414  pointis_[minEarAnglei1]
415  );
417  ({
418  edges_[minEarAnglei0],
419  edges_[minEarAnglei],
420  trii == triPoints.size() - 1 ? edges_[minEarAnglei1] : n + trii
421  });
422  forAll(triEdges[trii], triEdgei)
423  {
424  const label edgei = triEdges[trii][triEdgei];
425  edgeTris[edgei][edgeTris[edgei][0] != -1] = trii;
426  }
427 
428  // Remove the cut off point from the polygon
429  edges_[minEarAnglei0] = triEdges[trii][2];
430  for (label i = minEarAnglei; i < pointis_.size() - 1; ++ i)
431  {
432  pointis_[i] = pointis_[i + 1];
433  edges_[i] = edges_[i + 1];
434  angle_[i] = angle_[i + 1];
435  ear_[i] = ear_[i + 1];
436  }
437  pointis_.resize(pointis_.size() - 1);
438  edges_.resize(edges_.size() - 1);
439  angle_.resize(angle_.size() - 1);
440  ear_.resize(ear_.size() - 1);
441 
442  // Update connected points
443  if (minEarAnglei0 > minEarAnglei) -- minEarAnglei0;
444  angle_[minEarAnglei0] = angle(minEarAnglei0, points, normal);
445  ear_[minEarAnglei0] = ear(minEarAnglei0, points, normal);
446  if (minEarAnglei1 > minEarAnglei) -- minEarAnglei1;
447  angle_[minEarAnglei1] = angle(minEarAnglei1, points, normal);
448  ear_[minEarAnglei1] = ear(minEarAnglei1, points, normal);
449 
450  // Optimise the quality
451  if (optimal)
452  {
453  triPointsSet_.clear();
454  optimiseTriangulation
455  (
456  trii,
457  allPoints,
458  normal,
459  triPoints,
460  triEdges,
461  edgeTris
462  );
463  }
464  }
465 }
466 
467 
468 void Foam::polygonTriangulate::partitionTriangulate
469 (
470  const UList<point>& points,
471  const vector& normal,
472  const label spanPointi,
473  const label spanEdgei,
474  UList<triFace>& triPoints,
476  UList<labelPair>& edgeTris,
477  const bool simple,
478  const bool optimal
479 )
480 {
481  // Get the polygon size
482  const label n = points.size();
483 
484  // Get the points of the spanning triangle
485  const label piA = (spanEdgei + 1) % n;
486  const label piP = spanPointi;
487  const label piB = spanEdgei;
488 
489  // Determine the sizes of the polygons connected to either side of the
490  // spanning triangle
491  const label nA = (piP - piA + n) % n + 1;
492  const label nB = (piB - piP + n) % n + 1;
493 
494  // Get the edges of the spanning triangle
495  const label eiA = nA >= 3 ? n : piA;
496  const label eiB = nB >= 3 ? n + (nA >= 3) : piP;
497  const label eiE = piB;
498 
499  // Add the spanning triangle
500  triPoints[0] = triFace(piA, piP, piB);
501  triEdges[0] = FixedList<label, 3>({eiA, eiB, eiE});
502 
503  // Rotate the point list so that it can be passed to sub-triangulations
504  // without indirect addressing
505  inplaceRotateList(const_cast<UList<point>&>(points), - piA);
506 
507  // Triangulate the sub-polygon connected to edge A
508  if (nA >= 3)
509  {
510  // Subset the polygon and triangulate
511  SubList<triFace> triPointsA(triPoints, nA - 2, 1);
512  SubList<FixedList<label, 3>> triEdgesA(triEdges, nA - 2, 1);
513  SubList<labelPair> edgeTrisA(edgeTris, 2*nA - 3);
514  triangulate
515  (
516  SubList<point>(points, nA),
517  normal,
518  triPointsA,
519  triEdgesA,
520  edgeTrisA,
521  simple,
522  optimal
523  );
524 
525  // Map the point labels back to the full polygon
526  forAll(triPointsA, triAi)
527  {
528  forAll(triPointsA[triAi], triPointAi)
529  {
530  label& pointi = triPointsA[triAi][triPointAi];
531  pointi = (pointi + piA) % n;
532  }
533  }
534 
535  // Map the edge labels back to the full polygon
536  forAll(triEdgesA, triAi)
537  {
538  forAll(triEdgesA[triAi], triEdgeAi)
539  {
540  label& edgei = triEdgesA[triAi][triEdgeAi];
541  edgei =
542  edgei >= nA ? max(eiA, eiB) + edgei - nA + 1
543  : edgei == nA - 1 ? eiA
544  : (piA + edgei) % n;
545  }
546  }
547 
548  // Map the tri labels back to the full polygon
549  forAll(edgeTrisA, edgeAi)
550  {
551  forAll(edgeTrisA[edgeAi], edgeTriAi)
552  {
553  label& trii = edgeTrisA[edgeAi][edgeTriAi];
554  trii = trii == -1 ? -1 : trii + 1;
555  }
556  }
557  }
558 
559  // Triangulate the sub-polygon connected to edge B
560  if (nB >= 3)
561  {
562  // Subset the polygon and triangulate
563  SubList<triFace> triPointsB(triPoints, nB - 2, nA - 1);
564  SubList<FixedList<label, 3>> triEdgesB(triEdges, nB - 2, nA - 1);
565  SubList<labelPair> edgeTrisB(edgeTris, 2*nB - 3, 2*nA - 3);
566  triangulate
567  (
568  SubList<point>(points, nB, nA - 1),
569  normal,
570  triPointsB,
571  triEdgesB,
572  edgeTrisB,
573  simple,
574  optimal
575  );
576 
577  // Map the point labels back to the full polygon
578  forAll(triPointsB, triBi)
579  {
580  forAll(triPointsB[triBi], triPointBi)
581  {
582  label& pointi = triPointsB[triBi][triPointBi];
583  pointi = (pointi + piA + nA - 1) % n;
584  }
585  }
586 
587  // Map the edge labels back to the full polygon
588  forAll(triEdgesB, triBi)
589  {
590  forAll(triEdgesB[triBi], triEdgeBi)
591  {
592  label& edgei = triEdgesB[triBi][triEdgeBi];
593  edgei =
594  edgei >= nB ? max(eiA, eiB) + nA - 3 + edgei - nB + 1
595  : edgei == nB - 1 ? eiB
596  : (piA + nA - 1 + edgei) % n;
597  }
598  }
599 
600  // Map the tri labels back to the full polygon
601  forAll(edgeTrisB, edgeBi)
602  {
603  forAll(edgeTrisB[edgeBi], edgeTriBi)
604  {
605  label& trii = edgeTrisB[edgeBi][edgeTriBi];
606  trii = trii == -1 ? -1 : trii + nA - 1;
607  }
608  }
609  }
610 
611  // Rotate the point list back
612  inplaceRotateList(const_cast<UList<point>&>(points), piA);
613 
614  // Reorder the edge-tris
615  {
616  // Swap the A-internal-edges and B-outer-edges to get all outer
617  // edges and all internal edges in contiguous blocks
618  SubList<labelPair> l(edgeTris, nA + nB - 3, nA - 1);
619  inplaceRotateList(l, - nA + 2);
620  }
621  {
622  // Rotate the outer edges right by one so that the intersection
623  // edge (currently at the end) moves adjacent to the outer edges
624  SubList<labelPair> l(edgeTris, nA + nB - 3, nA + nB - 2);
625  inplaceRotateList(l, 1);
626  }
627  {
628  // Rotate the outer edges so that they are in sequence with the
629  // original point ordering
630  SubList<labelPair> l(edgeTris, n);
631  inplaceRotateList(l, - n + piA);
632  }
633  if (nA >= 3 && nB >= 3)
634  {
635  // Move edge-B leftwards adjacent to edge-A as this is created
636  // before any of the other internal edges
637  const label fromi = n + nA - 2, toi = n + 1;
638  const labelPair temp = edgeTris[fromi];
639  for (label i = fromi; i > toi; -- i)
640  {
641  edgeTris[i] = edgeTris[i-1];
642  }
643  edgeTris[toi] = temp;
644  }
645 
646  // Add associations to the intersection triangle
647  edgeTris[eiA][edgeTris[eiA][0] != -1] = 0;
648  edgeTris[eiB][edgeTris[eiB][0] != -1] = 0;
649  edgeTris[eiE] = labelPair(0, -1);
650 }
651 
652 
653 void Foam::polygonTriangulate::complexTriangulate
654 (
655  const UList<point>& points,
656  const vector& normal,
657  UList<triFace>& triPoints,
659  UList<labelPair>& edgeTris,
660  const bool optimal
661 )
662 {
663  // Get the polygon size
664  const label n = points.size();
665 
666  // Detect self intersections. When one is found, remove one of the
667  // intersected edges by adding a spanning triangle and recurse. Pick the
668  // spanning triangle that results in the smallest negative area.
669  for (label piA0 = 0; piA0 < n; ++ piA0)
670  {
671  const label piA1 = points.fcIndex(piA0);
672 
673  for (label piB0 = piA0 + 2; piB0 < (piA0 == 0 ? n - 1 : n); ++ piB0)
674  {
675  const label piB1 = points.fcIndex(piB0);
676 
677  if
678  (
680  (
681  edge(piA0, piA1),
682  edge(piB0, piB1),
683  points,
684  normal
685  )
686  )
687  {
688  // Get a bunch of unique spanning triangles that remove one of
689  // the intersected edges
690  FixedList<labelPair, 8> spanPointAndEdgeis
691  ({
692  labelPair(piA0, piB0),
693  labelPair(piA1, piB0),
694  labelPair(piB0, piA0),
695  labelPair(piB1, piA0),
696  labelPair(points.rcIndex(piA0), piA0),
697  labelPair(points.fcIndex(piA1), piA0),
698  labelPair(points.rcIndex(piB0), piB0),
699  labelPair(points.fcIndex(piB1), piB0)
700  });
701  forAll(spanPointAndEdgeis, spani)
702  {
703  for (label spanj = spani + 1; spanj < 8; ++ spanj)
704  {
705  if
706  (
707  spanPointAndEdgeis[spani]
708  == spanPointAndEdgeis[spanj]
709  )
710  {
711  spanPointAndEdgeis[spanj] = labelPair(-1, -1);
712  }
713  }
714  }
715 
716  // Find the spanning triangle that results in the smallest
717  // negative triangulation area
718  scalar spanSumNegA = - vGreat;
719  label spanPointi = -1, spanEdgei = -1;
720  forAll(spanPointAndEdgeis, spani)
721  {
722  if (spanPointAndEdgeis[spani] == labelPair(-1, -1))
723  continue;
724 
725  const label piP = spanPointAndEdgeis[spani].first();
726  const label ei = spanPointAndEdgeis[spani].second();
727  const label piA0 = points.rcIndex(ei);
728  const label piA = ei;
729  const label piB = points.fcIndex(ei);
730  const label piB1 = points.fcIndex(piB);
731 
732  // Only consider this spanning triangle if it decreases
733  // the number of intersections
734  const label netNIntersections =
735  - nIntersections(edge(piA, piB), points, normal)
736  + (piP == piA0 ? -1 : +1)
737  *nIntersections(edge(piP, piA), points, normal)
738  + (piP == piB1 ? -1 : +1)
739  *nIntersections(edge(piP, piB), points, normal);
740  if (netNIntersections >= 0) continue;
741 
742  // Triangulate
743  partitionTriangulate
744  (
745  points,
746  normal,
747  piP,
748  ei,
749  triPoints,
750  triEdges,
751  edgeTris,
752  false,
753  false
754  );
755 
756  // Compute the negative area. If it is smaller than the
757  // previous value then update the spanning triangle. If it
758  // is zero then break and use this triangle.
759  scalar sumNegA = 0;
760  forAll(triPoints, trii)
761  {
762  const scalar a = area(triPoints[trii], points, normal);
763  sumNegA += negPart(a);
764  }
765  if (sumNegA > spanSumNegA)
766  {
767  spanSumNegA = sumNegA;
768  spanPointi = piP;
769  spanEdgei = ei;
770  }
771  if (spanSumNegA == 0)
772  {
773  break;
774  }
775  }
776 
777  // If a suitable spanning triangle was found then partition the
778  // triangulation
779  if (spanPointi != -1)
780  {
781  partitionTriangulate
782  (
783  points,
784  normal,
785  spanPointi,
786  spanEdgei,
787  triPoints,
788  triEdges,
789  edgeTris,
790  false,
791  optimal
792  );
793 
794  return;
795  }
796  }
797  }
798  }
799 
800  // If there are no intersections then do simple polygon triangulation
801  simpleTriangulate(points, normal, triPoints, triEdges, edgeTris, optimal);
802 }
803 
804 
805 void Foam::polygonTriangulate::triangulate
806 (
807  const UList<point>& points,
808  const vector& normal,
809  UList<triFace>& triPoints,
811  UList<labelPair>& edgeTris,
812  const bool simple,
813  const bool optimal
814 )
815 {
816  if (simple)
817  {
818  simpleTriangulate
819  (
820  points,
821  normal,
822  triPoints,
823  triEdges,
824  edgeTris,
825  optimal
826  );
827  }
828  else
829  {
830  complexTriangulate
831  (
832  points,
833  normal,
834  triPoints,
835  triEdges,
836  edgeTris,
837  optimal
838  );
839  }
840 }
841 
842 
843 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
844 
846 :
847  pointis_(),
848  edges_(),
849  angle_(),
850  ear_(),
851 
852  triPointsSet_(),
853 
854  points_(),
855  triPoints_(),
856  triEdges_(),
857  edgeTris_()
858 {}
859 
860 
861 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
862 
864 {}
865 
866 
867 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
868 
869 const Foam::UList<Foam::triFace>& Foam::polygonTriangulate::triangulate
870 (
871  const UList<point>& points,
872  const vector& normal,
873  const bool simple,
874  const bool optimal
875 )
876 {
877  // Get the polygon size
878  const label n = points.size();
879 
880  // Resize workspace
881  triPoints_.resize(n - 2);
882  triEdges_.resize(n - 2);
883  edgeTris_.resize(2*n - 3);
884 
885  // Run
886  triangulate
887  (
888  points,
889  normal,
890  triPoints_,
891  triEdges_,
892  edgeTris_,
893  simple,
894  optimal
895  );
896 
897  return triPoints_;
898 }
899 
900 
901 const Foam::UList<Foam::triFace>& Foam::polygonTriangulate::triangulate
902 (
903  const UList<point>& points,
904  const bool simple,
905  const bool optimal
906 )
907 {
908  return
909  triangulate
910  (
911  points,
912  normalised(face::area(points)),
913  simple,
914  optimal
915  );
916 }
917 
918 
919 // ************************************************************************* //
dimensionedScalar sign(const dimensionedScalar &ds)
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void resize(const label)
Alter the addressed list size.
Definition: DynamicListI.H:216
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
Definition: UListI.H:65
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Tensor< Cmpt > T() const
Return transpose.
Definition: TensorI.H:331
Vector< Cmpt > x() const
Definition: TensorI.H:289
const UList< FixedList< label, 3 > > & triEdges() const
Get the triangles&#39; edges.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
dimensionedScalar sqrt(const dimensionedScalar &ds)
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
label findMin(const ListType &, const label start=0)
Find index of min element (and less than given element).
Vector< Cmpt > y() const
Definition: TensorI.H:296
static vector area(const PointField &ps)
Return vector area given face points.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
dimensionedScalar det(const dimensionedSphericalTensor &dt)
label commonVertex(const edge &a) const
Return common vertex.
Definition: edgeI.H:122
polygonTriangulate()
Construct null.
Templated 2D tensor derived from VectorSpace adding construction from 4 components, element access using xx(), xy(), yx() and yy() member functions and the iner-product (dot-product) and outer-product of two Vector2Ds (tensor-product) operators.
Definition: Tensor2D.H:56
const Cmpt & y() const
Definition: VectorI.H:81
Vector< Cmpt > z() const
Definition: TensorI.H:303
void inplaceReverseList(ListType &list)
Inplace reversal of a list using Swap.
A List obtained as a section of another List.
Definition: SubList.H:53
Foam::intersection.
Definition: intersection.H:49
Type sample01()
Advance the state and return a sample of a given type from a.
Definition: RandomI.H:70
void setSize(const label)
Alter the addressed list size.
Definition: DynamicListI.H:175
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
Form normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:413
dimensionedScalar cos(const dimensionedScalar &ds)
const pointField & points
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:68
face triFace(3)
static List< point > randomPolygon(Random &rndGen, const label n, const scalar error)
Generate a random polygon for testing purposes.
T & first()
Return the first element of the list.
Definition: FixedListI.H:217
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
const Cmpt & x() const
Definition: VectorI.H:75
Random number generator.
Definition: Random.H:57
dimensioned< scalar > magSqr(const dimensioned< Type > &)
dimensionedScalar sin(const dimensionedScalar &ds)
dimensionedScalar atan2(const dimensionedScalar &x, const dimensionedScalar &y)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
Vector< Cmpt > perpendicular(const Vector< Cmpt > &v)
Definition: VectorI.H:166
vector point
Point is a vector.
Definition: point.H:41
void inplaceRotateList(ListType< DataType > &list, label n)
Inplace reversal of a list using the Reversal Block Swapping algorithm.
A List with indirect addressing.
Definition: fvMatrix.H:106
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
void cmptMin(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
label size() const
Return the number of elements in the UList.
Definition: UListI.H:311
Tensor2D< scalar > tensor2D
Tensor2D or scalars.
Definition: tensor2D.H:49
dimensionedSymmTensor cof(const dimensionedSymmTensor &dt)
#define M(I)
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
dimensionedScalar negPart(const dimensionedScalar &ds)