38 namespace triIntersect
60 return bStar == 0 ? 0 :
max(a*
sign(
b), 0)/bStar;
92 if (groups[i] ==
group)
106 if (groups[i] ==
group)
117 if (groups[i] ==
group)
126 if (groups[i] ==
group)
130 phi =
max(phi,
n*
x[i]/(
n*
x[i] - 1));
136 phi =
max(phi,
n*xn/(
n*xn - 1));
141 if (groups[i] ==
group)
143 y[i] =
min(
max((1 - phi)*
x[i] + phi/
n, 0), 1);
168 const vector CuCtut = Cu + Ctu*t;
173 CuCtut & CuCtut, CuCtut & CvCtvt,
174 CvCtvt & CuCtut, CvCtvt & CvCtvt
182 const scalar detA =
det(
A);
214 ((
C ^ Ctu) & Ctv) + ((Ct ^ Cu) & Ctv) + ((Ct ^ Ctu) &
Cv),
215 ((
C ^ Cu) & Ctv) + ((
C ^ Ctu) &
Cv) + ((Ct ^ Cu) &
Cv),
226 if (
mag(tRoots[tRooti]) > great)
continue;
243 tuvs[nTuvs ++] = tuv;
248 for (
label i = 0; i < nTuvs; ++ i)
250 const vector tuvOld = tuvs[i];
256 for (
label i = 0; i < nTuvs - 1; ++ i)
258 for (
label j = 0; j < nTuvs - 1; ++ j)
260 if (tuvClippage[j] > tuvClippage[j + 1])
262 Swap(tuvs[j], tuvs[j + 1]);
263 Swap(tuvClippage[j], tuvClippage[j + 1]);
270 for (
label i = 0; i < nTuvs; ++ i)
272 const scalar t = tuvs[i].x(), u = tuvs[i].y(), v = tuvs[i].z();
274 const scalar magSqrF =
275 magSqr(
C + Ct*t + Cu*u +
Cv*v + Ctu*t*u + Ctv*t*v);
276 const scalar magSqrErrF =
287 if (magSqrF < magSqrErrF)
311 const triPointRef tgtTri(tgtPs[0], tgtPs[1], tgtPs[2]);
318 srcFwd[srcPi] = (srcNs[srcPi] & tgtN) <
maxDot;
329 const label tgtPi0 = tgtEi, tgtPi1 = (tgtEi + 1) % 3;
335 {tgtPs[tgtPi0], tgtPs[tgtPi1]},
338 srcInTgtEdge[srcPi][tgtEi] = o > 0 ? +1 : -1;
349 srcInTgtEdge[srcPi] = {-1, -1, -1};
381 const label srcPi0 = srcEi, srcPi1 = (srcEi + 1) % 3;
383 if (srcFwd[srcPi0] != srcFwd[srcPi1])
390 maxDot - (srcNs[srcPi0] & tgtN),
391 (srcNs[srcPi1] - srcNs[srcPi0]) & tgtN
393 const point srcP = (1 - srcT)*srcPs[srcPi0] + srcT*srcPs[srcPi1];
394 const vector srcN = (1 - srcT)*srcNs[srcPi0] + srcT*srcNs[srcPi1];
398 const label tgtPi0 = tgtEi, tgtPi1 = (tgtEi + 1) % 3;
404 {tgtPs[tgtPi0], tgtPs[tgtPi1]},
407 srcInTgtEdge[srcFwd[srcPi0] ? srcPi1 : srcPi0][tgtEi] =
416 srcInTgtTri[srcPi] =
count(srcInTgtEdge[srcPi], 1) == 3 ? +1 : -1;
434 const triPointRef tgtTri(tgtPs[0], tgtPs[1], tgtPs[2]);
443 const label srcPi0 = srcEi, srcPi1 = (srcEi + 1) % 3;
447 {srcPs[srcPi0], srcPs[srcPi1]},
448 {srcNs[srcPi0], srcNs[srcPi1]},
452 tgtInSrcEdge[tgtPi][srcEi] = o > 0 ? +1 : -1;
459 tgtInSrcTri[tgtPi] =
count(tgtInSrcEdge[tgtPi], 1) == 3 ? +1 : -1;
465 if (tgtInSrcTri[tgtPi] == 1)
470 const bool tgtFwd = (srcN & tgtN) <
maxDot;
471 tgtInSrcTri[tgtPi] = tgtFwd ? +1 : -1;
486 forAll(thisOtherPis, thisPi)
488 const label otherPi = thisOtherPis[thisPi];
492 const label otherEi0 = (otherPi + 2) % 3, otherEi1 = otherPi;
494 thisInOtherTri[thisPi] = 0;
496 thisInOtherEdge[thisPi][otherEi0] = 0;
497 thisInOtherEdge[thisPi][otherEi1] = 0;
515 order[nVisited ++] = i0;
571 auto branch = [&](
const bool isSrcEdge)
592 && (branch(
true) || branch(
false))
604 order[-- nVisited] = -1;
623 srcPointNormals.
resize(pointLocations.
size());
626 forAll(pointLocations, pointi)
628 const location& l = pointLocations[pointi];
636 srcPoints[pointi] = srcP;
637 srcPointNormals[pointi] = srcN;
638 tgtPoints[pointi] = tgtP;
647 srcPoints[pointi] = srcP;
648 srcPointNormals[pointi] = srcN;
659 tgtPoints[pointi] = tgtP;
663 const label srcPi0 = l.
srcEdgei(), srcPi1 = (srcPi0 + 1) % 3;
664 const label tgtPi0 = l.
tgtEdgei(), tgtPi1 = (tgtPi0 + 1) % 3;
669 {srcPs[srcPi0], srcPs[srcPi1]},
670 {srcNs[srcPi0], srcNs[srcPi1]},
671 {tgtPs[tgtPi0], tgtPs[tgtPi1]}
673 const scalar srcT = ts.
first(), tgtT = ts.
second();
676 (1 - srcT)*srcPs[srcPi0] + srcT*srcPs[srcPi1];
677 srcPointNormals[pointi] =
678 (1 - srcT)*srcNs[srcPi0] + srcT*srcNs[srcPi1];
680 (1 - tgtT)*tgtPs[tgtPi0] + tgtT*tgtPs[tgtPi1];
701 scalar lengthScale = 0;
702 for (
label i = 0; i < 3; ++ i)
704 lengthScale =
max(lengthScale,
mag(srcPs[i] - srcPs[(i + 1) % 3]));
707 const label nu = nEdge, nv = nNormal;
708 const scalar
u0 = 0, u1 = 1;
709 const scalar v0 = -lNormal/2*lengthScale, v1 = lNormal/2*lengthScale;
712 for (
label i = 0; i < 3; ++ i)
714 const point& p0 = srcPs[i], p1 = srcPs[(i + 1) % 3];
715 const vector& n0 = srcNs[i], n1 = srcNs[(i + 1) % 3];
717 for (
label iu = 0; iu <= nu; ++ iu)
719 const scalar u =
u0 + (u1 -
u0)*scalar(iu)/nu;
720 for (
label iv = 0; iv <= nv; ++ iv)
722 const scalar v = v0 + (v1 - v0)*scalar(iv)/nv;
723 const vector x = p0 + (p1 - p0)*u + (n0 + (n1 - n0)*u)*v;
724 ps[i*(nu + 1)*(nv + 1) + iu*(nv + 1) + iv] =
x;
730 for (
label i = 0; i < 3; ++ i)
732 for (
label iu = 0; iu < nu; ++ iu)
734 for (
label iv = 0; iv < nv; ++ iv)
736 fs[i*nu*nv + iu*nv + iv] =
739 i*(nu + 1)*(nv + 1) + (nv + 1)*iu + iv,
740 i*(nu + 1)*(nv + 1) + (nv + 1)*iu + iv + 1,
741 i*(nu + 1)*(nv + 1) + (nv + 1)*(iu + 1) + iv + 1,
742 i*(nu + 1)*(nv + 1) + (nv + 1)*(iu + 1) + iv
769 const tensor A(srcPs[1] - srcPs[0], srcNs[0], srcNs[1]);
770 const scalar detA =
det(
A);
774 const vector detAY =
T & (tgtP - srcPs[0]);
777 const scalar
offset = Yx*detAY.
y() - (1 - Yx)*detAY.
z();
788 const bool srcDirection
809 const tensor A(srcN, tgtPs[1] - tgtPs[0], srcN^(tgtPs[1] - tgtPs[0]));
810 const scalar detA =
det(
A);
814 const vector detAY =
T & (tgtPs[0] - srcP);
827 const bool tgtDirection
921 scalar tgtT0 = 0, tgtT1 = 1;
926 const scalar
s = o0 > o1 ? +1 : -1;
930 const scalar tgtT = (tgtT0 + tgtT1)/2;
931 const vector tgtP = (1 - tgtT)*tgtPs[0] + tgtT*tgtPs[1];
945 tgtT = (tgtT0 + tgtT1)/2;
948 const vector srcDP = srcPs[1] - srcPs[0];
949 const vector tgtP = (1 - tgtT)*tgtPs[0] + tgtT*tgtPs[1];
951 const tensor A(srcDP, srcNs[0], srcNs[1]);
952 const scalar detA =
det(
A);
956 const scalar magDetAYx =
sign(detA)*(Tx & (tgtP - srcPs[0]));
960 const vector srcN = (1 - srcTStar)*srcNs[0] + srcTStar*srcNs[1];
962 const vector srcDPPerpN = srcDP - (srcDP & srcN)*srcN;
967 (tgtP - srcPs[0]) & srcDPPerpN,
1002 const scalar srcT0 = 1 - srcTs.x() - srcTs.y();
1003 return srcT0*srcPNs[0] + srcTs.x()*srcPNs[1] + srcTs.y()*srcPNs[2];
1006 auto srcEdgePNs = [srcPN]
1013 return Pair<vector>(srcPN(srcPNs, srcTs0), srcPN(srcPNs, srcTs1));
1020 srcEdgePNs(srcPs, srcTs0, srcTs1),
1021 srcEdgePNs(srcNs, srcTs0, srcTs1),
1032 if (offsets[
findMin(offsets)] >= 0 || offsets[
findMax(offsets)] <= 0)
1042 srcPs[1] - srcPs[0],
1043 srcPs[2] - srcPs[0],
1044 srcNs[1] - srcNs[0],
1045 srcNs[2] - srcNs[0],
1058 const scalar
sign = offsets[
findMin(offsets)] > 0 ? +1 : -1;
1060 vector2D srcTsA(0, 0), srcTsB(1, 0), srcTsC(0, 1);
1064 const vector2D srcTsAB = (srcTsA + srcTsB)/2;
1065 const vector2D srcTsBC = (srcTsB + srcTsC)/2;
1066 const vector2D srcTsCA = (srcTsC + srcTsA)/2;
1068 const scalar oA =
sign*
offset(srcTsCA, srcTsAB);
1069 const scalar oB =
sign*
offset(srcTsAB, srcTsBC);
1070 const scalar oC =
sign*
offset(srcTsBC, srcTsCA);
1074 if (offsets[offsetMini] > 0)
1080 else if (offsetMini == 0)
1085 else if (offsetMini == 1)
1090 else if (offsetMini == 2)
1097 srcT = (srcTsA[0] + srcTsB[0] + srcTsC[0])/3;
1098 srcU = (srcTsA[1] + srcTsB[1] + srcTsC[1])/3;
1107 const label srcEi0 = (srcEi + 2) % 3, srcEi1 = (srcEi + 1) % 3;
1112 && offsets[srcEi0] >= 0
1113 && offsets[srcEi1] >= 0
1116 const label srcPi0 = srcEi, srcPi1 = (srcEi + 1) % 3;
1117 const label srcPiOpp = (srcEi + 2) % 3;
1125 srcPs[srcPi0] - tgtP,
1126 srcPs[srcPi1] - srcPs[srcPi0],
1127 srcPs[srcPi0] - srcPs[srcPiOpp],
1129 srcPs[srcPi1] - srcPs[srcPi0],
1130 srcNs[srcPi1] - srcNs[srcPi0],
1143 const vector2D srcTsOpp(srcPiOpp == 1, srcPiOpp == 2);
1144 const vector2D srcTs0(srcPi0 == 1, srcPi0 == 2);
1145 const vector2D srcTs1(srcPi1 == 1, srcPi1 == 2);
1147 scalar srcT0 = 0, srcT1 = 1;
1151 const scalar srcT = (srcT0 + srcT1)/2;
1152 const vector2D srcTs01(srcTs0*(1 - srcT) + srcTs1*srcT);
1154 const scalar o =
offset(srcTsOpp, srcTs01);
1166 srcT = (srcT0 + srcT1)/2;
1170 srcPs[srcPi0] - tgtP,
1171 srcPs[srcPi1] - srcPs[srcPi0],
1172 srcPs[srcPi0] - srcPs[srcPiOpp],
1174 srcPs[srcPi1] - srcPs[srcPi0],
1175 srcNs[srcPi1] - srcNs[srcPi0],
1183 y[srcPiOpp] = - srcU;
1184 y[srcPi0] = (1 + srcU)*(1 - srcT);
1185 y[srcPi1] = (1 + srcU)*srcT;
1193 const label srcEiOpp = (srcPi + 1) % 3;
1194 const label srcEi0 = (srcPi + 2) % 3, srcEi1 = srcPi;
1198 offsets[srcEiOpp] >= 0
1199 && offsets[srcEi0] <= 0
1200 && offsets[srcEi1] <= 0
1204 const label srcPi0 = (srcPi + 2) % 3, srcPi1 = (srcPi + 1) % 3;
1208 srcPs[srcPi] - srcPs[srcPi0],
1209 srcPs[srcPi] - srcPs[srcPi1],
1213 const scalar detA =
A.x() &
T0;
1224 y[srcPi] = 1 + srcT + srcU;
1231 <<
"Point " << tgtP <<
" could not be classified within triangle "
1232 << srcPs <<
" with projection normals " << srcNs <<
exit(
FatalError);
1233 return barycentric2D::uniform(
NaN);
1244 const tensor A(tgtPs[1] - tgtPs[0], tgtPs[2] - tgtPs[0], - srcN);
1245 const scalar detA =
det(
A);
1250 const vector detAY = (
T & (srcP - tgtPs[0])) +
vector(detA, 0, 0);
1255 maxMagDetAY/vGreat <
mag(detA)
1257 : detAY/maxMagDetAY*vGreat;
1277 const word& writePrefix
1280 const bool write = writePrefix != word::null;
1298 srcInTgt(srcPs, srcNs, srcOwns, tgtPs, tgtOwns, srcInTgtEdge, srcInTgtTri);
1299 thisIsOther(srcTgtPis, srcInTgtEdge, srcInTgtTri);
1300 tgtInSrc(srcPs, srcNs, srcOwns, tgtPs, tgtOwns, tgtInSrcEdge, tgtInSrcTri);
1301 thisIsOther(tgtSrcPis, tgtInSrcEdge, tgtInSrcTri);
1305 if (
count(srcTgtPis, -1) != 3)
1310 <<
indent <<
"srcInTgtEdge=" << srcInTgtEdge <<
endl;
1311 if (
count(tgtSrcPis, -1) != 3)
1316 <<
indent <<
"tgtInSrcEdge=" << tgtInSrcEdge <<
endl;
1320 bool noIntersection =
false;
1326 bool outside =
true;
1329 if (tgtInSrcEdge[tgtPi][srcEi] == 1)
1337 noIntersection =
true;
1345 if (
count(srcInTgtEdge, {-1, -1, -1}) == 3)
1347 noIntersection =
true;
1352 if (!noIntersection)
1357 auto addPointLocations = [&srcInTgtTri,&pointLocations]
1361 const bool add = true
1364 if (!pointLocations.
empty())
1370 const label srcEi0 =
1374 const label tgtEi0 =
1379 const label srcEi1 =
1383 const label tgtEi1 =
1386 : (l1.tgtPointi() + 2) % 3;
1396 label srcEj = srcEi0;
1398 srcEj = (srcEj + 2) % 3
1401 if (srcInTgtTri[srcEj] != 1)
1406 pointLocations.
append(location::srcPoint(srcEj));
1417 pointLocations.
append(l1);
1421 pointLocations.
append(l2);
1428 auto inTriToOut = [&addPointLocations,&srcInTgtEdge]
1431 const label tgtOutSrcEi1,
1432 const label tgtOutSrcPi1,
1439 : srcInTgtEdge[tgtOutSrcPi1][tgtEi] == 1
1440 ? (tgtOutSrcPi1 + 2*
reverse) % 3
1441 : (tgtOutSrcPi1 + 2*!
reverse) % 3;
1443 return addPointLocations(location::intersection(srcEi, tgtEi));
1448 auto isPointToOutEdge = [&addPointLocations,&srcInTgtEdge]
1451 const label tgtIsSrcPi0,
1452 const label tgtOutSrcEi1,
1456 const label srcEi0Next = (tgtIsSrcPi0 + 2*
reverse) % 3;
1457 const label srcEi0Opp = (tgtIsSrcPi0 + 1) % 3;
1461 srcInTgtEdge[(tgtIsSrcPi0 + 1) % 3][tgtEi] == -1
1462 && srcInTgtEdge[(tgtIsSrcPi0 + 2) % 3][tgtEi] == -1
1468 if (tgtOutSrcEi1 == srcEi0Next)
1473 if (tgtOutSrcEi1 == srcEi0Opp)
1478 location::intersection(srcEi0Opp, tgtEi)
1487 auto isPointToOutCorner = []
1490 const label tgtIsSrcPi0,
1491 const label tgtOutSrcPi1,
1495 return tgtOutSrcPi1 != (tgtIsSrcPi0 + 1 +
reverse) % 3;
1499 auto outEdgeToOutEdge = [&addPointLocations,&srcInTgtEdge]
1502 const label tgtOutSrcEi0,
1503 const label tgtOutSrcEi1
1506 const label srcPi = (5 - tgtOutSrcEi0 - tgtOutSrcEi1) % 3;
1510 (tgtOutSrcEi0 != (tgtOutSrcEi1 + 1) % 3)
1511 && (srcInTgtEdge[srcPi][tgtEi] != 1)
1519 (tgtOutSrcEi0 == (tgtOutSrcEi1 + 1) % 3)
1520 != (srcInTgtEdge[srcPi][tgtEi] == 1)
1526 location::intersection(tgtOutSrcEi0, tgtEi),
1527 location::intersection(tgtOutSrcEi1, tgtEi)
1536 auto outEdgeToOutCorner = [&addPointLocations,&srcInTgtEdge]
1539 const label tgtOutSrcEi0,
1540 const label tgtOutSrcPi1,
1544 if (tgtOutSrcEi0 != (tgtOutSrcPi1 + 1) % 3)
1549 const label srcPi1Prev = (tgtOutSrcPi1 + 1 + !
reverse) % 3;
1551 if (srcInTgtEdge[srcPi1Prev][tgtEi] == -1)
1556 const label srcPi1Next = (tgtOutSrcPi1 + 1 +
reverse) % 3;
1558 if (srcInTgtEdge[srcPi1Next][tgtEi] == 1)
1563 location l1 = location::intersection(tgtOutSrcEi0, tgtEi);
1566 srcInTgtEdge[tgtOutSrcPi1][tgtEi] == 1
1567 ? (tgtOutSrcPi1 + 2*
reverse) % 3
1568 : (tgtOutSrcPi1 + 2*!
reverse) % 3;
1570 location l2 = location::intersection(srcEi, tgtEi);
1577 return addPointLocations(l1, l2);
1581 auto outCornerToOutCorner = []
1584 const label tgtOutSrcPi0,
1585 const label tgtOutSrcPi1
1588 return tgtOutSrcPi0 != (tgtOutSrcPi1 + 2) % 3;
1594 const label tgtPi0 = tgtEi, tgtPi1 = (tgtEi + 1) % 3;
1596 const bool tgtInSrcTri0 = tgtInSrcTri[tgtPi0] == 1;
1597 const bool tgtInSrcTri1 = tgtInSrcTri[tgtPi1] == 1;
1599 const label tgtIsSrcPi0 =
1600 tgtInSrcTri[tgtPi0] == 0 ? tgtSrcPis[tgtPi0] : -1;
1601 const label tgtIsSrcPi1 =
1602 tgtInSrcTri[tgtPi1] == 0 ? tgtSrcPis[tgtPi1] : -1;
1604 const label tgtOutSrcEi0 =
1605 count(tgtInSrcEdge[tgtPi0], -1) == 1
1608 const label tgtOutSrcEi1 =
1609 count(tgtInSrcEdge[tgtPi1], -1) == 1
1613 const label tgtOutSrcPi0 =
1614 count(tgtInSrcEdge[tgtPi0], -1) == 2
1615 ? (
findIndex(tgtInSrcEdge[tgtPi0], 1) + 2) % 3
1617 const label tgtOutSrcPi1 =
1618 count(tgtInSrcEdge[tgtPi1], -1) == 2
1619 ? (
findIndex(tgtInSrcEdge[tgtPi1], 1) + 2) % 3
1625 pointLocations.
append(location::tgtPoint(tgtPi0));
1627 if (tgtIsSrcPi0 != -1)
1633 location::srcTgtPoint(tgtIsSrcPi0, tgtPi0)
1637 pointLocations.
clear();
1645 (tgtInSrcTri0 && tgtInSrcTri1)
1646 || (tgtOutSrcEi0 != -1 && tgtOutSrcEi0 == tgtOutSrcEi1)
1647 || (tgtOutSrcPi0 != -1 && tgtOutSrcPi0 == tgtOutSrcPi1)
1655 (tgtInSrcTri0 && tgtIsSrcPi1 != -1)
1656 || (tgtIsSrcPi0 != -1 && tgtInSrcTri1)
1662 else if (tgtInSrcTri0 && (tgtOutSrcEi1 != -1 || tgtOutSrcPi1 != -1))
1666 if (!inTriToOut(tgtEi, tgtOutSrcEi1, tgtOutSrcPi1, 0))
1668 pointLocations.
clear();
1672 else if ((tgtOutSrcEi0 != -1 || tgtOutSrcPi0 != -1) && tgtInSrcTri1)
1675 if (!inTriToOut(tgtEi, tgtOutSrcEi0, tgtOutSrcPi0, 1))
1677 pointLocations.
clear();
1681 else if (tgtIsSrcPi0 != -1 && tgtIsSrcPi1 != -1)
1685 if (tgtIsSrcPi0 != (tgtIsSrcPi1 + 1) % 3)
1687 pointLocations.
clear();
1691 else if (tgtIsSrcPi0 != -1 && tgtOutSrcEi1 != -1)
1695 if (!isPointToOutEdge(tgtEi, tgtIsSrcPi0, tgtOutSrcEi1, 0))
1697 pointLocations.
clear();
1701 else if (tgtOutSrcEi0 != -1 && tgtIsSrcPi1 != -1)
1704 if (!isPointToOutEdge(tgtEi, tgtIsSrcPi1, tgtOutSrcEi0, 1))
1706 pointLocations.
clear();
1710 else if (tgtIsSrcPi0 != -1 && tgtOutSrcPi1 != -1)
1714 if (!isPointToOutCorner(tgtEi, tgtIsSrcPi0, tgtOutSrcPi1, 0))
1716 pointLocations.
clear();
1720 else if (tgtOutSrcPi0 != -1 && tgtIsSrcPi1 != -1)
1723 if (!isPointToOutCorner(tgtEi, tgtIsSrcPi1, tgtOutSrcPi0, 1))
1725 pointLocations.
clear();
1729 else if (tgtOutSrcEi0 != -1 && tgtOutSrcEi1 != -1)
1732 if (!outEdgeToOutEdge(tgtEi, tgtOutSrcEi0, tgtOutSrcEi1))
1734 pointLocations.
clear();
1738 else if (tgtOutSrcEi0 != -1 && tgtOutSrcPi1 != -1)
1742 if (!outEdgeToOutCorner(tgtEi, tgtOutSrcEi0, tgtOutSrcPi1, 0))
1744 pointLocations.
clear();
1748 else if (tgtOutSrcPi0 != -1 && tgtOutSrcEi1 != -1)
1751 if (!outEdgeToOutCorner(tgtEi, tgtOutSrcEi1, tgtOutSrcPi0, 1))
1753 pointLocations.
clear();
1757 else if (tgtOutSrcPi0 != -1 && tgtOutSrcPi1 != -1)
1760 if (!outCornerToOutCorner(tgtEi, tgtOutSrcPi0, tgtOutSrcPi1))
1762 pointLocations.
clear();
1770 pointLocations.
clear();
1776 if (!pointLocations.
empty())
1782 if (!addPointLocations(l,
location(),
false))
1784 pointLocations.
clear();
1799 if (pointLocations.
empty() &&
count(srcInTgtTri, 1) == 3)
1803 pointLocations.
append(location::srcPoint(srcPi));
static const Foam::dimensionedScalar A("A", Foam::dimPressure, 611.21)
static const Foam::dimensionedScalar B("B", Foam::dimless, 18.678)
scalar Cv(const scalar p, const scalar T) const
#define forAll(list, i)
Loop across all elements in list.
Templated 2D Barycentric derived from VectorSpace. Has 3 components, one of which is redundant.
Graphite solid properties.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
void resize(const label)
Alter the addressed list size.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
void clear()
Clear the addressed list, i.e. set the size to zero.
A 1D vector of objects of type <T> with a fixed size <Size>.
void size(const label)
Override size to be inconsistent with allocated storage.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
const Type & second() const
Return second.
const Type & first() const
Return first.
Templated storage for the roots of polynomial equations, plus flags to indicate the nature of the roo...
void type(const direction i, const rootType t)
Set the type of the i-th root.
Templated 2D tensor derived from VectorSpace adding construction from 4 components,...
A 2-tuple for storing two objects of different types.
T * first()
Return the first entry.
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
T & first()
Return the first element of the list.
bool empty() const
Return true if the UList is empty (ie, size() is zero)
T & last()
Return the last element of the list.
static Form uniform(const Cmpt &s)
Return a VectorSpace with all elements = s.
Cubic equation of the form a*x^3 + b*x^2 + c*x + d = 0.
Roots< 3 > roots() const
Get the roots.
A face is a list of labels corresponding to mesh vertices.
Selector class for relaxation factors, solver type and solution.
bool isTgtNotSrcPoint() const
Return whether the location is a target point and not a source.
bool isTgtPoint() const
Return whether the location is a target point.
label tgtEdgei() const
Return the target edge index.
bool isSrcNotTgtPoint() const
Return whether the location is a source point and not a target.
label srcPointi() const
Return the source point index.
label tgtPointi() const
Return the target point index.
bool isSrcAndTgtPoint() const
Return whether the location is a source point and a target point.
bool isIntersection() const
Return whether the location is an intersection.
bool isSrcPoint() const
Return whether the location is a source point.
label srcEdgei() const
Return the source edge index.
A triangle primitive used to calculate face areas and swept volumes.
vector normal() const
Return unit normal.
A class for handling words, derived from string.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
const char *const group
Group name for atomic constants.
barycentric2D srcPointTgtTriIntersection(const point &srcP, const vector &srcN, const FixedList< point, 3 > &tgtPs)
Calculate the intersection of a projected source point with a target.
scalar srcPointTgtEdgeOffset(const point &srcP, const vector &srcN, const Pair< point > &tgtPs)
Calculate the signed offset of a projected source point in relation to a.
void srcInTgt(const FixedList< point, 3 > &srcPs, const FixedList< vector, 3 > &srcNs, const FixedList< bool, 3 > &srcOwns, const FixedList< point, 3 > &tgtPs, const FixedList< bool, 3 > &tgtOwns, FixedList< FixedList< label, 3 >, 3 > &srcInTgtEdge, FixedList< label, 3 > &srcInTgtTri)
Calculate whether the points of the given source triangle project inside or.
barycentric2D srcTriTgtPointIntersection(const FixedList< point, 3 > &srcPs, const FixedList< vector, 3 > &srcNs, const point &tgtP)
Calculate the intersection of a target point with a source triangle's.
Tuple2< bool, vector > solveProjection(const vector &C, const vector &Ct, const vector &Cu, const vector &Cv, const vector &Ctu, const vector &Ctv, const FixedList< label, 3 > groups)
Solve a projection equation.
void tgtInSrc(const FixedList< point, 3 > &srcPs, const FixedList< vector, 3 > &srcNs, const FixedList< bool, 3 > &srcOwns, const FixedList< point, 3 > &tgtPs, const FixedList< bool, 3 > &tgtOwns, FixedList< FixedList< label, 3 >, 3 > &tgtInSrcEdge, FixedList< label, 3 > &tgtInSrcTri)
Calculate whether the points of the given target triangle project inside or.
void writeTriProjection(const word &name, const FixedList< point, 3 > &srcTriPs, const FixedList< vector, 3 > &srcTriNs, const label nEdge=20, const label nNormal=20, const scalar lNormal=0.5)
Write a VTK file of a triangle projection.
scalar protectedDivideAndClip01(const scalar a, const scalar b)
Divide two numbers, protect the result from overflowing, and clip the.
void thisIsOther(const FixedList< label, 3 > &thisOtherPis, FixedList< FixedList< label, 3 >, 3 > &thisInOtherEdge, FixedList< label, 3 > &thisInOtherTri)
Override results of the srcInTgt/tgtInSrc calculations with explicit.
Ostream & operator<<(Ostream &os, const FixedList< FixedList< Type, 3 >, 3 > &l)
Print 3x3 FixedListLists on one line.
vector clipped01(const vector x, const FixedList< label, 3 > groups)
Clip the given vector between values of 0 and 1, and also clip one minus.
const scalar maxDot
The maximum dot product between a source point normal and a target plane.
scalar srcEdgeTgtPointOffset(const Pair< point > &srcPs, const Pair< vector > &srcNs, const point &tgtP)
Calculate the signed offset of a target point in relation to a projected.
vector solveProjectionGivenT(const vector &C, const vector &Ct, const vector &Cu, const vector &Cv, const vector &Ctu, const vector &Ctv, const FixedList< label, 3 > groups, const scalar t)
Solve a projection equation given a value of the t variable.
Type tgtTriInterpolate(const barycentric2D &y, const FixedList< Type, 3 > &tgtPsis)
Use the coordinates obtained from srcPointTgtTriIntersection to interpolate.
void writePolygon(const word &name, const PointField &ps)
Write a VTK file of a polygon.
scalar protectedDivide(const scalar a, const scalar b)
Divide two numbers and protect the result from overflowing.
bool orderLocations(const UList< location > &locations, bool isSrcEdge, const label i0, label &nVisited, boolList &visited, labelList &order)
Order intersection locations into a polygon.
Pair< scalar > srcEdgeTgtEdgeIntersection(const Pair< point > &srcPs, const Pair< vector > &srcNs, const Pair< point > &tgtPs)
Calculate the intersection of a target edge with a source edge's.
void generateGeometryForLocations(const FixedList< point, 3 > &srcPs, const FixedList< vector, 3 > &srcNs, const FixedList< point, 3 > &tgtPs, DynamicList< point > &srcPoints, DynamicList< vector > &srcPointNormals, DynamicList< point > &tgtPoints, const DynamicList< location > &pointLocations)
Type srcTriInterpolate(const barycentric2D &y, const FixedList< Type, 3 > &tgtPsis)
Use the coordinates obtained from srcTriTgtPointIntersection to interpolate.
void intersectTris(const FixedList< point, 3 > &srcPs, const FixedList< vector, 3 > &srcNs, const FixedList< bool, 3 > &srcOwns, const FixedList< label, 3 > &srcTgtPis, const FixedList< point, 3 > &tgtPs, const FixedList< bool, 3 > &tgtOwns, const FixedList< label, 3 > &tgtSrcPis, DynamicList< point > &srcPoints, DynamicList< vector > &srcPointNormals, DynamicList< point > &tgtPoints, DynamicList< location > &pointLocations, const bool debug, const word &writePrefix=word::null)
Construct the intersection of a source triangle's projected volume and a.
void write(std::ostream &os, const bool binary, List< floatScalar > &fField)
Write floats ascii or binary.
Vector2D< scalar > vector2D
vector2D obtained from generic Vector2D
dimensionedScalar det(const dimensionedSphericalTensor &dt)
errorManipArg< error, int > exit(error &err, const int errNo=1)
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
List< label > labelList
A List of labels.
label log2(label i)
Return the log base 2 by successive bit-shifting of the given label.
dimensionedScalar sign(const dimensionedScalar &ds)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
label findMax(const ListType &, const label start=0)
Find index of max element (and larger than given element).
Ostream & endl(Ostream &os)
Add newline and flush stream.
void inplaceReverseList(ListType &list)
Inplace reversal of a list using Swap.
Barycentric2D< scalar > barycentric2D
A scalar version of the templated Barycentric2D.
Ostream & incrIndent(Ostream &os)
Increment the indent level.
void reverse(UList< T > &, const label n)
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
dimensionedSymmTensor cof(const dimensionedSymmTensor &dt)
Vector< scalar > vector
A scalar version of the templated Vector.
List< labelList > labelListList
A List of labelList.
void add(FieldField< Field1, typename typeOfSum< Type1, Type2 >::type > &f, const FieldField< Field1, Type1 > &f1, const FieldField< Field2, Type2 > &f2)
dimensioned< scalar > mag(const dimensioned< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
dimensionSet cmptMag(const dimensionSet &)
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
label findMin(const ListType &, const label start=0)
Find index of min element (and less than given element).
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
label count(const ListType &l, typename ListType::const_reference x)
Count the number of occurrences of a value in a list.
Cmpt cmptSum(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
void offset(label &lst, const label o)
word name(const complex &)
Return a string representation of a complex.
Ostream & indent(Ostream &os)
Indent stream.
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
dimensionedScalar cos(const dimensionedScalar &ds)
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
Functions with which to perform an intersection of a pair of triangles; the source and target....