tetrahedronI.H
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2013 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 "tetPoints.H"
29 #include "plane.H"
30 
31 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
32 
33 template<class Point, class PointRef>
35 (
36  const Point& a,
37  const Point& b,
38  const Point& c,
39  const Point& d
40 )
41 :
42  a_(a),
43  b_(b),
44  c_(c),
45  d_(d)
46 {}
47 
48 
49 template<class Point, class PointRef>
51 (
52  const UList<Point>& points,
53  const FixedList<label, 4>& indices
54 )
55 :
56  a_(points[indices[0]]),
57  b_(points[indices[1]]),
58  c_(points[indices[2]]),
59  d_(points[indices[3]])
60 {}
61 
62 
63 template<class Point, class PointRef>
65 {
66  is >> *this;
67 }
68 
69 
70 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
71 
72 template<class Point, class PointRef>
73 inline const Point& Foam::tetrahedron<Point, PointRef>::a() const
74 {
75  return a_;
76 }
77 
78 
79 template<class Point, class PointRef>
80 inline const Point& Foam::tetrahedron<Point, PointRef>::b() const
81 {
82  return b_;
83 }
84 
85 
86 template<class Point, class PointRef>
87 inline const Point& Foam::tetrahedron<Point, PointRef>::c() const
88 {
89  return c_;
90 }
91 
92 
93 template<class Point, class PointRef>
94 inline const Point& Foam::tetrahedron<Point, PointRef>::d() const
95 {
96  return d_;
97 }
98 
99 
100 template<class Point, class PointRef>
102 (
103  const label faceI
104 ) const
105 {
106  // Warning. Ordering of faces needs to be the same for a tetrahedron
107  // class, a tetrahedron cell shape model and a tetCell
108 
109  if (faceI == 0)
110  {
111  return triPointRef(b_, c_, d_);
112  }
113  else if (faceI == 1)
114  {
115  return triPointRef(a_, d_, c_);
116  }
117  else if (faceI == 2)
118  {
119  return triPointRef(a_, b_, d_);
120  }
121  else if (faceI == 3)
122  {
123  return triPointRef(a_, c_, b_);
124  }
125  else
126  {
127  FatalErrorIn("tetrahedron::tri(const label faceI) const")
128  << "index out of range 0 -> 3. faceI = " << faceI
129  << abort(FatalError);
130  return triPointRef(b_, c_, d_);
131  }
132 }
133 
134 
135 template<class Point, class PointRef>
137 {
138  return triangle<Point, PointRef>(b_, c_, d_).normal();
139 }
140 
141 
142 template<class Point, class PointRef>
144 {
145  return triangle<Point, PointRef>(a_, d_, c_).normal();
146 }
147 
148 
149 template<class Point, class PointRef>
151 {
152  return triangle<Point, PointRef>(a_, b_, d_).normal();
153 }
154 
155 
156 template<class Point, class PointRef>
158 {
159  return triangle<Point, PointRef>(a_, c_, b_).normal();
160 }
161 
162 
163 template<class Point, class PointRef>
165 {
166  return 0.25*(a_ + b_ + c_ + d_);
167 }
168 
169 
170 template<class Point, class PointRef>
171 inline Foam::scalar Foam::tetrahedron<Point, PointRef>::mag() const
172 {
173  return (1.0/6.0)*(((b_ - a_) ^ (c_ - a_)) & (d_ - a_));
174 }
175 
176 
177 template<class Point, class PointRef>
179 {
180  vector a = b_ - a_;
181  vector b = c_ - a_;
182  vector c = d_ - a_;
183 
184  scalar lambda = magSqr(c) - (a & c);
185  scalar mu = magSqr(b) - (a & b);
186 
187  vector ba = b ^ a;
188  vector ca = c ^ a;
189 
190  vector num = lambda*ba - mu*ca;
191  scalar denom = (c & ba);
192 
193  if (Foam::mag(denom) < ROOTVSMALL)
194  {
195  // Degenerate tetrahedron, returning centre instead of circumCentre.
196 
197  return centre();
198  }
199 
200  return a_ + 0.5*(a + num/denom);
201 }
202 
203 
204 template<class Point, class PointRef>
206 {
207  vector a = b_ - a_;
208  vector b = c_ - a_;
209  vector c = d_ - a_;
210 
211  scalar lambda = magSqr(c) - (a & c);
212  scalar mu = magSqr(b) - (a & b);
213 
214  vector ba = b ^ a;
215  vector ca = c ^ a;
216 
217  vector num = lambda*ba - mu*ca;
218  scalar denom = (c & ba);
219 
220  if (Foam::mag(denom) < ROOTVSMALL)
221  {
222  // Degenerate tetrahedron, returning GREAT for circumRadius.
223  return GREAT;
224  }
225 
226  return Foam::mag(0.5*(a + num/denom));
227 }
228 
229 
230 template<class Point, class PointRef>
232 {
233  return
234  mag()
235  /(
236  8.0/(9.0*sqrt(3.0))
237  *pow3(min(circumRadius(), GREAT))
238  + ROOTVSMALL
239  );
240 }
241 
242 
243 template<class Point, class PointRef>
245 (
246  Random& rndGen
247 ) const
248 {
249  // Adapted from
250  // http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html
251 
252  scalar s = rndGen.scalar01();
253  scalar t = rndGen.scalar01();
254  scalar u = rndGen.scalar01();
255 
256  if (s + t > 1.0)
257  {
258  s = 1.0 - s;
259  t = 1.0 - t;
260  }
261 
262  if (t + u > 1.0)
263  {
264  scalar tmp = u;
265  u = 1.0 - s - t;
266  t = 1.0 - tmp;
267  }
268  else if (s + t + u > 1.0)
269  {
270  scalar tmp = u;
271  u = s + t + u - 1.0;
272  s = 1.0 - t - tmp;
273  }
274 
275  return (1 - s - t - u)*a_ + s*b_ + t*c_ + u*d_;
276 }
277 
278 
279 template<class Point, class PointRef>
281 (
282  cachedRandom& rndGen
283 ) const
284 {
285  // Adapted from
286  // http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html
287 
288  scalar s = rndGen.sample01<scalar>();
289  scalar t = rndGen.sample01<scalar>();
290  scalar u = rndGen.sample01<scalar>();
291 
292  if (s + t > 1.0)
293  {
294  s = 1.0 - s;
295  t = 1.0 - t;
296  }
297 
298  if (t + u > 1.0)
299  {
300  scalar tmp = u;
301  u = 1.0 - s - t;
302  t = 1.0 - tmp;
303  }
304  else if (s + t + u > 1.0)
305  {
306  scalar tmp = u;
307  u = s + t + u - 1.0;
308  s = 1.0 - t - tmp;
309  }
310 
311  return (1 - s - t - u)*a_ + s*b_ + t*c_ + u*d_;
312 }
313 
314 
315 template<class Point, class PointRef>
317 (
318  const point& pt,
319  List<scalar>& bary
320 ) const
321 {
322  // From:
323  // http://en.wikipedia.org/wiki/Barycentric_coordinate_system_(mathematics)
324 
325  vector e0(a_ - d_);
326  vector e1(b_ - d_);
327  vector e2(c_ - d_);
328 
329  tensor t
330  (
331  e0.x(), e1.x(), e2.x(),
332  e0.y(), e1.y(), e2.y(),
333  e0.z(), e1.z(), e2.z()
334  );
335 
336  scalar detT = det(t);
337 
338  if (Foam::mag(detT) < SMALL)
339  {
340  // Degenerate tetrahedron, returning 1/4 barycentric coordinates.
341 
342  bary = List<scalar>(4, 0.25);
343 
344  return detT;
345  }
346 
347  vector res = inv(t, detT) & (pt - d_);
348 
349  bary.setSize(4);
350 
351  bary[0] = res.x();
352  bary[1] = res.y();
353  bary[2] = res.z();
354  bary[3] = (1.0 - res.x() - res.y() - res.z());
355 
356  return detT;
357 }
358 
359 
360 template<class Point, class PointRef>
362 (
363  const point& p
364 ) const
365 {
366  // Adapted from:
367  // Real-time collision detection, Christer Ericson, 2005, p142-144
368 
369  // Assuming initially that the point is inside all of the
370  // halfspaces, so closest to itself.
371 
372  point closestPt = p;
373 
374  scalar minOutsideDistance = VGREAT;
375 
376  bool inside = true;
377 
378  if (((p - b_) & Sa()) >= 0)
379  {
380  // p is outside halfspace plane of tri
381  pointHit info = triangle<Point, PointRef>(b_, c_, d_).nearestPoint(p);
382 
383  inside = false;
384 
385  if (info.distance() < minOutsideDistance)
386  {
387  closestPt = info.rawPoint();
388 
389  minOutsideDistance = info.distance();
390  }
391  }
392 
393  if (((p - a_) & Sb()) >= 0)
394  {
395  // p is outside halfspace plane of tri
396  pointHit info = triangle<Point, PointRef>(a_, d_, c_).nearestPoint(p);
397 
398  inside = false;
399 
400  if (info.distance() < minOutsideDistance)
401  {
402  closestPt = info.rawPoint();
403 
404  minOutsideDistance = info.distance();
405  }
406  }
407 
408  if (((p - a_) & Sc()) >= 0)
409  {
410  // p is outside halfspace plane of tri
411  pointHit info = triangle<Point, PointRef>(a_, b_, d_).nearestPoint(p);
412 
413  inside = false;
414 
415  if (info.distance() < minOutsideDistance)
416  {
417  closestPt = info.rawPoint();
418 
419  minOutsideDistance = info.distance();
420  }
421  }
422 
423  if (((p - a_) & Sd()) >= 0)
424  {
425  // p is outside halfspace plane of tri
426  pointHit info = triangle<Point, PointRef>(a_, c_, b_).nearestPoint(p);
427 
428  inside = false;
429 
430  if (info.distance() < minOutsideDistance)
431  {
432  closestPt = info.rawPoint();
433 
434  minOutsideDistance = info.distance();
435  }
436  }
437 
438  // If the point is inside, then the distance to the closest point
439  // is zero
440  if (inside)
441  {
442  minOutsideDistance = 0;
443  }
444 
445  return pointHit
446  (
447  inside,
448  closestPt,
449  minOutsideDistance,
450  !inside
451  );
452 }
453 
454 
455 template<class Point, class PointRef>
457 {
458  // For robustness, assuming that the point is in the tet unless
459  // "definitively" shown otherwise by obtaining a positive dot
460  // product greater than a tolerance of SMALL.
461 
462  // The tet is defined: tet(Cc, tetBasePt, pA, pB) where the normal
463  // vectors and base points for the half-space planes are:
464  // area[0] = Sa();
465  // area[1] = Sb();
466  // area[2] = Sc();
467  // area[3] = Sd();
468  // planeBase[0] = tetBasePt = b_
469  // planeBase[1] = ptA = c_
470  // planeBase[2] = tetBasePt = b_
471  // planeBase[3] = tetBasePt = b_
472 
473  vector n = vector::zero;
474 
475  {
476  // 0, a
477  const point& basePt = b_;
478 
479  n = Sa();
480  n /= (Foam::mag(n) + VSMALL);
481 
482  if (((pt - basePt) & n) > SMALL)
483  {
484  return false;
485  }
486  }
487 
488  {
489  // 1, b
490  const point& basePt = c_;
491 
492  n = Sb();
493  n /= (Foam::mag(n) + VSMALL);
494 
495  if (((pt - basePt) & n) > SMALL)
496  {
497  return false;
498  }
499  }
500 
501  {
502  // 2, c
503  const point& basePt = b_;
504 
505  n = Sc();
506  n /= (Foam::mag(n) + VSMALL);
507 
508  if (((pt - basePt) & n) > SMALL)
509  {
510  return false;
511  }
512  }
513 
514  {
515  // 3, d
516  const point& basePt = b_;
517 
518  n = Sd();
519  n /= (Foam::mag(n) + VSMALL);
520 
521  if (((pt - basePt) & n) > SMALL)
522  {
523  return false;
524  }
525  }
526 
527  return true;
528 }
529 
530 
531 template<class Point, class PointRef>
533 (
534  const tetPoints&
535 )
536 {}
537 
538 
539 template<class Point, class PointRef>
541 :
542  vol_(0.0)
543 {}
544 
545 
546 template<class Point, class PointRef>
548 (
549  const tetPoints& tet
550 )
551 {
552  vol_ += tet.tet().mag();
553 }
554 
555 
556 template<class Point, class PointRef>
558 (
559  tetIntersectionList& tets,
560  label& nTets
561 )
562 :
563  tets_(tets),
564  nTets_(nTets)
565 {}
566 
567 
568 template<class Point, class PointRef>
570 (
571  const tetPoints& tet
572 )
573 {
574  tets_[nTets_++] = tet;
575 }
576 
577 
578 template<class Point, class PointRef>
580 (
581  const FixedList<scalar, 4>& d,
582  const tetPoints& t,
583  const label negI,
584  const label posI
585 )
586 {
587  return
588  (d[posI]*t[negI] - d[negI]*t[posI])
589  / (-d[negI]+d[posI]);
590 }
591 
592 
593 template<class Point, class PointRef>
594 template<class TetOp>
596 (
597  const FixedList<point, 6>& points,
598  TetOp& op
599 )
600 {
601  op(tetPoints(points[1], points[3], points[2], points[0]));
602  op(tetPoints(points[1], points[2], points[3], points[4]));
603  op(tetPoints(points[4], points[2], points[3], points[5]));
604 }
605 
606 
607 template<class Point, class PointRef>
608 template<class AboveTetOp, class BelowTetOp>
611 (
612  const plane& pl,
613  const tetPoints& tet,
614  AboveTetOp& aboveOp,
615  BelowTetOp& belowOp
616 )
617 {
618  // Distance to plane
620  label nPos = 0;
621  forAll(tet, i)
622  {
623  d[i] = ((tet[i]-pl.refPoint()) & pl.normal());
624  if (d[i] > 0)
625  {
626  nPos++;
627  }
628  }
629 
630  if (nPos == 4)
631  {
632  aboveOp(tet);
633  }
634  else if (nPos == 3)
635  {
636  // Sliced into below tet and above prism. Prism gets split into
637  // two tets.
638 
639  // Find the below tet
640  label i0 = -1;
641  forAll(d, i)
642  {
643  if (d[i] <= 0)
644  {
645  i0 = i;
646  break;
647  }
648  }
649 
650  label i1 = d.fcIndex(i0);
651  label i2 = d.fcIndex(i1);
652  label i3 = d.fcIndex(i2);
653 
654  point p01 = planeIntersection(d, tet, i0, i1);
655  point p02 = planeIntersection(d, tet, i0, i2);
656  point p03 = planeIntersection(d, tet, i0, i3);
657 
658  // i0 = tetCell vertex 0: p01,p02,p03 outwards pointing triad
659  // ,, 1 : ,, inwards pointing triad
660  // ,, 2 : ,, outwards pointing triad
661  // ,, 3 : ,, inwards pointing triad
662 
663  //Pout<< "Split 3pos tet " << tet << " d:" << d << " into" << nl;
664 
665  if (i0 == 0 || i0 == 2)
666  {
667  tetPoints t(tet[i0], p01, p02, p03);
668  //Pout<< " belowtet:" << t << " around i0:" << i0 << endl;
669  //checkTet(t, "nPos 3, belowTet i0==0 or 2");
670  belowOp(t);
671 
672  // Prism
674  p[0] = tet[i1];
675  p[1] = tet[i3];
676  p[2] = tet[i2];
677  p[3] = p01;
678  p[4] = p03;
679  p[5] = p02;
680  //Pout<< " aboveprism:" << p << endl;
681  decomposePrism(p, aboveOp);
682  }
683  else
684  {
685  tetPoints t(p01, p02, p03, tet[i0]);
686  //Pout<< " belowtet:" << t << " around i0:" << i0 << endl;
687  //checkTet(t, "nPos 3, belowTet i0==1 or 3");
688  belowOp(t);
689 
690  // Prism
692  p[0] = tet[i3];
693  p[1] = tet[i1];
694  p[2] = tet[i2];
695  p[3] = p03;
696  p[4] = p01;
697  p[5] = p02;
698  //Pout<< " aboveprism:" << p << endl;
699  decomposePrism(p, aboveOp);
700  }
701  }
702  else if (nPos == 2)
703  {
704  // Tet cut into two prisms. Determine the positive one.
705  label pos0 = -1;
706  label pos1 = -1;
707  forAll(d, i)
708  {
709  if (d[i] > 0)
710  {
711  if (pos0 == -1)
712  {
713  pos0 = i;
714  }
715  else
716  {
717  pos1 = i;
718  }
719  }
720  }
721 
722  //Pout<< "Split 2pos tet " << tet << " d:" << d
723  // << " around pos0:" << pos0 << " pos1:" << pos1
724  // << " neg0:" << neg0 << " neg1:" << neg1 << " into" << nl;
725 
726  const edge posEdge(pos0, pos1);
727 
728  if (posEdge == edge(0, 1))
729  {
730  point p02 = planeIntersection(d, tet, 0, 2);
731  point p03 = planeIntersection(d, tet, 0, 3);
732  point p12 = planeIntersection(d, tet, 1, 2);
733  point p13 = planeIntersection(d, tet, 1, 3);
734  // Split the resulting prism
735  {
737  p[0] = tet[0];
738  p[1] = p02;
739  p[2] = p03;
740  p[3] = tet[1];
741  p[4] = p12;
742  p[5] = p13;
743  //Pout<< " 01 aboveprism:" << p << endl;
744  decomposePrism(p, aboveOp);
745  }
746  {
748  p[0] = tet[2];
749  p[1] = p02;
750  p[2] = p12;
751  p[3] = tet[3];
752  p[4] = p03;
753  p[5] = p13;
754  //Pout<< " 01 belowprism:" << p << endl;
755  decomposePrism(p, belowOp);
756  }
757  }
758  else if (posEdge == edge(1, 2))
759  {
760  point p01 = planeIntersection(d, tet, 0, 1);
761  point p13 = planeIntersection(d, tet, 1, 3);
762  point p02 = planeIntersection(d, tet, 0, 2);
763  point p23 = planeIntersection(d, tet, 2, 3);
764  // Split the resulting prism
765  {
767  p[0] = tet[1];
768  p[1] = p01;
769  p[2] = p13;
770  p[3] = tet[2];
771  p[4] = p02;
772  p[5] = p23;
773  //Pout<< " 12 aboveprism:" << p << endl;
774  decomposePrism(p, aboveOp);
775  }
776  {
778  p[0] = tet[3];
779  p[1] = p23;
780  p[2] = p13;
781  p[3] = tet[0];
782  p[4] = p02;
783  p[5] = p01;
784  //Pout<< " 12 belowprism:" << p << endl;
785  decomposePrism(p, belowOp);
786  }
787  }
788  else if (posEdge == edge(2, 0))
789  {
790  point p01 = planeIntersection(d, tet, 0, 1);
791  point p03 = planeIntersection(d, tet, 0, 3);
792  point p12 = planeIntersection(d, tet, 1, 2);
793  point p23 = planeIntersection(d, tet, 2, 3);
794  // Split the resulting prism
795  {
797  p[0] = tet[2];
798  p[1] = p12;
799  p[2] = p23;
800  p[3] = tet[0];
801  p[4] = p01;
802  p[5] = p03;
803  //Pout<< " 20 aboveprism:" << p << endl;
804  decomposePrism(p, aboveOp);
805  }
806  {
808  p[0] = tet[1];
809  p[1] = p12;
810  p[2] = p01;
811  p[3] = tet[3];
812  p[4] = p23;
813  p[5] = p03;
814  //Pout<< " 20 belowprism:" << p << endl;
815  decomposePrism(p, belowOp);
816  }
817  }
818  else if (posEdge == edge(0, 3))
819  {
820  point p01 = planeIntersection(d, tet, 0, 1);
821  point p02 = planeIntersection(d, tet, 0, 2);
822  point p13 = planeIntersection(d, tet, 1, 3);
823  point p23 = planeIntersection(d, tet, 2, 3);
824  // Split the resulting prism
825  {
827  p[0] = tet[3];
828  p[1] = p23;
829  p[2] = p13;
830  p[3] = tet[0];
831  p[4] = p02;
832  p[5] = p01;
833  //Pout<< " 03 aboveprism:" << p << endl;
834  decomposePrism(p, aboveOp);
835  }
836  {
838  p[0] = tet[2];
839  p[1] = p23;
840  p[2] = p02;
841  p[3] = tet[1];
842  p[4] = p13;
843  p[5] = p01;
844  //Pout<< " 03 belowprism:" << p << endl;
845  decomposePrism(p, belowOp);
846  }
847  }
848  else if (posEdge == edge(1, 3))
849  {
850  point p01 = planeIntersection(d, tet, 0, 1);
851  point p12 = planeIntersection(d, tet, 1, 2);
852  point p03 = planeIntersection(d, tet, 0, 3);
853  point p23 = planeIntersection(d, tet, 2, 3);
854  // Split the resulting prism
855  {
857  p[0] = tet[1];
858  p[1] = p12;
859  p[2] = p01;
860  p[3] = tet[3];
861  p[4] = p23;
862  p[5] = p03;
863  //Pout<< " 13 aboveprism:" << p << endl;
864  decomposePrism(p, aboveOp);
865  }
866  {
868  p[0] = tet[2];
869  p[1] = p12;
870  p[2] = p23;
871  p[3] = tet[0];
872  p[4] = p01;
873  p[5] = p03;
874  //Pout<< " 13 belowprism:" << p << endl;
875  decomposePrism(p, belowOp);
876  }
877  }
878  else if (posEdge == edge(2, 3))
879  {
880  point p02 = planeIntersection(d, tet, 0, 2);
881  point p12 = planeIntersection(d, tet, 1, 2);
882  point p03 = planeIntersection(d, tet, 0, 3);
883  point p13 = planeIntersection(d, tet, 1, 3);
884  // Split the resulting prism
885  {
887  p[0] = tet[2];
888  p[1] = p02;
889  p[2] = p12;
890  p[3] = tet[3];
891  p[4] = p03;
892  p[5] = p13;
893  //Pout<< " 23 aboveprism:" << p << endl;
894  decomposePrism(p, aboveOp);
895  }
896  {
898  p[0] = tet[0];
899  p[1] = p02;
900  p[2] = p03;
901  p[3] = tet[1];
902  p[4] = p12;
903  p[5] = p13;
904  //Pout<< " 23 belowprism:" << p << endl;
905  decomposePrism(p, belowOp);
906  }
907  }
908  else
909  {
910  FatalErrorIn("tetSliceWithPlane(..)")
911  << "Missed edge:" << posEdge
912  << abort(FatalError);
913  }
914  }
915  else if (nPos == 1)
916  {
917  // Find the positive tet
918  label i0 = -1;
919  forAll(d, i)
920  {
921  if (d[i] > 0)
922  {
923  i0 = i;
924  break;
925  }
926  }
927 
928  label i1 = d.fcIndex(i0);
929  label i2 = d.fcIndex(i1);
930  label i3 = d.fcIndex(i2);
931 
932  point p01 = planeIntersection(d, tet, i0, i1);
933  point p02 = planeIntersection(d, tet, i0, i2);
934  point p03 = planeIntersection(d, tet, i0, i3);
935 
936  //Pout<< "Split 1pos tet " << tet << " d:" << d << " into" << nl;
937 
938  if (i0 == 0 || i0 == 2)
939  {
940  tetPoints t(tet[i0], p01, p02, p03);
941  //Pout<< " abovetet:" << t << " around i0:" << i0 << endl;
942  //checkTet(t, "nPos 1, aboveTets i0==0 or 2");
943  aboveOp(t);
944 
945  // Prism
947  p[0] = tet[i1];
948  p[1] = tet[i3];
949  p[2] = tet[i2];
950  p[3] = p01;
951  p[4] = p03;
952  p[5] = p02;
953  //Pout<< " belowprism:" << p << endl;
954  decomposePrism(p, belowOp);
955  }
956  else
957  {
958  tetPoints t(p01, p02, p03, tet[i0]);
959  //Pout<< " abovetet:" << t << " around i0:" << i0 << endl;
960  //checkTet(t, "nPos 1, aboveTets i0==1 or 3");
961  aboveOp(t);
962 
963  // Prism
965  p[0] = tet[i3];
966  p[1] = tet[i1];
967  p[2] = tet[i2];
968  p[3] = p03;
969  p[4] = p01;
970  p[5] = p02;
971  //Pout<< " belowprism:" << p << endl;
972  decomposePrism(p, belowOp);
973  }
974  }
975  else // nPos == 0
976  {
977  belowOp(tet);
978  }
979 }
980 
981 
982 template<class Point, class PointRef>
983 template<class AboveTetOp, class BelowTetOp>
985 (
986  const plane& pl,
987  AboveTetOp& aboveOp,
988  BelowTetOp& belowOp
989 ) const
990 {
991  tetSliceWithPlane(pl, tetPoints(a_, b_, c_, d_), aboveOp, belowOp);
992 }
993 
994 
995 // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
996 
997 template<class Point, class PointRef>
998 inline Foam::Istream& Foam::operator>>
999 (
1000  Istream& is,
1002 )
1003 {
1004  is.readBegin("tetrahedron");
1005  is >> t.a_ >> t.b_ >> t.c_ >> t.d_;
1006  is.readEnd("tetrahedron");
1007 
1008  is.check("Istream& operator>>(Istream&, tetrahedron&)");
1009 
1010  return is;
1011 }
1012 
1013 
1014 template<class Point, class PointRef>
1015 inline Foam::Ostream& Foam::operator<<
1016 (
1017  Ostream& os,
1019 )
1020 {
1021  os << nl
1023  << t.a_ << token::SPACE
1024  << t.b_ << token::SPACE
1025  << t.c_ << token::SPACE
1026  << t.d_
1027  << token::END_LIST;
1028 
1029  return os;
1030 }
1031 
1032 
1033 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
const Point & rawPoint() const
Return point with no checking.
Definition: PointHit.H:158
dimensionedScalar pow3(const dimensionedScalar &ds)
Simple random number generator.
Definition: Random.H:49
tetPointRef tet() const
Return the tetrahedron.
Definition: tetPoints.H:80
const dimensionedScalar mu
Atomic mass unit.
Definition: createFields.H:13
const Point & b() const
Definition: tetrahedronI.H:80
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.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE ))
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:28
dimensioned< scalar > mag(const dimensioned< Type > &)
const Point & d() const
Definition: tetrahedronI.H:94
tetrahedron(const Point &a, const Point &b, const Point &c, const Point &d)
Construct from points.
Definition: tetrahedronI.H:35
dimensioned< scalar > magSqr(const dimensioned< Type > &)
triPointRef tri(const label faceI) const
Return i-th face.
Definition: tetrahedronI.H:102
const vector & normal() const
Return plane normal.
Definition: plane.C:262
vector Sc() const
Definition: tetrahedronI.H:150
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
Definition: Istream.H:57
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
scalar barycentric(const point &pt, List< scalar > &bary) const
Calculate the barycentric coordinates of the given.
Definition: tetrahedronI.H:317
bool inside(const point &pt) const
Return true if point is inside tetrahedron.
Definition: tetrahedronI.H:456
void sliceWithPlane(const plane &pl, AboveTetOp &aboveOp, BelowTetOp &belowOp) const
Decompose tet into tets above and below plane.
Definition: tetrahedronI.H:985
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:57
PointHit< point > pointHit
Definition: pointHit.H:41
scalar mag() const
Return volume.
Definition: tetrahedronI.H:171
Point centre() const
Return centre (centroid)
Definition: tetrahedronI.H:164
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
scalar quality() const
Return quality: Ratio of tetrahedron and circum-sphere.
Definition: tetrahedronI.H:231
label n
Useful combination of include files which define Sin, Sout and Serr and the use of IO streams general...
static const char nl
Definition: Ostream.H:260
void setSize(const label)
Reset size of List.
Definition: List.C:318
const Cmpt & y() const
Definition: VectorI.H:71
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
volScalarField & p
Definition: createFields.H:51
vector Sb() const
Definition: tetrahedronI.H:143
dimensionedScalar det(const dimensionedSphericalTensor &dt)
#define forAll(list, i)
Definition: UList.H:421
storeOp(tetIntersectionList &, label &)
Definition: tetrahedronI.H:558
dimensionedScalar lambda(laminarTransport.lookup("lambda"))
const Point & a() const
Return vertices.
Definition: tetrahedronI.H:73
const Cmpt & x() const
Definition: VectorI.H:65
const Point & c() const
Definition: tetrahedronI.H:87
vector Sd() const
Definition: tetrahedronI.H:157
A tetrahedron primitive.
Definition: tetrahedron.H:62
const Cmpt & z() const
Definition: VectorI.H:77
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Istream & readBegin(const char *funcName)
Definition: Istream.C:88
Point circumCentre() const
Return circum-centre.
Definition: tetrahedronI.H:178
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
Definition: plane.H:60
scalar distance() const
Return distance to hit.
Definition: PointHit.H:139
Point randomPoint(Random &rndGen) const
Return a random point in the tetrahedron from a.
Definition: tetrahedronI.H:245
scalar scalar01()
Scalar [0..1] (so including 0,1)
Definition: Random.C:67
error FatalError
Random number generator.
Definition: cachedRandom.H:63
const point & refPoint() const
Return or return plane base point.
Definition: plane.C:269
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
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: FixedListI.H:115
A normal distribution model.
vector Sa() const
Return face normal.
Definition: tetrahedronI.H:136
Type sample01()
Return a sample whose components lie in the range 0-1.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
const dimensionedScalar c
Speed of light in a vacuum.
scalar circumRadius() const
Return circum-radius.
Definition: tetrahedronI.H:205
pointHit nearestPoint(const point &p) const
Return nearest point to p on tetrahedron. Is p itself.
Definition: tetrahedronI.H:362
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
A class for managing temporary objects.
Definition: PtrList.H:118
Tet storage. Null constructable (unfortunately tetrahedron<point, point> is not)
Definition: tetPoints.H:50