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