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-2024 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>
179 {
180  const vector a = b_ - a_;
181  const vector b = c_ - a_;
182  const vector c = d_ - a_;
183 
184  const vector ba = b ^ a;
185  const vector ca = c ^ a;
186 
187  const scalar lambda = magSqr(c) - (a & c);
188  const scalar mu = magSqr(b) - (a & b);
189 
190  const vector num = lambda*ba - mu*ca;
191  const scalar denom = (c & ba);
192 
193  if (Foam::mag(denom) < rootVSmall)
194  {
195  // Degenerate. Return a point far away.
196  static const scalar sqrt3 = sqrt(scalar(3));
197  return Tuple2<Point, scalar>(point::uniform(great), sqrt3*great);
198  }
199  else
200  {
201  const vector v = (a + num/denom)/2;
202  return Tuple2<Point, scalar>(a_ + v, Foam::mag(v));
203  }
204 }
205 
206 
207 template<class Point, class PointRef>
209 {
210  const scalar r = circumSphere().second();
211  static const scalar sqrt3 = sqrt(scalar(3));
212  return mag()/((8.0/27.0)*sqrt3*pow3(min(r, great)) + rootVSmall);
213 }
214 
215 
216 template<class Point, class PointRef>
218 (
220 ) const
221 {
222  return barycentricToPoint(barycentric01(rndGen));
223 }
224 
225 
226 template<class Point, class PointRef>
228 (
229  const barycentric& bary
230 ) const
231 {
232  return bary[0]*a_ + bary[1]*b_ + bary[2]*c_ + bary[3]*d_;
233 }
234 
235 
236 template<class Point, class PointRef>
238 (
239  const point& pt
240 ) const
241 {
242  barycentric bary;
243  pointToBarycentric(pt, bary);
244  return bary;
245 }
246 
247 
248 template<class Point, class PointRef>
250 (
251  const point& pt,
252  barycentric& bary
253 ) const
254 {
255  // Reference:
256  // http://en.wikipedia.org/wiki/Barycentric_coordinate_system_(mathematics)
257 
258  vector e0(a_ - d_);
259  vector e1(b_ - d_);
260  vector e2(c_ - d_);
261 
262  tensor t
263  (
264  e0.x(), e1.x(), e2.x(),
265  e0.y(), e1.y(), e2.y(),
266  e0.z(), e1.z(), e2.z()
267  );
268 
269  scalar detT = det(t);
270 
271  if (Foam::mag(detT) < small)
272  {
273  // Degenerate tetrahedron, returning 1/4 barycentric coordinates
274 
275  bary = barycentric(0.25, 0.25, 0.25, 0.25);
276 
277  return detT;
278  }
279 
280  vector res = inv(t, detT) & (pt - d_);
281 
282  bary[0] = res.x();
283  bary[1] = res.y();
284  bary[2] = res.z();
285  bary[3] = 1 - cmptSum(res);
286 
287  return detT;
288 }
289 
290 
291 template<class Point, class PointRef>
293 (
294  const point& p
295 ) const
296 {
297  // Adapted from:
298  // Real-time collision detection, Christer Ericson, 2005, p142-144
299 
300  // Assuming initially that the point is inside all of the
301  // halfspaces, so closest to itself.
302 
303  point closestPt = p;
304 
305  scalar minOutsideDistance = vGreat;
306 
307  bool inside = true;
308 
309  if (((p - b_) & Sa()) >= 0)
310  {
311  // p is outside halfspace plane of tri
312  pointHit info = triangle<Point, PointRef>(b_, c_, d_).nearestPoint(p);
313 
314  inside = false;
315 
316  if (info.distance() < minOutsideDistance)
317  {
318  closestPt = info.rawPoint();
319 
320  minOutsideDistance = info.distance();
321  }
322  }
323 
324  if (((p - a_) & Sb()) >= 0)
325  {
326  // p is outside halfspace plane of tri
327  pointHit info = triangle<Point, PointRef>(a_, d_, c_).nearestPoint(p);
328 
329  inside = false;
330 
331  if (info.distance() < minOutsideDistance)
332  {
333  closestPt = info.rawPoint();
334 
335  minOutsideDistance = info.distance();
336  }
337  }
338 
339  if (((p - a_) & Sc()) >= 0)
340  {
341  // p is outside halfspace plane of tri
342  pointHit info = triangle<Point, PointRef>(a_, b_, d_).nearestPoint(p);
343 
344  inside = false;
345 
346  if (info.distance() < minOutsideDistance)
347  {
348  closestPt = info.rawPoint();
349 
350  minOutsideDistance = info.distance();
351  }
352  }
353 
354  if (((p - a_) & Sd()) >= 0)
355  {
356  // p is outside halfspace plane of tri
357  pointHit info = triangle<Point, PointRef>(a_, c_, b_).nearestPoint(p);
358 
359  inside = false;
360 
361  if (info.distance() < minOutsideDistance)
362  {
363  closestPt = info.rawPoint();
364 
365  minOutsideDistance = info.distance();
366  }
367  }
368 
369  // If the point is inside, then the distance to the closest point
370  // is zero
371  if (inside)
372  {
373  minOutsideDistance = 0;
374  }
375 
376  return pointHit
377  (
378  inside,
379  closestPt,
380  minOutsideDistance,
381  !inside
382  );
383 }
384 
385 
386 template<class Point, class PointRef>
388 {
389  // For robustness, assuming that the point is in the tet unless
390  // "definitively" shown otherwise by obtaining a positive dot
391  // product greater than a tolerance of small.
392 
393  // The tet is defined: tet(Cc, tetBasePt, pA, pB) where the area
394  // vectors and base points for the half-space planes are:
395  // area[0] = Sa();
396  // area[1] = Sb();
397  // area[2] = Sc();
398  // area[3] = Sd();
399  // planeBase[0] = tetBasePt = b_
400  // planeBase[1] = ptA = c_
401  // planeBase[2] = tetBasePt = b_
402  // planeBase[3] = tetBasePt = b_
403 
404  vector n = Zero;
405 
406  {
407  // 0, a
408  const point& basePt = b_;
409 
410  n = Sa();
411  n /= (Foam::mag(n) + vSmall);
412 
413  if (((pt - basePt) & n) > small)
414  {
415  return false;
416  }
417  }
418 
419  {
420  // 1, b
421  const point& basePt = c_;
422 
423  n = Sb();
424  n /= (Foam::mag(n) + vSmall);
425 
426  if (((pt - basePt) & n) > small)
427  {
428  return false;
429  }
430  }
431 
432  {
433  // 2, c
434  const point& basePt = b_;
435 
436  n = Sc();
437  n /= (Foam::mag(n) + vSmall);
438 
439  if (((pt - basePt) & n) > small)
440  {
441  return false;
442  }
443  }
444 
445  {
446  // 3, d
447  const point& basePt = b_;
448 
449  n = Sd();
450  n /= (Foam::mag(n) + vSmall);
451 
452  if (((pt - basePt) & n) > small)
453  {
454  return false;
455  }
456  }
457 
458  return true;
459 }
460 
461 
462 template<class Point, class PointRef>
464 {
465  return
466  boundBox
467  (
468  min(a(), min(b(), min(c(), d()))),
469  max(a(), max(b(), max(c(), d())))
470  );
471 }
472 
473 
474 // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
475 
476 template<class Point, class PointRef>
477 inline Foam::Istream& Foam::operator>>
478 (
479  Istream& is,
481 )
482 {
483  is.readBegin("tetrahedron");
484  is >> t.a_ >> t.b_ >> t.c_ >> t.d_;
485  is.readEnd("tetrahedron");
486 
487  is.check("Istream& operator>>(Istream&, tetrahedron&)");
488 
489  return is;
490 }
491 
492 
493 template<class Point, class PointRef>
494 inline Foam::Ostream& Foam::operator<<
495 (
496  Ostream& os,
497  const tetrahedron<Point, PointRef>& t
498 )
499 {
500  os << nl
502  << t.a_ << token::SPACE
503  << t.b_ << token::SPACE
504  << t.c_ << token::SPACE
505  << t.d_
506  << token::END_LIST;
507 
508  return os;
509 }
510 
511 
512 // ************************************************************************* //
Useful combination of include files which define Sin, Sout and Serr and the use of IO streams general...
label n
virtual bool check(const char *operation) const
Check IOstream status for given operation.
Definition: IOstream.C:92
An Istream is an abstract base class for all input systems (streams, files, token lists etc)....
Definition: Istream.H:60
Istream & readEnd(const char *funcName)
Definition: Istream.C:103
Istream & readBegin(const char *funcName)
Definition: Istream.C:86
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:57
scalar distance() const
Return distance to hit.
Definition: PointHit.H:139
const Point & rawPoint() const
Return point with no checking.
Definition: PointHit.H:158
A 2-tuple for storing two objects of different types.
Definition: Tuple2.H:66
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: UList.H:74
static Form uniform(const Cmpt &s)
Return a VectorSpace with all elements = s.
Definition: VectorSpaceI.H:168
const Cmpt & z() const
Definition: VectorI.H:87
const Cmpt & y() const
Definition: VectorI.H:81
const Cmpt & x() const
Definition: VectorI.H:75
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:59
Random number generator.
A tetrahedron primitive.
Definition: tetrahedron.H:82
barycentric pointToBarycentric(const point &pt) const
Calculate the barycentric coordinates from the given point.
Definition: tetrahedronI.H:238
vector Sd() const
Definition: tetrahedronI.H:156
triPointRef tri(const label facei) const
Return i-th face.
Definition: tetrahedronI.H:101
boundBox bounds() const
Calculate the bounding box.
Definition: tetrahedronI.H:463
const Point & a() const
Return vertices.
Definition: tetrahedronI.H:72
Point centre() const
Return centre (centroid)
Definition: tetrahedronI.H:163
vector Sc() const
Definition: tetrahedronI.H:149
bool inside(const point &pt) const
Return true if point is inside tetrahedron.
Definition: tetrahedronI.H:387
vector Sb() const
Definition: tetrahedronI.H:142
Tuple2< Point, scalar > circumSphere() const
Return the circum centre and radius.
Definition: tetrahedronI.H:178
const Point & c() const
Definition: tetrahedronI.H:86
scalar mag() const
Return volume.
Definition: tetrahedronI.H:170
pointHit nearestPoint(const point &p) const
Return nearest point to p on tetrahedron. Is p itself.
Definition: tetrahedronI.H:293
Point randomPoint(randomGenerator &rndGen) const
Return a random point in the tetrahedron from a.
Definition: tetrahedronI.H:218
Point barycentricToPoint(const barycentric &bary) const
Calculate the point from the given barycentric coordinates.
Definition: tetrahedronI.H:228
tetrahedron(const Point &a, const Point &b, const Point &c, const Point &d)
Construct from points.
Definition: tetrahedronI.H:34
const Point & b() const
Definition: tetrahedronI.H:79
scalar quality() const
Return quality: Ratio of tetrahedron and circum-sphere.
Definition: tetrahedronI.H:208
const Point & d() const
Definition: tetrahedronI.H:93
vector Sa() const
Return face normal.
Definition: tetrahedronI.H:135
@ BEGIN_LIST
Definition: token.H:109
@ END_LIST
Definition: token.H:110
A triangle primitive used to calculate face areas and swept volumes.
Definition: triangle.H:80
vector area() const
Return vector area.
Definition: triangleI.H:96
pointHit nearestPoint(const point &p) const
Return nearest point to p on triangle.
Definition: triangleI.H:645
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
const pointField & points
volScalarField & b
Definition: createFields.H:25
dimensionedScalar lambda(viscosity->lookup("lambda"))
const dimensionedScalar mu
Atomic mass unit.
const dimensionedScalar c
Speed of light in a vacuum.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
static const zero Zero
Definition: zero.H:97
Barycentric< scalar > barycentric
A scalar version of the templated Barycentric.
Definition: barycentric.H:45
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
dimensionedScalar pow3(const dimensionedScalar &ds)
errorManip< error > abort(error &err)
Definition: errorManip.H:131
barycentric barycentric01(randomGenerator &rndGen)
Generate a random barycentric coordinate within the unit tetrahedron.
Definition: barycentric.C:31
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
PointHit< point > pointHit
Definition: pointHit.H:41
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
error FatalError
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
Cmpt cmptSum(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:507
static const char nl
Definition: Ostream.H:266
dimensioned< scalar > magSqr(const dimensioned< Type > &)
randomGenerator rndGen(653213)
volScalarField & p