34 const Foam::scalar Foam::faceTriangulation::edgeRelTol = 1
e-6;
38 Foam::label Foam::faceTriangulation::right(
const label, label i)
45 Foam::label Foam::faceTriangulation::left(
const label size, label i)
47 return i ? i-1 : size-1;
60 auto& edges = tedges.ref();
67 vector vec(nextPt - thisPt);
69 edges[i] = vec.normalise();
77 void Foam::faceTriangulation::calcHalfAngle
89 scalar
sin = (e0 ^ e1) & normal;
91 if (
sin < -ROOTVSMALL)
111 const point& rayOrigin,
122 const vector y = normal ^ rayDir;
124 posOnEdge = plane(rayOrigin,
y).normalIntersect(p1, (p2-p1));
127 if ((posOnEdge < 0) || (posOnEdge > 1))
134 point intersectPt = p1 + posOnEdge * (p2 - p1);
136 if (((intersectPt - rayOrigin) & rayDir) < 0)
144 result.setPoint(intersectPt);
145 result.setDistance(
mag(intersectPt - rayOrigin));
154 bool Foam::faceTriangulation::triangleContainsPoint
164 const scalar area12Pt =
triPointRef(p1, p2, pt).areaNormal() &
n;
167 if ((area01Pt > 0) && (area12Pt > 0) && (area20Pt > 0))
171 else if ((area01Pt < 0) && (area12Pt < 0) && (area20Pt < 0))
183 void Foam::faceTriangulation::findDiagonal
189 const label startIndex,
197 const vector& rightE = edges[right(
f.size(), startIndex)];
198 const vector leftE = -edges[left(
f.size(), startIndex)];
201 scalar cosHalfAngle = GREAT;
202 scalar sinHalfAngle = GREAT;
203 calcHalfAngle(normal, rightE, leftE, cosHalfAngle, sinHalfAngle);
207 + sinHalfAngle*(normal ^ rightE)
218 label faceVertI =
f.fcIndex(startIndex);
222 scalar minPosOnEdge = GREAT;
224 for (label i = 0; i <
f.size() - 2; i++)
238 if (inter.hit() && inter.distance() < minInter.distance())
241 minIndex = faceVertI;
242 minPosOnEdge = posOnEdge;
245 faceVertI =
f.fcIndex(faceVertI);
260 const label leftIndex = minIndex;
261 const label rightIndex =
f.fcIndex(minIndex);
267 if (
mag(minPosOnEdge) < edgeRelTol &&
f.fcIndex(startIndex) != leftIndex)
276 mag(minPosOnEdge - 1) < edgeRelTol
277 &&
f.fcIndex(rightIndex) != startIndex
294 scalar maxCos = -GREAT;
297 faceVertI =
f.fcIndex(
f.fcIndex(startIndex));
298 for (label i = 0; i <
f.size() - 3; i++)
304 (faceVertI == leftIndex)
305 || (faceVertI == rightIndex)
306 || (triangleContainsPoint(normal, startPt, leftPt, rightPt, pt))
311 vector edgePt0 = (pt - startPt);
314 scalar
cos = rayDir & edgePt0;
318 minIndex = faceVertI;
321 faceVertI =
f.fcIndex(faceVertI);
330 if (
f.fcIndex(startIndex) != leftIndex)
350 Foam::label Foam::faceTriangulation::findStart
357 const label size =
f.size();
359 scalar minCos = GREAT;
364 const vector& rightEdge = edges[right(size, fp)];
365 const vector leftEdge = -edges[left(size, fp)];
367 if (((rightEdge ^ leftEdge) & normal) < ROOTVSMALL)
369 scalar
cos = rightEdge & leftEdge;
385 const vector& rightEdge = edges[right(size, fp)];
386 const vector leftEdge = -edges[left(size, fp)];
388 scalar
cos = rightEdge & leftEdge;
405 bool Foam::faceTriangulation::split
414 const label size =
f.size();
419 <<
"Illegal face:" <<
f
420 <<
" with points " << UIndirectList<point>(
points,
f)
428 triFace& tri = operator[](triI++);
440 tmp<vectorField> tedges(calcEdges(
f,
points));
443 label startIndex = findStart(
f, edges, normal);
462 if (index1 != -1 && index2 != -1)
469 startIndex =
f.fcIndex(startIndex);
472 if (index1 == -1 || index2 == -1)
479 scalar maxCos = -GREAT;
483 const vector& rightEdge = edges[right(size, fp)];
484 const vector leftEdge = -edges[left(size, fp)];
486 scalar
cos = rightEdge & leftEdge;
495 <<
"Cannot find valid diagonal on face " <<
f
496 <<
" with points " << UIndirectList<point>(
points,
f)
498 <<
"Returning naive triangulation starting from "
499 <<
f[maxIndex] <<
" which might not be correct for a"
500 <<
" concave or warped face" <<
endl;
503 label fp =
f.fcIndex(maxIndex);
505 for (label i = 0; i < size-2; i++)
507 label nextFp =
f.fcIndex(fp);
509 triFace& tri = operator[](triI++);
510 tri[0] =
f[maxIndex];
522 <<
"Cannot find valid diagonal on face " <<
f
523 <<
" with points " << UIndirectList<point>(
points,
f)
525 <<
"Returning empty triFaceList" <<
endl;
540 diff = index2 - index1;
545 diff = index2 + size - index1;
548 label nPoints1 =
diff + 1;
549 label nPoints2 = size -
diff + 1;
551 if (nPoints1 == size || nPoints2 == size)
554 <<
"Illegal split of face:" <<
f
555 <<
" with points " << UIndirectList<point>(
points,
f)
556 <<
" at indices " << index1 <<
" and " << index2
562 face face1(nPoints1);
564 label faceVertI = index1;
565 for (
int i = 0; i < nPoints1; i++)
567 face1[i] =
f[faceVertI];
568 faceVertI =
f.fcIndex(faceVertI);
572 face face2(nPoints2);
575 for (
int i = 0; i < nPoints2; i++)
577 face2[i] =
f[faceVertI];
578 faceVertI =
f.fcIndex(faceVertI);
619 bool valid =
split(fallBack,
points,
f, avgNormal, triI);