37 template<
class Po
int,
class Po
intRef>
38 pointHit triangle<Point, PointRef>::nearestPoint
40 const Point& baseVertex,
47 const vector D(baseVertex - P);
50 const scalar a = E0 & E0;
51 const scalar
b = E0 & E1;
52 const scalar c = E1 & E1;
55 const scalar
d = E0 & D;
56 const scalar
e = E1 & D;
57 const scalar
f = D & D;
60 const scalar
det = a*c - b*
b;
77 s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
83 t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
90 t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
97 s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
102 const scalar invDet = 1/
det;
114 const scalar tmp0 = b +
d;
115 const scalar tmp1 = c +
e;
119 const scalar numer = tmp1 - tmp0;
120 const scalar denom = a-2*b+c;
121 s = (numer >= denom ? 1 : numer/denom);
128 t = (tmp1 <= 0 ? 1 : (e >= 0 ? 0 : - e/c));
134 const scalar tmp0 = b +
d;
135 const scalar tmp1 = c +
e;
139 const scalar numer = tmp1 - tmp0;
140 const scalar denom = a-2*b+c;
141 s = (numer >= denom ? 1 : numer/denom);
148 s = (tmp1 <= 0 ? 1 : (d >= 0 ? 0 : - d/a));
154 const scalar numer = c+e-(b+
d);
161 const scalar denom = a-2*b+c;
162 s = (numer >= denom ? 1 : numer/denom);
177 baseVertex + s*E0 + t*E1,
180 Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f)
189 template<
class Po
int,
class Po
intRef>
203 template<
class Po
int,
class Po
intRef>
209 is >> a_ >> b_ >> c_;
215 is.
check(
"triangle::triangle(Istream& is)");
221 template<
class Po
int,
class Po
intRef>
227 template<
class Po
int,
class Po
intRef>
233 template<
class Po
int,
class Po
intRef>
240 template<
class Po
int,
class Po
intRef>
243 return (1.0/3.0)*(a_ + b_ + c_);
247 template<
class Po
int,
class Po
intRef>
254 template<
class Po
int,
class Po
intRef>
257 return 0.5*((b_ - a_)^(c_ - a_));
261 template<
class Po
int,
class Po
intRef>
264 scalar d1 = (c_ - a_)&(b_ - a_);
265 scalar d2 = -(c_ - b_)&(b_ - a_);
266 scalar d3 = (c_ - a_)&(c_ - b_);
272 scalar c = c1 + c2 + c3;
276 ((c2 + c3)*a_ + (c3 + c1)*b_ + (c1 + c2)*c_)/(2*c)
281 template<
class Po
int,
class Po
intRef>
284 scalar d1 = (c_ - a_) & (b_ - a_);
285 scalar d2 = - (c_ - b_) & (b_ - a_);
286 scalar d3 = (c_ - a_) & (c_ - b_);
288 scalar denom = d2*d3 + d3*d1 + d1*d2;
296 scalar a = (d1 + d2)*(d2 + d3)*(d3 + d1) / denom;
303 template<
class Po
int,
class Po
intRef>
317 template<
class Po
int,
class Po
intRef>
322 ((t.a_ - a_) & ((b_ - a_)^(c_ - a_)))
323 + ((t.b_ - b_) & ((c_ - b_)^(t.a_ - b_)))
324 + ((c_ - t.c_) & ((t.b_ - t.c_)^(t.a_ - t.c_)))
326 + ((t.a_ - a_) & ((b_ - a_)^(c_ - a_)))
327 + ((b_ - t.b_) & ((t.a_ - t.b_)^(t.c_ - t.b_)))
328 + ((c_ - t.c_) & ((b_ - t.c_)^(t.a_ - t.c_)))
333 template<
class Po
int,
class Po
intRef>
386 hit = fastInter.
hit();
396 pInter = p + (q1&v)*q1;
401 scalar
dist = q1 & (pInter -
p);
403 const scalar planarPointTol =
419 && ((q1 &
normal()) < -VSMALL)
435 inter.
setPoint(nearestPoint(a_, E0, E1, p).rawPoint());
449 template<
class Po
int,
class Po
intRef>
458 const vector edge1 = b_ - a_;
459 const vector edge2 = c_ - a_;
462 const vector pVec = dir ^ edge2;
465 const scalar det = edge1 & pVec;
473 if (det < ROOTVSMALL)
482 if (det > -ROOTVSMALL && det < ROOTVSMALL)
489 const scalar inv_det = 1.0 /
det;
492 const vector tVec = orig-a_;
495 const scalar u = (tVec & pVec)*inv_det;
497 if (u < -tol || u > 1.0+tol)
504 const vector qVec = tVec ^ edge1;
507 const scalar v = (dir & qVec) * inv_det;
509 if (v < -tol || u + v > 1.0+tol)
516 const scalar t = (edge2 & qVec) * inv_det;
525 intersection.
setPoint(a_ + u*edge1 + v*edge2);
532 template<
class Po
int,
class Po
intRef>
543 return nearestPoint(a_, E0, E1, p);
547 template<
class Po
int,
class Po
intRef>
556 const vector E0 = b_ - a_;
557 const vector E1 = c_ - a_;
558 const vector n = 0.5*(E0 ^ E1);
566 if ((magX >= magY) && (magX >= magZ))
570 else if ((magY >= magX) && (magY >= magZ))
580 label i1 = (i0 + 1) % 3;
581 label i2 = (i1 + 1) % 3;
590 scalar det = v2*u1 - u2*
v1;
592 scalar u0 = p[i1] - a_[i1];
593 scalar v0 = p[i2] - a_[i2];
604 alpha = (v0 - beta*
v2)/v1;
606 hit = ((alpha >= 0) && ((alpha + beta) <= 1));
610 beta = (v0*u1 - u0*
v1)/det;
612 alpha = (u0 - beta*u2)/u1;
614 hit = ((alpha >= 0) && ((alpha + beta) <= 1));
630 alpha =
max(0.0,
min(1.0, alpha));
631 beta =
max(0.0,
min(1.0, beta));
651 else if ((alpha >= 0) && (alpha <= 1) && (beta >= 0) && (beta <= 1))
671 else if ((beta >= 0) && (beta <= 1))
691 else if ((alpha >= 0) && (alpha <= 1))
707 template<
class po
int,
class po
intRef>
713 is >> t.a_ >> t.b_ >> t.c_;
719 is.
check(
"Istream& operator>>(Istream&, triangle&)");
725 template<
class Po
int,
class Po
intRef>
726 inline Ostream& operator<<(Ostream& os, const triangle<Point, PointRef>& t)