All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
tetrahedronI.H
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
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.
15 
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
19  for more details.
20 
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/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "triangle.H"
27 #include "IOstreams.H"
28 #include "plane.H"
29 
30 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
31 
32 template<class Point, class PointRef>
34 (
35  const Point& a,
36  const Point& b,
37  const Point& c,
38  const Point& d
39 )
40 :
41  a_(a),
42  b_(b),
43  c_(c),
44  d_(d)
45 {}
46 
47 
48 template<class Point, class PointRef>
50 (
51  const UList<Point>& points,
52  const FixedList<label, 4>& indices
53 )
54 :
55  a_(points[indices[0]]),
56  b_(points[indices[1]]),
57  c_(points[indices[2]]),
58  d_(points[indices[3]])
59 {}
60 
61 
62 template<class Point, class PointRef>
64 {
65  is >> *this;
66 }
67 
68 
69 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
70 
71 template<class Point, class PointRef>
72 inline const Point& Foam::tetrahedron<Point, PointRef>::a() const
73 {
74  return a_;
75 }
76 
77 
78 template<class Point, class PointRef>
79 inline const Point& Foam::tetrahedron<Point, PointRef>::b() const
80 {
81  return b_;
82 }
83 
84 
85 template<class Point, class PointRef>
86 inline const Point& Foam::tetrahedron<Point, PointRef>::c() const
87 {
88  return c_;
89 }
90 
91 
92 template<class Point, class PointRef>
93 inline const Point& Foam::tetrahedron<Point, PointRef>::d() const
94 {
95  return d_;
96 }
97 
98 
99 template<class Point, class PointRef>
101 (
102  const label facei
103 ) const
104 {
105  // Warning. Ordering of faces needs to be the same for a tetrahedron
106  // class, a tetrahedron cell shape model and a tetCell
107 
108  if (facei == 0)
109  {
110  return triPointRef(b_, c_, d_);
111  }
112  else if (facei == 1)
113  {
114  return triPointRef(a_, d_, c_);
115  }
116  else if (facei == 2)
117  {
118  return triPointRef(a_, b_, d_);
119  }
120  else if (facei == 3)
121  {
122  return triPointRef(a_, c_, b_);
123  }
124  else
125  {
127  << "index out of range 0 -> 3. facei = " << facei
128  << abort(FatalError);
129  return triPointRef(b_, c_, d_);
130  }
131 }
132 
133 
134 template<class Point, class PointRef>
136 {
137  return triangle<Point, PointRef>(b_, c_, d_).area();
138 }
139 
140 
141 template<class Point, class PointRef>
143 {
144  return triangle<Point, PointRef>(a_, d_, c_).area();
145 }
146 
147 
148 template<class Point, class PointRef>
150 {
151  return triangle<Point, PointRef>(a_, b_, d_).area();
152 }
153 
154 
155 template<class Point, class PointRef>
157 {
158  return triangle<Point, PointRef>(a_, c_, b_).area();
159 }
160 
161 
162 template<class Point, class PointRef>
164 {
165  return 0.25*(a_ + b_ + c_ + d_);
166 }
167 
168 
169 template<class Point, class PointRef>
170 inline Foam::scalar Foam::tetrahedron<Point, PointRef>::mag() const
171 {
172  return (1.0/6.0)*(((b_ - a_) ^ (c_ - a_)) & (d_ - a_));
173 }
174 
175 
176 template<class Point, class PointRef>
178 {
179  vector a = b_ - a_;
180  vector b = c_ - a_;
181  vector c = d_ - a_;
182 
183  scalar lambda = magSqr(c) - (a & c);
184  scalar mu = magSqr(b) - (a & b);
185 
186  vector ba = b ^ a;
187  vector ca = c ^ a;
188 
189  vector num = lambda*ba - mu*ca;
190  scalar denom = (c & ba);
191 
192  if (Foam::mag(denom) < rootVSmall)
193  {
194  // Degenerate tetrahedron, returning centre instead of circumCentre.
195 
196  return centre();
197  }
198 
199  return a_ + 0.5*(a + num/denom);
200 }
201 
202 
203 template<class Point, class PointRef>
205 {
206  vector a = b_ - a_;
207  vector b = c_ - a_;
208  vector c = d_ - a_;
209 
210  scalar lambda = magSqr(c) - (a & c);
211  scalar mu = magSqr(b) - (a & b);
212 
213  vector ba = b ^ a;
214  vector ca = c ^ a;
215 
216  vector num = lambda*ba - mu*ca;
217  scalar denom = (c & ba);
218 
219  if (Foam::mag(denom) < rootVSmall)
220  {
221  // Degenerate tetrahedron, returning great for circumRadius.
222  return great;
223  }
224 
225  return Foam::mag(0.5*(a + num/denom));
226 }
227 
228 
229 template<class Point, class PointRef>
231 {
232  return
233  mag()
234  /(
235  8.0/(9.0*sqrt(3.0))
236  *pow3(min(circumRadius(), great))
237  + rootVSmall
238  );
239 }
240 
241 
242 template<class Point, class PointRef>
244 (
245  Random& rndGen
246 ) const
247 {
248  return barycentricToPoint(barycentric01(rndGen));
249 }
250 
251 
252 template<class Point, class PointRef>
254 (
255  const barycentric& bary
256 ) const
257 {
258  return bary[0]*a_ + bary[1]*b_ + bary[2]*c_ + bary[3]*d_;
259 }
260 
261 
262 template<class Point, class PointRef>
264 (
265  const point& pt
266 ) const
267 {
268  barycentric bary;
269  pointToBarycentric(pt, bary);
270  return bary;
271 }
272 
273 
274 template<class Point, class PointRef>
276 (
277  const point& pt,
278  barycentric& bary
279 ) const
280 {
281  // Reference:
282  // http://en.wikipedia.org/wiki/Barycentric_coordinate_system_(mathematics)
283 
284  vector e0(a_ - d_);
285  vector e1(b_ - d_);
286  vector e2(c_ - d_);
287 
288  tensor t
289  (
290  e0.x(), e1.x(), e2.x(),
291  e0.y(), e1.y(), e2.y(),
292  e0.z(), e1.z(), e2.z()
293  );
294 
295  scalar detT = det(t);
296 
297  if (Foam::mag(detT) < small)
298  {
299  // Degenerate tetrahedron, returning 1/4 barycentric coordinates
300 
301  bary = barycentric(0.25, 0.25, 0.25, 0.25);
302 
303  return detT;
304  }
305 
306  vector res = inv(t, detT) & (pt - d_);
307 
308  bary[0] = res.x();
309  bary[1] = res.y();
310  bary[2] = res.z();
311  bary[3] = 1 - cmptSum(res);
312 
313  return detT;
314 }
315 
316 
317 template<class Point, class PointRef>
319 (
320  const point& p
321 ) const
322 {
323  // Adapted from:
324  // Real-time collision detection, Christer Ericson, 2005, p142-144
325 
326  // Assuming initially that the point is inside all of the
327  // halfspaces, so closest to itself.
328 
329  point closestPt = p;
330 
331  scalar minOutsideDistance = vGreat;
332 
333  bool inside = true;
334 
335  if (((p - b_) & Sa()) >= 0)
336  {
337  // p is outside halfspace plane of tri
338  pointHit info = triangle<Point, PointRef>(b_, c_, d_).nearestPoint(p);
339 
340  inside = false;
341 
342  if (info.distance() < minOutsideDistance)
343  {
344  closestPt = info.rawPoint();
345 
346  minOutsideDistance = info.distance();
347  }
348  }
349 
350  if (((p - a_) & Sb()) >= 0)
351  {
352  // p is outside halfspace plane of tri
353  pointHit info = triangle<Point, PointRef>(a_, d_, c_).nearestPoint(p);
354 
355  inside = false;
356 
357  if (info.distance() < minOutsideDistance)
358  {
359  closestPt = info.rawPoint();
360 
361  minOutsideDistance = info.distance();
362  }
363  }
364 
365  if (((p - a_) & Sc()) >= 0)
366  {
367  // p is outside halfspace plane of tri
368  pointHit info = triangle<Point, PointRef>(a_, b_, d_).nearestPoint(p);
369 
370  inside = false;
371 
372  if (info.distance() < minOutsideDistance)
373  {
374  closestPt = info.rawPoint();
375 
376  minOutsideDistance = info.distance();
377  }
378  }
379 
380  if (((p - a_) & Sd()) >= 0)
381  {
382  // p is outside halfspace plane of tri
383  pointHit info = triangle<Point, PointRef>(a_, c_, b_).nearestPoint(p);
384 
385  inside = false;
386 
387  if (info.distance() < minOutsideDistance)
388  {
389  closestPt = info.rawPoint();
390 
391  minOutsideDistance = info.distance();
392  }
393  }
394 
395  // If the point is inside, then the distance to the closest point
396  // is zero
397  if (inside)
398  {
399  minOutsideDistance = 0;
400  }
401 
402  return pointHit
403  (
404  inside,
405  closestPt,
406  minOutsideDistance,
407  !inside
408  );
409 }
410 
411 
412 template<class Point, class PointRef>
414 {
415  // For robustness, assuming that the point is in the tet unless
416  // "definitively" shown otherwise by obtaining a positive dot
417  // product greater than a tolerance of small.
418 
419  // The tet is defined: tet(Cc, tetBasePt, pA, pB) where the area
420  // vectors and base points for the half-space planes are:
421  // area[0] = Sa();
422  // area[1] = Sb();
423  // area[2] = Sc();
424  // area[3] = Sd();
425  // planeBase[0] = tetBasePt = b_
426  // planeBase[1] = ptA = c_
427  // planeBase[2] = tetBasePt = b_
428  // planeBase[3] = tetBasePt = b_
429 
430  vector n = Zero;
431 
432  {
433  // 0, a
434  const point& basePt = b_;
435 
436  n = Sa();
437  n /= (Foam::mag(n) + vSmall);
438 
439  if (((pt - basePt) & n) > small)
440  {
441  return false;
442  }
443  }
444 
445  {
446  // 1, b
447  const point& basePt = c_;
448 
449  n = Sb();
450  n /= (Foam::mag(n) + vSmall);
451 
452  if (((pt - basePt) & n) > small)
453  {
454  return false;
455  }
456  }
457 
458  {
459  // 2, c
460  const point& basePt = b_;
461 
462  n = Sc();
463  n /= (Foam::mag(n) + vSmall);
464 
465  if (((pt - basePt) & n) > small)
466  {
467  return false;
468  }
469  }
470 
471  {
472  // 3, d
473  const point& basePt = b_;
474 
475  n = Sd();
476  n /= (Foam::mag(n) + vSmall);
477 
478  if (((pt - basePt) & n) > small)
479  {
480  return false;
481  }
482  }
483 
484  return true;
485 }
486 
487 
488 // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
489 
490 template<class Point, class PointRef>
491 inline Foam::Istream& Foam::operator>>
492 (
493  Istream& is,
495 )
496 {
497  is.readBegin("tetrahedron");
498  is >> t.a_ >> t.b_ >> t.c_ >> t.d_;
499  is.readEnd("tetrahedron");
500 
501  is.check("Istream& operator>>(Istream&, tetrahedron&)");
502 
503  return is;
504 }
505 
506 
507 template<class Point, class PointRef>
508 inline Foam::Ostream& Foam::operator<<
509 (
510  Ostream& os,
512 )
513 {
514  os << nl
515  << token::BEGIN_LIST
516  << t.a_ << token::SPACE
517  << t.b_ << token::SPACE
518  << t.c_ << token::SPACE
519  << t.d_
520  << token::END_LIST;
521 
522  return os;
523 }
524 
525 
526 // ************************************************************************* //
Istream & readBegin(const char *funcName)
Definition: Istream.C:86
A tetrahedron primitive.
Definition: tetrahedron.H:61
A triangle primitive used to calculate face areas and swept volumes.
Definition: triangle.H:57
barycentric barycentric01(Random &rndGen)
Generate a random barycentric coordinate within the unit tetrahedron.
Definition: barycentric.C:31
error FatalError
scalar circumRadius() const
Return circum-radius.
Definition: tetrahedronI.H:204
tetrahedron(const Point &a, const Point &b, const Point &c, const Point &d)
Construct from points.
Definition: tetrahedronI.H:34
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
virtual bool check(const char *operation) const
Check IOstream status for given operation.
Definition: IOstream.C:92
Point circumCentre() const
Return circum-centre.
Definition: tetrahedronI.H:177
const Point & c() const
Definition: tetrahedronI.H:86
vector Sd() const
Definition: tetrahedronI.H:156
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
Definition: Istream.H:57
const Point & a() const
Return vertices.
Definition: tetrahedronI.H:72
Cmpt cmptSum(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:500
dimensionedScalar sqrt(const dimensionedScalar &ds)
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
dimensionedScalar lambda(viscosity->lookup("lambda"))
const Cmpt & z() const
Definition: VectorI.H:87
Barycentric< scalar > barycentric
A scalar version of the templated Barycentric.
Definition: barycentric.H:45
dimensionedScalar det(const dimensionedSphericalTensor &dt)
const dimensionedScalar c
Speed of light in a vacuum.
const Cmpt & y() const
Definition: VectorI.H:81
Point barycentricToPoint(const barycentric &bary) const
Calculate the point from the given barycentric coordinates.
Definition: tetrahedronI.H:254
Useful combination of include files which define Sin, Sout and Serr and the use of IO streams general...
const Point & d() const
Definition: tetrahedronI.H:93
Istream & readEnd(const char *funcName)
Definition: Istream.C:103
const Point & b() const
Definition: tetrahedronI.H:79
scalar quality() const
Return quality: Ratio of tetrahedron and circum-sphere.
Definition: tetrahedronI.H:230
vector Sb() const
Definition: tetrahedronI.H:142
scalar mag() const
Return volume.
Definition: tetrahedronI.H:170
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
vector Sa() const
Return face normal.
Definition: tetrahedronI.H:135
static const zero Zero
Definition: zero.H:97
Point centre() const
Return centre (centroid)
Definition: tetrahedronI.H:163
bool inside(const point &pt) const
Return true if point is inside tetrahedron.
Definition: tetrahedronI.H:413
errorManip< error > abort(error &err)
Definition: errorManip.H:131
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
const Cmpt & x() const
Definition: VectorI.H:75
Random number generator.
Definition: Random.H:57
dimensioned< scalar > magSqr(const dimensioned< Type > &)
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:54
const Point & rawPoint() const
Return point with no checking.
Definition: PointHit.H:158
static const char nl
Definition: Ostream.H:260
const dimensionedScalar mu
Atomic mass unit.
triPointRef tri(const label facei) const
Return i-th face.
Definition: tetrahedronI.H:101
dimensionedScalar pow3(const dimensionedScalar &ds)
barycentric pointToBarycentric(const point &pt) const
Calculate the barycentric coordinates from the given point.
Definition: tetrahedronI.H:264
pointHit nearestPoint(const point &p) const
Return nearest point to p on tetrahedron. Is p itself.
Definition: tetrahedronI.H:319
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
PointHit< point > pointHit
Definition: pointHit.H:41
Point randomPoint(Random &rndGen) const
Return a random point in the tetrahedron from a.
Definition: tetrahedronI.H:244
scalar distance() const
Return distance to hit.
Definition: PointHit.H:139
volScalarField & p
vector Sc() const
Definition: tetrahedronI.H:149