src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
   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 "triPointRef.H"
  29
  30 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  31
  32 template<class Point, class PointRef>
  33 inline Foam::tetrahedron<Point, PointRef>::tetrahedron
  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>
  49 inline Foam::tetrahedron<Point, PointRef>::tetrahedron
  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>
  63 inline Foam::tetrahedron<Point, PointRef>::tetrahedron(Istream& is)
  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>
 100 inline Foam::vector Foam::tetrahedron<Point, PointRef>::Sa() const
 101 {
 102     return triangle<Point, PointRef>(b_, c_, d_).normal();
 103 }
 104
 105
 106 template<class Point, class PointRef>
 107 inline Foam::vector Foam::tetrahedron<Point, PointRef>::Sb() const
 108 {
 109     return triangle<Point, PointRef>(a_, d_, c_).normal();
 110 }
 111
 112
 113 template<class Point, class PointRef>
 114 inline Foam::vector Foam::tetrahedron<Point, PointRef>::Sc() const
 115 {
 116     return triangle<Point, PointRef>(a_, b_, d_).normal();
 117 }
 118
 119
 120 template<class Point, class PointRef>
 121 inline Foam::vector Foam::tetrahedron<Point, PointRef>::Sd() const
 122 {
 123     return triangle<Point, PointRef>(a_, c_, b_).normal();
 124 }
 125
 126
 127 template<class Point, class PointRef>
 128 inline Point Foam::tetrahedron<Point, PointRef>::centre() const
 129 {
 130     return 0.25*(a_ + b_ + c_ + d_);
 131 }
 132
 133
 134 template<class Point, class PointRef>
 135 inline Foam::scalar Foam::tetrahedron<Point, PointRef>::mag() const
 136 {
 137     return (1.0/6.0)*(((b_ - a_) ^ (c_ - a_)) & (d_ - a_));
 138 }
 139
 140
 141 template<class Point, class PointRef>
 142 inline Point Foam::tetrahedron<Point, PointRef>::circumCentre() const
 143 {
 144     vector a = b_ - a_;
 145     vector b = c_ - a_;
 146     vector c = d_ - a_;
 147
 148     scalar lambda = magSqr(c) - (a & c);
 149     scalar mu = magSqr(b) - (a & b);
 150
 151     vector ba = b ^ a;
 152     vector ca = c ^ a;
 153
 154     vector num = lambda*ba - mu*ca;
 155     scalar denom = (c & ba);
 156
 157     if (Foam::mag(denom) < ROOTVSMALL)
 158     {
 159         // Degenerate tetrahedron, returning centre instead of circumCentre.
 160
 161         return centre();
 162     }
 163
 164     return a_ + 0.5*(a + num/denom);
 165 }
 166
 167
 168 template<class Point, class PointRef>
 169 inline Foam::scalar Foam::tetrahedron<Point, PointRef>::circumRadius() const
 170 {
 171     vector a = b_ - a_;
 172     vector b = c_ - a_;
 173     vector c = d_ - a_;
 174
 175     scalar lambda = magSqr(c) - (a & c);
 176     scalar mu = magSqr(b) - (a & b);
 177
 178     vector ba = b ^ a;
 179     vector ca = c ^ a;
 180
 181     vector num = lambda*ba - mu*ca;
 182     scalar denom = (c & ba);
 183
 184     if (Foam::mag(denom) < ROOTVSMALL)
 185     {
 186         // Degenerate tetrahedron, returning GREAT for circumRadius.
 187         return GREAT;
 188     }
 189
 190     return Foam::mag(0.5*(a + num/denom));
 191 }
 192
 193
 194 template<class Point, class PointRef>
 195 inline Foam::scalar Foam::tetrahedron<Point, PointRef>::quality() const
 196 {
 197     return
 198         mag()
 199        /(
 200             8.0/(9.0*sqrt(3.0))
 201            *pow3(min(circumRadius(), GREAT))
 202           + ROOTVSMALL
 203         );
 204 }
 205
 206
 207 template<class Point, class PointRef>
 208 inline Point Foam::tetrahedron<Point, PointRef>::randomPoint
 209 (
 210     Random& rndGen
 211 ) const
 212 {
 213     // Adapted from
 214     // http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html
 215
 216     scalar s = rndGen.scalar01();
 217     scalar t = rndGen.scalar01();
 218     scalar u = rndGen.scalar01();
 219
 220     if (s + t > 1.0)
 221     {
 222         s = 1.0 - s;
 223         t = 1.0 - t;
 224     }
 225
 226     if (t + u > 1.0)
 227     {
 228         scalar tmp = u;
 229         u = 1.0 - s - t;
 230         t = 1.0 - tmp;
 231     }
 232     else if (s + t + u > 1.0)
 233     {
 234         scalar tmp = u;
 235         u = s + t + u - 1.0;
 236         s = 1.0 - t - tmp;
 237     }
 238
 239     return (1 - s - t - u)*a_ + s*b_ + t*c_ + u*d_;
 240 }
 241
 242
 243 template<class Point, class PointRef>
 244 inline Point Foam::tetrahedron<Point, PointRef>::randomPoint
 245 (
 246     cachedRandom& rndGen
 247 ) const
 248 {
 249     // Adapted from
 250     // http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html
 251
 252     scalar s = rndGen.sample01<scalar>();
 253     scalar t = rndGen.sample01<scalar>();
 254     scalar u = rndGen.sample01<scalar>();
 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>
 280 Foam::scalar Foam::tetrahedron<Point, PointRef>::barycentric
 281 (
 282     const point& pt,
 283     List<scalar>& bary
 284 ) const
 285 {
 286     // From:
 287     // http://en.wikipedia.org/wiki/Barycentric_coordinate_system_(mathematics)
 288
 289     vector e0(a_ - d_);
 290     vector e1(b_ - d_);
 291     vector e2(c_ - d_);
 292
 293     tensor t
 294     (
 295         e0.x(), e1.x(), e2.x(),
 296         e0.y(), e1.y(), e2.y(),
 297         e0.z(), e1.z(), e2.z()
 298     );
 299
 300     scalar detT = det(t);
 301
 302     if (Foam::mag(detT) < SMALL)
 303     {
 304         // Degenerate tetrahedron, returning 1/4 barycentric coordinates.
 305
 306         bary = List<scalar>(4, 0.25);
 307
 308         return detT;
 309     }
 310
 311     vector res = inv(t, detT) & (pt - d_);
 312
 313     bary.setSize(4);
 314
 315     bary[0] = res.x();
 316     bary[1] = res.y();
 317     bary[2] = res.z();
 318     bary[3] = (1.0 - res.x() - res.y() - res.z());
 319
 320     return detT;
 321 }
 322
 323
 324 template<class Point, class PointRef>
 325 inline Foam::pointHit Foam::tetrahedron<Point, PointRef>::nearestPoint
 326 (
 327     const point& p
 328 ) const
 329 {
 330     // Adapted from:
 331     // Real-time collision detection, Christer Ericson, 2005, p142-144
 332
 333     // Assuming initially that the point is inside all of the
 334     // halfspaces, so closest to itself.
 335
 336     point closestPt = p;
 337
 338     scalar minOutsideDistance = VGREAT;
 339
 340     bool inside = true;
 341
 342     if (((p - b_) & Sa()) >= 0)
 343     {
 344         // p is outside halfspace plane of tri
 345         pointHit info = triangle<Point, PointRef>(b_, c_, d_).nearestPoint(p);
 346
 347         inside = false;
 348
 349         if (info.distance() < minOutsideDistance)
 350         {
 351             closestPt = info.rawPoint();
 352
 353             minOutsideDistance = info.distance();
 354         }
 355     }
 356
 357     if (((p - a_) & Sb()) >= 0)
 358     {
 359         // p is outside halfspace plane of tri
 360         pointHit info = triangle<Point, PointRef>(a_, d_, c_).nearestPoint(p);
 361
 362         inside = false;
 363
 364         if (info.distance() < minOutsideDistance)
 365         {
 366             closestPt = info.rawPoint();
 367
 368             minOutsideDistance = info.distance();
 369         }
 370     }
 371
 372     if (((p - a_) & Sc()) >= 0)
 373     {
 374         // p is outside halfspace plane of tri
 375         pointHit info = triangle<Point, PointRef>(a_, b_, d_).nearestPoint(p);
 376
 377         inside = false;
 378
 379         if (info.distance() < minOutsideDistance)
 380         {
 381             closestPt = info.rawPoint();
 382
 383             minOutsideDistance = info.distance();
 384         }
 385     }
 386
 387     if (((p - a_) & Sd()) >= 0)
 388     {
 389         // p is outside halfspace plane of tri
 390         pointHit info = triangle<Point, PointRef>(a_, c_, b_).nearestPoint(p);
 391
 392         inside = false;
 393
 394         if (info.distance() < minOutsideDistance)
 395         {
 396             closestPt = info.rawPoint();
 397
 398             minOutsideDistance = info.distance();
 399         }
 400     }
 401
 402     // If the point is inside, then the distance to the closest point
 403     // is zero
 404     if (inside)
 405     {
 406         minOutsideDistance = 0;
 407     }
 408
 409     return pointHit
 410     (
 411         inside,
 412         closestPt,
 413         minOutsideDistance,
 414         !inside
 415     );
 416 }
 417
 418
 419 template<class Point, class PointRef>
 420 bool Foam::tetrahedron<Point, PointRef>::inside(const point& pt) const
 421 {
 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.
 425
 426     // The tet is defined: tet(Cc, tetBasePt, pA, pB) where the normal
 427     // vectors and base points for the half-space planes are:
 428     // area[0] = Sa();
 429     // area[1] = Sb();
 430     // area[2] = Sc();
 431     // area[3] = Sd();
 432     // planeBase[0] = tetBasePt = b_
 433     // planeBase[1] = ptA       = c_
 434     // planeBase[2] = tetBasePt = b_
 435     // planeBase[3] = tetBasePt = b_
 436
 437     vector n = vector::zero;
 438
 439     {
 440         // 0, a
 441         const point& basePt = b_;
 442
 443         n = Sa();
 444         n /= (Foam::mag(n) + VSMALL);
 445
 446         if (((pt - basePt) & n) > SMALL)
 447         {
 448             return false;
 449         }
 450     }
 451
 452     {
 453         // 1, b
 454         const point& basePt = c_;
 455
 456         n = Sb();
 457         n /= (Foam::mag(n) + VSMALL);
 458
 459         if (((pt - basePt) & n) > SMALL)
 460         {
 461             return false;
 462         }
 463     }
 464
 465     {
 466         // 2, c
 467         const point& basePt = b_;
 468
 469         n = Sc();
 470         n /= (Foam::mag(n) + VSMALL);
 471
 472         if (((pt - basePt) & n) > SMALL)
 473         {
 474             return false;
 475         }
 476     }
 477
 478     {
 479         // 3, d
 480         const point& basePt = b_;
 481
 482         n = Sd();
 483         n /= (Foam::mag(n) + VSMALL);
 484
 485         if (((pt - basePt) & n) > SMALL)
 486         {
 487             return false;
 488         }
 489     }
 490
 491     return true;
 492 }
 493
 494
 495 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 496
 497 template<class Point, class PointRef>
 498 inline Foam::Istream& Foam::operator>>
 499 (
 500     Istream& is,
 501     tetrahedron<Point, PointRef>& t
 502 )
 503 {
 504     is.readBegin("tetrahedron");
 505     is  >> t.a_ >> t.b_ >> t.c_ >> t.d_;
 506     is.readEnd("tetrahedron");
 507
 508     is.check("Istream& operator>>(Istream&, tetrahedron&)");
 509
 510     return is;
 511 }
 512
 513
 514 template<class Point, class PointRef>
 515 inline Foam::Ostream& Foam::operator<<
 516 (
 517     Ostream& os,
 518     const tetrahedron<Point, PointRef>& t
 519 )
 520 {
 521     os  << nl
 522         << token::BEGIN_LIST
 523         << t.a_ << token::SPACE
 524         << t.b_ << token::SPACE
 525         << t.c_ << token::SPACE
 526         << t.d_
 527         << token::END_LIST;
 528
 529     return os;
 530 }
 531
 532
 533 // ************************************************************************* //