faceTriangulation.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "faceTriangulation.H"
27 #include "plane.H"
28 
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 const Foam::scalar Foam::faceTriangulation::edgeRelTol = 1e-6;
33 
34 
35 // Edge to the right of face vertex i
36 Foam::label Foam::faceTriangulation::right(const label, label i)
37 {
38  return i;
39 }
40 
41 
42 // Edge to the left of face vertex i
43 Foam::label Foam::faceTriangulation::left(const label size, label i)
44 {
45  return i ? i-1 : size-1;
46 }
47 
48 
49 // Calculate (normalized) edge vectors.
50 // edges[i] gives edge between point i+1 and i.
51 Foam::tmp<Foam::vectorField> Foam::faceTriangulation::calcEdges
52 (
53  const face& f,
54  const pointField& points
55 )
56 {
57  tmp<vectorField> tedges(new vectorField(f.size()));
58  vectorField& edges = tedges.ref();
59 
60  forAll(f, i)
61  {
62  point thisPt = points[f[i]];
63  point nextPt = points[f[f.fcIndex(i)]];
64 
65  vector vec(nextPt - thisPt);
66  vec /= mag(vec) + VSMALL;
67 
68  edges[i] = vec;
69  }
70 
71  return tedges;
72 }
73 
74 
75 // Calculates half angle components of angle from e0 to e1
76 void Foam::faceTriangulation::calcHalfAngle
77 (
78  const vector& normal,
79  const vector& e0,
80  const vector& e1,
81  scalar& cosHalfAngle,
82  scalar& sinHalfAngle
83 )
84 {
85  // truncate cos to +-1 to prevent negative numbers
86  scalar cos = max(-1, min(1, e0 & e1));
87 
88  scalar sin = (e0 ^ e1) & normal;
89 
90  if (sin < -ROOTVSMALL)
91  {
92  // 3rd or 4th quadrant
93  cosHalfAngle = - Foam::sqrt(0.5*(1 + cos));
94  sinHalfAngle = Foam::sqrt(0.5*(1 - cos));
95  }
96  else
97  {
98  // 1st or 2nd quadrant
99  cosHalfAngle = Foam::sqrt(0.5*(1 + cos));
100  sinHalfAngle = Foam::sqrt(0.5*(1 - cos));
101  }
102 }
103 
104 
105 // Calculate intersection point between edge p1-p2 and ray (in 2D).
106 // Return true and intersection point if intersection between p1 and p2.
107 Foam::pointHit Foam::faceTriangulation::rayEdgeIntersect
108 (
109  const vector& normal,
110  const point& rayOrigin,
111  const vector& rayDir,
112  const point& p1,
113  const point& p2,
114  scalar& posOnEdge
115 )
116 {
117  // Start off from miss
118  pointHit result(p1);
119 
120  // Construct plane normal to rayDir and intersect
121  const vector y = normal ^ rayDir;
122 
123  posOnEdge = plane(rayOrigin, y).normalIntersect(p1, (p2-p1));
124 
125  // Check intersection to left of p1 or right of p2
126  if ((posOnEdge < 0) || (posOnEdge > 1))
127  {
128  // Miss
129  }
130  else
131  {
132  // Check intersection behind rayOrigin
133  point intersectPt = p1 + posOnEdge * (p2 - p1);
134 
135  if (((intersectPt - rayOrigin) & rayDir) < 0)
136  {
137  // Miss
138  }
139  else
140  {
141  // Hit
142  result.setHit();
143  result.setPoint(intersectPt);
144  result.setDistance(mag(intersectPt - rayOrigin));
145  }
146  }
147  return result;
148 }
149 
150 
151 // Return true if triangle given its three points (anticlockwise ordered)
152 // contains point
153 bool Foam::faceTriangulation::triangleContainsPoint
154 (
155  const vector& n,
156  const point& p0,
157  const point& p1,
158  const point& p2,
159  const point& pt
160 )
161 {
162  scalar area01Pt = triPointRef(p0, p1, pt).normal() & n;
163  scalar area12Pt = triPointRef(p1, p2, pt).normal() & n;
164  scalar area20Pt = triPointRef(p2, p0, pt).normal() & n;
165 
166  if ((area01Pt > 0) && (area12Pt > 0) && (area20Pt > 0))
167  {
168  return true;
169  }
170  else if ((area01Pt < 0) && (area12Pt < 0) && (area20Pt < 0))
171  {
173  return false;
174  }
175  else
176  {
177  return false;
178  }
179 }
180 
181 
182 // Starting from startIndex find diagonal. Return in index1, index2.
183 // Index1 always startIndex except when convex polygon
184 void Foam::faceTriangulation::findDiagonal
185 (
186  const pointField& points,
187  const face& f,
188  const vectorField& edges,
189  const vector& normal,
190  const label startIndex,
191  label& index1,
192  label& index2
193 )
194 {
195  const point& startPt = points[f[startIndex]];
196 
197  // Calculate angle at startIndex
198  const vector& rightE = edges[right(f.size(), startIndex)];
199  const vector leftE = -edges[left(f.size(), startIndex)];
200 
201  // Construct ray which bisects angle
202  scalar cosHalfAngle = GREAT;
203  scalar sinHalfAngle = GREAT;
204  calcHalfAngle(normal, rightE, leftE, cosHalfAngle, sinHalfAngle);
205  vector rayDir
206  (
207  cosHalfAngle*rightE
208  + sinHalfAngle*(normal ^ rightE)
209  );
210  // rayDir should be normalized already but is not due to rounding errors
211  // so normalize.
212  rayDir /= mag(rayDir) + VSMALL;
213 
214 
215  //
216  // Check all edges (apart from rightE and leftE) for nearest intersection
217  //
218 
219  label faceVertI = f.fcIndex(startIndex);
220 
221  pointHit minInter(false, Zero, GREAT, true);
222  label minIndex = -1;
223  scalar minPosOnEdge = GREAT;
224 
225  for (label i = 0; i < f.size() - 2; i++)
226  {
227  scalar posOnEdge;
228  pointHit inter =
229  rayEdgeIntersect
230  (
231  normal,
232  startPt,
233  rayDir,
234  points[f[faceVertI]],
235  points[f[f.fcIndex(faceVertI)]],
236  posOnEdge
237  );
238 
239  if (inter.hit() && inter.distance() < minInter.distance())
240  {
241  minInter = inter;
242  minIndex = faceVertI;
243  minPosOnEdge = posOnEdge;
244  }
245 
246  faceVertI = f.fcIndex(faceVertI);
247  }
248 
249 
250  if (minIndex == -1)
251  {
252  //WarningInFunction
253  // << "Could not find intersection starting from " << f[startIndex]
254  // << " for face " << f << endl;
255 
256  index1 = -1;
257  index2 = -1;
258  return;
259  }
260 
261  const label leftIndex = minIndex;
262  const label rightIndex = f.fcIndex(minIndex);
263 
264  // Now ray intersects edge from leftIndex to rightIndex.
265  // Check for intersection being one of the edge points. Make sure never
266  // to return two consecutive points.
267 
268  if (mag(minPosOnEdge) < edgeRelTol && f.fcIndex(startIndex) != leftIndex)
269  {
270  index1 = startIndex;
271  index2 = leftIndex;
272 
273  return;
274  }
275  if
276  (
277  mag(minPosOnEdge - 1) < edgeRelTol
278  && f.fcIndex(rightIndex) != startIndex
279  )
280  {
281  index1 = startIndex;
282  index2 = rightIndex;
283 
284  return;
285  }
286 
287  // Select visible vertex that minimizes
288  // angle to bisection. Visibility checking by checking if inside triangle
289  // formed by startIndex, leftIndex, rightIndex
290 
291  const point& leftPt = points[f[leftIndex]];
292  const point& rightPt = points[f[rightIndex]];
293 
294  minIndex = -1;
295  scalar maxCos = -GREAT;
296 
297  // all vertices except for startIndex and ones to left and right of it.
298  faceVertI = f.fcIndex(f.fcIndex(startIndex));
299  for (label i = 0; i < f.size() - 3; i++)
300  {
301  const point& pt = points[f[faceVertI]];
302 
303  if
304  (
305  (faceVertI == leftIndex)
306  || (faceVertI == rightIndex)
307  || (triangleContainsPoint(normal, startPt, leftPt, rightPt, pt))
308  )
309  {
310  // pt inside triangle (so perhaps visible)
311  // Select based on minimal angle (so guaranteed visible).
312  vector edgePt0 = pt - startPt;
313  edgePt0 /= mag(edgePt0);
314 
315  scalar cos = rayDir & edgePt0;
316  if (cos > maxCos)
317  {
318  maxCos = cos;
319  minIndex = faceVertI;
320  }
321  }
322  faceVertI = f.fcIndex(faceVertI);
323  }
324 
325  if (minIndex == -1)
326  {
327  // no vertex found. Return startIndex and one of the intersected edge
328  // endpoints.
329  index1 = startIndex;
330 
331  if (f.fcIndex(startIndex) != leftIndex)
332  {
333  index2 = leftIndex;
334  }
335  else
336  {
337  index2 = rightIndex;
338  }
339 
340  return;
341  }
342 
343  index1 = startIndex;
344  index2 = minIndex;
345 }
346 
347 
348 // Find label of vertex to start splitting from. Is:
349 // 1] flattest concave angle
350 // 2] flattest convex angle if no concave angles.
351 Foam::label Foam::faceTriangulation::findStart
352 (
353  const face& f,
354  const vectorField& edges,
355  const vector& normal
356 )
357 {
358  const label size = f.size();
359 
360  scalar minCos = GREAT;
361  label minIndex = -1;
362 
363  forAll(f, fp)
364  {
365  const vector& rightEdge = edges[right(size, fp)];
366  const vector leftEdge = -edges[left(size, fp)];
367 
368  if (((rightEdge ^ leftEdge) & normal) < ROOTVSMALL)
369  {
370  scalar cos = rightEdge & leftEdge;
371  if (cos < minCos)
372  {
373  minCos = cos;
374  minIndex = fp;
375  }
376  }
377  }
378 
379  if (minIndex == -1)
380  {
381  // No concave angle found. Get flattest convex angle
382  minCos = GREAT;
383 
384  forAll(f, fp)
385  {
386  const vector& rightEdge = edges[right(size, fp)];
387  const vector leftEdge = -edges[left(size, fp)];
388 
389  scalar cos = rightEdge & leftEdge;
390  if (cos < minCos)
391  {
392  minCos = cos;
393  minIndex = fp;
394  }
395  }
396  }
397 
398  return minIndex;
399 }
400 
401 
402 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
403 
404 // Split face f into triangles. Handles all simple (convex & concave)
405 // polygons.
406 bool Foam::faceTriangulation::split
407 (
408  const bool fallBack,
409  const pointField& points,
410  const face& f,
411  const vector& normal,
412  label& triI
413 )
414 {
415  const label size = f.size();
416 
417  if (size <= 2)
418  {
420  << "Illegal face:" << f
421  << " with points " << UIndirectList<point>(points, f)()
422  << endl;
423 
424  return false;
425  }
426  else if (size == 3)
427  {
428  // Triangle. Just copy.
429  triFace& tri = operator[](triI++);
430  tri[0] = f[0];
431  tri[1] = f[1];
432  tri[2] = f[2];
433 
434  return true;
435  }
436  else
437  {
438  // General case. Start splitting for -flattest concave angle
439  // -or flattest convex angle if no concave angles.
440 
441  tmp<vectorField> tedges(calcEdges(f, points));
442  const vectorField& edges = tedges();
443 
444  label startIndex = findStart(f, edges, normal);
445 
446  // Find diagonal to split face across
447  label index1 = -1;
448  label index2 = -1;
449 
450  forAll(f, iter)
451  {
452  findDiagonal
453  (
454  points,
455  f,
456  edges,
457  normal,
458  startIndex,
459  index1,
460  index2
461  );
462 
463  if (index1 != -1 && index2 != -1)
464  {
465  // Found correct diagonal
466  break;
467  }
468 
469  // Try splitting from next startingIndex.
470  startIndex = f.fcIndex(startIndex);
471  }
472 
473  if (index1 == -1 || index2 == -1)
474  {
475  if (fallBack)
476  {
477  // Do naive triangulation. Find smallest angle to start
478  // triangulating from.
479  label maxIndex = -1;
480  scalar maxCos = -GREAT;
481 
482  forAll(f, fp)
483  {
484  const vector& rightEdge = edges[right(size, fp)];
485  const vector leftEdge = -edges[left(size, fp)];
486 
487  scalar cos = rightEdge & leftEdge;
488  if (cos > maxCos)
489  {
490  maxCos = cos;
491  maxIndex = fp;
492  }
493  }
494 
496  << "Cannot find valid diagonal on face " << f
497  << " with points " << UIndirectList<point>(points, f)()
498  << nl
499  << "Returning naive triangulation starting from "
500  << f[maxIndex] << " which might not be correct for a"
501  << " concave or warped face" << endl;
502 
503 
504  label fp = f.fcIndex(maxIndex);
505 
506  for (label i = 0; i < size-2; i++)
507  {
508  label nextFp = f.fcIndex(fp);
509 
510  triFace& tri = operator[](triI++);
511  tri[0] = f[maxIndex];
512  tri[1] = f[fp];
513  tri[2] = f[nextFp];
514 
515  fp = nextFp;
516  }
517 
518  return true;
519  }
520  else
521  {
523  << "Cannot find valid diagonal on face " << f
524  << " with points " << UIndirectList<point>(points, f)()
525  << nl
526  << "Returning empty triFaceList" << endl;
527 
528  return false;
529  }
530  }
531 
532 
533  // Split into two subshapes.
534  // face1: index1 to index2
535  // face2: index2 to index1
536 
537  // Get sizes of the two subshapes
538  label diff = 0;
539  if (index2 > index1)
540  {
541  diff = index2 - index1;
542  }
543  else
544  {
545  // folded round
546  diff = index2 + size - index1;
547  }
548 
549  label nPoints1 = diff + 1;
550  label nPoints2 = size - diff + 1;
551 
552  if (nPoints1 == size || nPoints2 == size)
553  {
555  << "Illegal split of face:" << f
556  << " with points " << UIndirectList<point>(points, f)()
557  << " at indices " << index1 << " and " << index2
558  << abort(FatalError);
559  }
560 
561 
562  // Collect face1 points
563  face face1(nPoints1);
564 
565  label faceVertI = index1;
566  for (int i = 0; i < nPoints1; i++)
567  {
568  face1[i] = f[faceVertI];
569  faceVertI = f.fcIndex(faceVertI);
570  }
571 
572  // Collect face2 points
573  face face2(nPoints2);
574 
575  faceVertI = index2;
576  for (int i = 0; i < nPoints2; i++)
577  {
578  face2[i] = f[faceVertI];
579  faceVertI = f.fcIndex(faceVertI);
580  }
581 
582  // Decompose the split faces
583  //Pout<< "Split face:" << f << " into " << face1 << " and " << face2
584  // << endl;
585  //string oldPrefix(Pout.prefix());
586  //Pout.prefix() = " " + oldPrefix;
587 
588  bool splitOk =
589  split(fallBack, points, face1, normal, triI)
590  && split(fallBack, points, face2, normal, triI);
591 
592  //Pout.prefix() = oldPrefix;
593 
594  return splitOk;
595  }
596 }
597 
598 
599 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
600 
601 // Null constructor
603 :
604  triFaceList()
605 {}
606 
607 
608 // Construct from components
610 (
611  const pointField& points,
612  const face& f,
613  const bool fallBack
614 )
615 :
616  triFaceList(f.size()-2)
617 {
618  vector avgNormal = f.normal(points);
619  avgNormal /= mag(avgNormal) + VSMALL;
620 
621  label triI = 0;
622 
623  bool valid = split(fallBack, points, f, avgNormal, triI);
624 
625  if (!valid)
626  {
627  setSize(0);
628  }
629 }
630 
631 
632 // Construct from components
634 (
635  const pointField& points,
636  const face& f,
637  const vector& n,
638  const bool fallBack
639 )
640 :
641  triFaceList(f.size()-2)
642 {
643  label triI = 0;
644 
645  bool valid = split(fallBack, points, f, n, triI);
646 
647  if (!valid)
648  {
649  setSize(0);
650  }
651 }
652 
653 
654 // Construct from Istream
656 :
657  triFaceList(is)
658 {}
659 
660 
661 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
T & operator[](const label)
Return element of UList.
Definition: UListI.H:167
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
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
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
Definition: Istream.H:57
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
vector normal(const pointField &) const
Vector normal; magnitude is equal to area of face.
Definition: face.C:552
faceTriangulation()
Construct null.
List< triFace > triFaceList
list of triFaces
Definition: triFaceList.H:42
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
const dimensionedScalar e
Elementary charge.
Definition: doubleFloat.H:78
dimensionedScalar cos(const dimensionedScalar &ds)
const pointField & points
face triFace(3)
static const zero Zero
Definition: zero.H:91
errorManip< error > abort(error &err)
Definition: errorManip.H:131
static const char nl
Definition: Ostream.H:262
labelList f(nPoints)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void setSize(const label)
Reset size of List.
Definition: List.C:295
vector point
Point is a vector.
Definition: point.H:41
#define WarningInFunction
Report a warning using Foam::Warning.
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
Field< vector > vectorField
Specialisation of Field<T> for vector.
PointHit< point > pointHit
Definition: pointHit.H:41
A class for managing temporary objects.
Definition: PtrList.H:54