39 { triangulationMode::tmFan,
"fan" },
40 { triangulationMode::tmMesh,
"mesh" },
43 Foam::scalar Foam::faceAreaIntersect::tol = 1
e-6;
47 void Foam::faceAreaIntersect::triSliceWithPlane
51 FixedList<triPoints, 10>& tris,
57 FixedList<scalar, 3> d;
67 d[i] = pln.signedDistance(tri[i]);
69 if (
mag(d[i]) < tol*len)
95 || ((nPos == 2) && (nCoPlanar == 1))
96 || ((nPos == 1) && (nCoPlanar == 2))
111 else if ((nPos == 2) && (nCoPlanar == 0))
130 label i1 = d.fcIndex(i0);
131 label i2 = d.fcIndex(i1);
134 point p01 = planeIntersection(d, tri, i0, i1);
135 point p02 = planeIntersection(d, tri, i0, i2);
139 setTriPoints(tri[i1], tri[i2], p02, nTris, tris);
140 setTriPoints(tri[i1], p02, p01, nTris, tris);
163 label i1 = d.fcIndex(i0);
164 label i2 = d.fcIndex(i1);
167 point p01 = planeIntersection(d, tri, i1, i0);
168 point p02 = planeIntersection(d, tri, i2, i0);
171 setTriPoints(tri[i0], p01, p02, nTris, tris);
194 point p01 = planeIntersection(d, tri, i1, i0);
197 if (d.fcIndex(i0) == i1)
199 setTriPoints(tri[i0], p01, tri[i2], nTris, tris);
203 setTriPoints(tri[i0], tri[i2], p01, nTris, tris);
222 Foam::scalar Foam::faceAreaIntersect::triangleIntersect
224 const triPoints& src,
232 FixedList<triPoints, 10> workTris1;
233 label nWorkTris1 = 0;
235 FixedList<triPoints, 10> workTris2;
236 label nWorkTris2 = 0;
241 const scalar srcArea(triArea(src));
242 if (srcArea < ROOTVSMALL)
248 const scalar t =
sqrt(srcArea);
255 scalar
s =
mag(tgt1 - tgt0);
264 const vector n0((tgt0 - tgt1)^(-
s*
n));
265 const scalar magSqrN0(
magSqr(n0));
266 if (magSqrN0 < ROOTVSMALL)
275 plane pl0(tgt0, n0/
Foam::sqrt(magSqrN0),
false);
276 triSliceWithPlane(src, pl0, workTris1, nWorkTris1, t);
290 scalar
s =
mag(tgt2 - tgt1);
295 const vector n1((tgt1 - tgt2)^(-
s*
n));
296 const scalar magSqrN1(
magSqr(n1));
298 if (magSqrN1 < ROOTVSMALL)
307 plane pl1(tgt1, n1/
Foam::sqrt(magSqrN1),
false);
311 for (
label i = 0; i < nWorkTris1; ++i)
313 triSliceWithPlane(workTris1[i], pl1, workTris2, nWorkTris2, t);
328 scalar
s =
mag(tgt2 - tgt0);
333 const vector n2((tgt2 - tgt0)^(-
s*
n));
334 const scalar magSqrN2(
magSqr(n2));
336 if (magSqrN2 < ROOTVSMALL)
345 plane pl2(tgt2, n2/
Foam::sqrt(magSqrN2),
false);
349 for (
label i = 0; i < nWorkTris2; ++i)
351 triSliceWithPlane(workTris2[i], pl2, workTris1, nWorkTris1, t);
362 for (
label i = 0; i < nWorkTris1; ++i)
364 area += triArea(workTris1[i]);
366 if (cacheTriangulation_)
368 triangles_.append(workTris1[i]);
388 const bool cacheTriangulation
396 cacheTriangulation_(cacheTriangulation),
397 triangles_(cacheTriangulation ? 10 : 0)
411 faceTris.
resize(
f.nTriangles());
415 case triangulationMode::tmFan:
417 for (
label i = 0; i <
f.nTriangles(); ++i)
419 faceTris[i] =
face(3);
420 faceTris[i][0] =
f[0];
421 faceTris[i][1] =
f[i + 1];
422 faceTris[i][2] =
f[i + 2];
427 case triangulationMode::tmMesh:
429 const label nFaceTris =
f.nTriangles();
431 label nFaceTris1 = 0;
432 const label nFaceTris2 =
f.triangles(
points, nFaceTris1, faceTris);
434 if (nFaceTris != nFaceTris1 || nFaceTris != nFaceTris2)
437 <<
"The numbers of reported triangles in the face do not "
438 <<
"match that generated by the triangulation"
454 if (cacheTriangulation_)
460 scalar totalArea = 0.0;
461 for (
const face& triA : trisA_)
463 triPoints tpA = getTriPoints(pointsA_, triA,
false);
465 for (
const face& triB : trisB_)
503 const scalar threshold
507 scalar totalArea = 0.0;
508 for (
const face& triA : trisA_)
510 const triPoints tpA = getTriPoints(pointsA_, triA,
false);
512 for (
const face& triB : trisB_)
539 if (totalArea > threshold)