32 const Foam::scalar Foam::faceTriangulation::edgeRelTol = 1
e-6;
45 return i ? i-1 : size-1;
62 point thisPt = points[f[i]];
63 point nextPt = points[f[f.fcIndex(i)]];
65 vector vec(nextPt - thisPt);
66 vec /=
mag(vec) + vSmall;
76 void Foam::faceTriangulation::calcHalfAngle
86 scalar cos =
max(-1,
min(1, e0 & e1));
88 scalar sin = (e0 ^ e1) & normal;
90 if (sin < -rootVSmall)
110 const point& rayOrigin,
121 const vector y = normal ^ rayDir;
123 posOnEdge = plane(rayOrigin, y).normalIntersect(p1, (p2-p1));
126 if ((posOnEdge < 0) || (posOnEdge > 1))
133 point intersectPt = p1 + posOnEdge * (p2 - p1);
135 if (((intersectPt - rayOrigin) & rayDir) < 0)
143 result.setPoint(intersectPt);
144 result.setDistance(
mag(intersectPt - rayOrigin));
153 bool Foam::faceTriangulation::triangleContainsPoint
166 if ((area01Pt > 0) && (area12Pt > 0) && (area20Pt > 0))
170 else if ((area01Pt < 0) && (area12Pt < 0) && (area20Pt < 0))
184 void Foam::faceTriangulation::findDiagonal
190 const label startIndex,
195 const point& startPt = points[f[startIndex]];
198 const vector& rightE = edges[right(f.size(), startIndex)];
199 const vector leftE = -edges[left(f.size(), startIndex)];
202 scalar cosHalfAngle = great;
203 scalar sinHalfAngle = great;
204 calcHalfAngle(normal, rightE, leftE, cosHalfAngle, sinHalfAngle);
208 + sinHalfAngle*(normal ^ rightE)
212 rayDir /=
mag(rayDir) + vSmall;
219 label faceVertI = f.fcIndex(startIndex);
223 scalar minPosOnEdge = great;
225 for (
label i = 0; i < f.size() - 2; i++)
234 points[f[faceVertI]],
235 points[f[f.fcIndex(faceVertI)]],
239 if (inter.hit() && inter.distance() < minInter.distance())
242 minIndex = faceVertI;
243 minPosOnEdge = posOnEdge;
246 faceVertI = f.fcIndex(faceVertI);
261 const label leftIndex = minIndex;
262 const label rightIndex = f.fcIndex(minIndex);
268 if (
mag(minPosOnEdge) < edgeRelTol && f.fcIndex(startIndex) != leftIndex)
277 mag(minPosOnEdge - 1) < edgeRelTol
278 && f.fcIndex(rightIndex) != startIndex
291 const point& leftPt = points[f[leftIndex]];
292 const point& rightPt = points[f[rightIndex]];
295 scalar maxCos = -great;
298 faceVertI = f.fcIndex(f.fcIndex(startIndex));
299 for (
label i = 0; i < f.size() - 3; i++)
301 const point& pt = points[f[faceVertI]];
305 (faceVertI == leftIndex)
306 || (faceVertI == rightIndex)
307 || (triangleContainsPoint(normal, startPt, leftPt, rightPt, pt))
312 vector edgePt0 = pt - startPt;
313 edgePt0 /=
mag(edgePt0);
315 scalar cos = rayDir & edgePt0;
319 minIndex = faceVertI;
322 faceVertI = f.fcIndex(faceVertI);
331 if (f.fcIndex(startIndex) != leftIndex)
358 const label size = f.size();
360 scalar minCos = great;
365 const vector& rightEdge = edges[right(size, fp)];
366 const vector leftEdge = -edges[left(size, fp)];
368 if (((rightEdge ^ leftEdge) & normal) < rootVSmall)
370 scalar cos = rightEdge & leftEdge;
386 const vector& rightEdge = edges[right(size, fp)];
387 const vector leftEdge = -edges[left(size, fp)];
389 scalar cos = rightEdge & leftEdge;
406 bool Foam::faceTriangulation::split
415 const label size = f.size();
420 <<
"Illegal face:" << f
421 <<
" with points " << UIndirectList<point>(
points,
f)()
441 tmp<vectorField> tedges(calcEdges(f, points));
444 label startIndex = findStart(f, edges, normal);
463 if (index1 != -1 && index2 != -1)
470 startIndex = f.fcIndex(startIndex);
473 if (index1 == -1 || index2 == -1)
480 scalar maxCos = -great;
484 const vector& rightEdge = edges[right(size, fp)];
485 const vector leftEdge = -edges[left(size, fp)];
487 scalar cos = rightEdge & leftEdge;
496 <<
"Cannot find valid diagonal on face " << f
497 <<
" with points " << UIndirectList<point>(
points,
f)()
499 <<
"Returning naive triangulation starting from " 500 << f[maxIndex] <<
" which might not be correct for a" 501 <<
" concave or warped face" <<
endl;
504 label fp = f.fcIndex(maxIndex);
506 for (
label i = 0; i < size-2; i++)
508 label nextFp = f.fcIndex(fp);
511 tri[0] = f[maxIndex];
523 <<
"Cannot find valid diagonal on face " << f
524 <<
" with points " << UIndirectList<point>(
points,
f)()
526 <<
"Returning empty triFaceList" <<
endl;
541 diff = index2 - index1;
546 diff = index2 + size - index1;
549 label nPoints1 = diff + 1;
550 label nPoints2 = size - diff + 1;
552 if (nPoints1 == size || nPoints2 == size)
555 <<
"Illegal split of face:" << f
556 <<
" with points " << UIndirectList<point>(
points,
f)()
557 <<
" at indices " << index1 <<
" and " << index2
563 face face1(nPoints1);
565 label faceVertI = index1;
566 for (
int i = 0; i < nPoints1; i++)
568 face1[i] = f[faceVertI];
569 faceVertI = f.fcIndex(faceVertI);
573 face face2(nPoints2);
576 for (
int i = 0; i < nPoints2; i++)
578 face2[i] = f[faceVertI];
579 faceVertI = f.fcIndex(faceVertI);
589 split(fallBack, points, face1, normal, triI)
590 && split(fallBack, points, face2, normal, triI);
620 bool valid = split(fallBack, points, f, avgNormal, triI);
641 bool valid = split(fallBack, points, f, n, triI);
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
A face is a list of labels corresponding to mesh vertices.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
T & operator[](const label)
Return element of UList.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
void size(const label)
Override size to be inconsistent with allocated storage.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Vector< scalar > vector
A scalar version of the templated Vector.
faceTriangulation()
Construct null.
vector normal(const pointField &) const
Return unit normal.
List< triFace > triFaceList
list of triFaces
vectorField pointField
pointField is a vectorField.
dimensionedScalar cos(const dimensionedScalar &ds)
const dimensionedScalar & e
Elementary charge.
errorManip< error > abort(error &err)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void setSize(const label)
Reset size of List.
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
triangle< point, const point & > triPointRef
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
PointHit< point > pointHit
A class for managing temporary objects.