1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
11 OpenFOAM is free software: you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by
13 the Free Software Foundation, either version 3 of the License, or
14 (at your option) any later version.
16 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
24 \*---------------------------------------------------------------------------*/
27 #include "IOstreams.H"
28 #include "triPointRef.H"
30 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
32 template<class Point, class PointRef>
33 inline Foam::tetrahedron<Point, PointRef>::tetrahedron
48 template<class Point, class PointRef>
49 inline Foam::tetrahedron<Point, PointRef>::tetrahedron
51 const UList<Point>& points,
52 const FixedList<label, 4>& indices
55 a_(points[indices[0]]),
56 b_(points[indices[1]]),
57 c_(points[indices[2]]),
58 d_(points[indices[3]])
62 template<class Point, class PointRef>
63 inline Foam::tetrahedron<Point, PointRef>::tetrahedron(Istream& is)
69 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
71 template<class Point, class PointRef>
72 inline const Point& Foam::tetrahedron<Point, PointRef>::a() const
78 template<class Point, class PointRef>
79 inline const Point& Foam::tetrahedron<Point, PointRef>::b() const
85 template<class Point, class PointRef>
86 inline const Point& Foam::tetrahedron<Point, PointRef>::c() const
92 template<class Point, class PointRef>
93 inline const Point& Foam::tetrahedron<Point, PointRef>::d() const
99 template<class Point, class PointRef>
100 inline Foam::vector Foam::tetrahedron<Point, PointRef>::Sa() const
102 return triangle<Point, PointRef>(b_, c_, d_).normal();
106 template<class Point, class PointRef>
107 inline Foam::vector Foam::tetrahedron<Point, PointRef>::Sb() const
109 return triangle<Point, PointRef>(a_, d_, c_).normal();
113 template<class Point, class PointRef>
114 inline Foam::vector Foam::tetrahedron<Point, PointRef>::Sc() const
116 return triangle<Point, PointRef>(a_, b_, d_).normal();
120 template<class Point, class PointRef>
121 inline Foam::vector Foam::tetrahedron<Point, PointRef>::Sd() const
123 return triangle<Point, PointRef>(a_, c_, b_).normal();
127 template<class Point, class PointRef>
128 inline Point Foam::tetrahedron<Point, PointRef>::centre() const
130 return 0.25*(a_ + b_ + c_ + d_);
134 template<class Point, class PointRef>
135 inline Foam::scalar Foam::tetrahedron<Point, PointRef>::mag() const
137 return (1.0/6.0)*(((b_ - a_) ^ (c_ - a_)) & (d_ - a_));
141 template<class Point, class PointRef>
142 inline Point Foam::tetrahedron<Point, PointRef>::circumCentre() const
148 scalar lambda = magSqr(c) - (a & c);
149 scalar mu = magSqr(b) - (a & b);
154 vector num = lambda*ba - mu*ca;
155 scalar denom = (c & ba);
157 if (Foam::mag(denom) < ROOTVSMALL)
159 // Degenerate tetrahedron, returning centre instead of circumCentre.
164 return a_ + 0.5*(a + num/denom);
168 template<class Point, class PointRef>
169 inline Foam::scalar Foam::tetrahedron<Point, PointRef>::circumRadius() const
175 scalar lambda = magSqr(c) - (a & c);
176 scalar mu = magSqr(b) - (a & b);
181 vector num = lambda*ba - mu*ca;
182 scalar denom = (c & ba);
184 if (Foam::mag(denom) < ROOTVSMALL)
186 // Degenerate tetrahedron, returning GREAT for circumRadius.
190 return Foam::mag(0.5*(a + num/denom));
194 template<class Point, class PointRef>
195 inline Foam::scalar Foam::tetrahedron<Point, PointRef>::quality() const
201 *pow3(min(circumRadius(), GREAT))
207 template<class Point, class PointRef>
208 inline Point Foam::tetrahedron<Point, PointRef>::randomPoint
214 // http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html
216 scalar s = rndGen.scalar01();
217 scalar t = rndGen.scalar01();
218 scalar u = rndGen.scalar01();
232 else if (s + t + u > 1.0)
239 return (1 - s - t - u)*a_ + s*b_ + t*c_ + u*d_;
243 template<class Point, class PointRef>
244 inline Point Foam::tetrahedron<Point, PointRef>::randomPoint
250 // http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html
252 scalar s = rndGen.sample01<scalar>();
253 scalar t = rndGen.sample01<scalar>();
254 scalar u = rndGen.sample01<scalar>();
268 else if (s + t + u > 1.0)
275 return (1 - s - t - u)*a_ + s*b_ + t*c_ + u*d_;
279 template<class Point, class PointRef>
280 Foam::scalar Foam::tetrahedron<Point, PointRef>::barycentric
287 // http://en.wikipedia.org/wiki/Barycentric_coordinate_system_(mathematics)
295 e0.x(), e1.x(), e2.x(),
296 e0.y(), e1.y(), e2.y(),
297 e0.z(), e1.z(), e2.z()
300 scalar detT = det(t);
302 if (Foam::mag(detT) < SMALL)
304 // Degenerate tetrahedron, returning 1/4 barycentric coordinates.
306 bary = List<scalar>(4, 0.25);
311 vector res = inv(t, detT) & (pt - d_);
318 bary[3] = (1.0 - res.x() - res.y() - res.z());
324 template<class Point, class PointRef>
325 inline Foam::pointHit Foam::tetrahedron<Point, PointRef>::nearestPoint
331 // Real-time collision detection, Christer Ericson, 2005, p142-144
333 // Assuming initially that the point is inside all of the
334 // halfspaces, so closest to itself.
338 scalar minOutsideDistance = VGREAT;
342 if (((p - b_) & Sa()) >= 0)
344 // p is outside halfspace plane of tri
345 pointHit info = triangle<Point, PointRef>(b_, c_, d_).nearestPoint(p);
349 if (info.distance() < minOutsideDistance)
351 closestPt = info.rawPoint();
353 minOutsideDistance = info.distance();
357 if (((p - a_) & Sb()) >= 0)
359 // p is outside halfspace plane of tri
360 pointHit info = triangle<Point, PointRef>(a_, d_, c_).nearestPoint(p);
364 if (info.distance() < minOutsideDistance)
366 closestPt = info.rawPoint();
368 minOutsideDistance = info.distance();
372 if (((p - a_) & Sc()) >= 0)
374 // p is outside halfspace plane of tri
375 pointHit info = triangle<Point, PointRef>(a_, b_, d_).nearestPoint(p);
379 if (info.distance() < minOutsideDistance)
381 closestPt = info.rawPoint();
383 minOutsideDistance = info.distance();
387 if (((p - a_) & Sd()) >= 0)
389 // p is outside halfspace plane of tri
390 pointHit info = triangle<Point, PointRef>(a_, c_, b_).nearestPoint(p);
394 if (info.distance() < minOutsideDistance)
396 closestPt = info.rawPoint();
398 minOutsideDistance = info.distance();
402 // If the point is inside, then the distance to the closest point
406 minOutsideDistance = 0;
419 template<class Point, class PointRef>
420 bool Foam::tetrahedron<Point, PointRef>::inside(const point& pt) const
422 // For robustness, assuming that the point is in the tet unless
423 // "definitively" shown otherwise by obtaining a positive dot
424 // product greater than a tolerance of SMALL.
426 // The tet is defined: tet(Cc, tetBasePt, pA, pB) where the normal
427 // vectors and base points for the half-space planes are:
432 // planeBase[0] = tetBasePt = b_
433 // planeBase[1] = ptA = c_
434 // planeBase[2] = tetBasePt = b_
435 // planeBase[3] = tetBasePt = b_
437 vector n = vector::zero;
441 const point& basePt = b_;
444 n /= (Foam::mag(n) + VSMALL);
446 if (((pt - basePt) & n) > SMALL)
454 const point& basePt = c_;
457 n /= (Foam::mag(n) + VSMALL);
459 if (((pt - basePt) & n) > SMALL)
467 const point& basePt = b_;
470 n /= (Foam::mag(n) + VSMALL);
472 if (((pt - basePt) & n) > SMALL)
480 const point& basePt = b_;
483 n /= (Foam::mag(n) + VSMALL);
485 if (((pt - basePt) & n) > SMALL)
495 // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
497 template<class Point, class PointRef>
498 inline Foam::Istream& Foam::operator>>
501 tetrahedron<Point, PointRef>& t
504 is.readBegin("tetrahedron");
505 is >> t.a_ >> t.b_ >> t.c_ >> t.d_;
506 is.readEnd("tetrahedron");
508 is.check("Istream& operator>>(Istream&, tetrahedron&)");
514 template<class Point, class PointRef>
515 inline Foam::Ostream& Foam::operator<<
518 const tetrahedron<Point, PointRef>& t
523 << t.a_ << token::SPACE
524 << t.b_ << token::SPACE
525 << t.c_ << token::SPACE
533 // ************************************************************************* //