33 const Foam::scalar Foam::faceTriangulation::edgeRelTol = 1
e-6;
37 Foam::label Foam::faceTriangulation::right(
const label, label i)
44 Foam::label Foam::faceTriangulation::left(
const label size, label i)
46 return i ? i-1 : size-1;
59 auto& edges = tedges.ref();
66 vector vec(nextPt - thisPt);
68 edges[i] = vec.normalise();
76 void Foam::faceTriangulation::calcHalfAngle
86 scalar
cos =
clamp((e0 & e1), -1, 1);
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)
142 result.hitPoint(intersectPt);
143 result.setDistance(
mag(intersectPt - rayOrigin));
152 bool Foam::faceTriangulation::triangleContainsPoint
165 if ((area01Pt > 0) && (area12Pt > 0) && (area20Pt > 0))
169 else if ((area01Pt < 0) && (area12Pt < 0) && (area20Pt < 0))
181 void Foam::faceTriangulation::findDiagonal
187 const label startIndex,
195 const vector& rightE = edges[right(
f.
size(), startIndex)];
196 const vector leftE = -edges[left(
f.
size(), startIndex)];
199 scalar cosHalfAngle = GREAT;
200 scalar sinHalfAngle = GREAT;
201 calcHalfAngle(normal, rightE, leftE, cosHalfAngle, sinHalfAngle);
205 + sinHalfAngle*(normal ^ rightE)
216 label faceVertI =
f.
fcIndex(startIndex);
220 scalar minPosOnEdge = GREAT;
222 for (label i = 0; i <
f.
size() - 2; i++)
236 if (inter.hit() && inter.distance() < minInter.distance())
239 minIndex = faceVertI;
240 minPosOnEdge = posOnEdge;
258 const label leftIndex = minIndex;
259 const label rightIndex =
f.
fcIndex(minIndex);
265 if (
mag(minPosOnEdge) < edgeRelTol &&
f.
fcIndex(startIndex) != leftIndex)
274 mag(minPosOnEdge - 1) < edgeRelTol
275 &&
f.
fcIndex(rightIndex) != startIndex
292 scalar maxCos = -GREAT;
296 for (label i = 0; i <
f.
size() - 3; i++)
302 (faceVertI == leftIndex)
303 || (faceVertI == rightIndex)
304 || (triangleContainsPoint(normal, startPt, leftPt, rightPt, pt))
309 vector edgePt0 = (pt - startPt);
312 scalar
cos = rayDir & edgePt0;
316 minIndex = faceVertI;
328 if (
f.
fcIndex(startIndex) != leftIndex)
348 Foam::label Foam::faceTriangulation::findStart
355 const label size =
f.
size();
357 scalar minCos = GREAT;
362 const vector& rightEdge = edges[right(size, fp)];
363 const vector leftEdge = -edges[left(size, fp)];
365 if (((rightEdge ^ leftEdge) & normal) < ROOTVSMALL)
367 scalar
cos = rightEdge & leftEdge;
383 const vector& rightEdge = edges[right(size, fp)];
384 const vector leftEdge = -edges[left(size, fp)];
386 scalar
cos = rightEdge & leftEdge;
403 bool Foam::faceTriangulation::split
412 const label size =
f.
size();
417 <<
"Illegal face:" <<
f 418 <<
" with points " << UIndirectList<point>(
points,
f)
426 triFace& tri = operator[](triI++);
438 tmp<vectorField> tedges(calcEdges(
f,
points));
441 label startIndex = findStart(
f, edges, normal);
460 if (index1 != -1 && index2 != -1)
470 if (index1 == -1 || index2 == -1)
477 scalar maxCos = -GREAT;
481 const vector& rightEdge = edges[right(size, fp)];
482 const vector leftEdge = -edges[left(size, fp)];
484 scalar
cos = rightEdge & leftEdge;
493 <<
"Cannot find valid diagonal on face " <<
f 494 <<
" with points " << UIndirectList<point>(
points,
f)
496 <<
"Returning naive triangulation starting from " 497 <<
f[maxIndex] <<
" which might not be correct for a" 498 <<
" concave or warped face" <<
endl;
503 for (label i = 0; i < size-2; i++)
507 triFace& tri = operator[](triI++);
508 tri[0] =
f[maxIndex];
520 <<
"Cannot find valid diagonal on face " <<
f 521 <<
" with points " << UIndirectList<point>(
points,
f)
523 <<
"Returning empty triFaceList" <<
endl;
538 diff = index2 - index1;
543 diff = index2 + size - index1;
546 label nPoints1 =
diff + 1;
547 label nPoints2 = size -
diff + 1;
549 if (nPoints1 == size || nPoints2 == size)
552 <<
"Illegal split of face:" <<
f 553 <<
" with points " << UIndirectList<point>(
points,
f)
554 <<
" at indices " << index1 <<
" and " << index2
560 face face1(nPoints1);
562 label faceVertI = index1;
563 for (
int i = 0; i < nPoints1; i++)
565 face1[i] =
f[faceVertI];
570 face face2(nPoints2);
573 for (
int i = 0; i < nPoints2; i++)
575 face2[i] =
f[faceVertI];
617 bool valid = split(fallBack,
points,
f, avgNormal, triI);
638 bool valid = split(fallBack,
points,
f,
n, triI);
List< triFace > triFaceList
List of triFace.
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
void size(const label n)
Older name for setAddressableSize.
A face is a list of labels corresponding to mesh vertices.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#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...
constexpr char nl
The newline '\n' character (0x0a)
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Vector< Cmpt > & normalise(const scalar tol=ROOTVSMALL)
Inplace normalise the vector by its magnitude.
faceTriangulation()
Construct null.
#define forAll(list, i)
Loop across all elements in list.
label fcIndex(const label i) const noexcept
The forward circular index. The next index in the list which returns to the first at the end of the l...
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
dimensionedScalar cos(const dimensionedScalar &ds)
void clear()
Clear the list, i.e. set size to zero.
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
Describes the interaction of a object and a (templated) point. It carries the info of a successful hi...
errorManip< error > abort(error &err)
dimensionedScalar sin(const dimensionedScalar &ds)
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
static bool split(const std::string &line, std::string &key, std::string &val)
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
Field< vector > vectorField
Specialisation of Field<T> for vector.
A class for managing temporary objects.
const volScalarField & p0
PointHit< point > pointHit
A PointHit with a 3D point.
dimensionSet clamp(const dimensionSet &a, const dimensionSet &range)
vector areaNormal() const
The area normal - with magnitude equal to area of triangle.
static constexpr const zero Zero
Global zero (0)