src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 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 "IOstreams.H"
  27 #include "pointHit.H"
  28 #include "mathematicalConstants.H"
  29
  30 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
  31
  32 template<class Point, class PointRef>
  33 inline Foam::triangle<Point, PointRef>::triangle
  34 (
  35     const Point& a,
  36     const Point& b,
  37     const Point& c
  38 )
  39 :
  40     a_(a),
  41     b_(b),
  42     c_(c)
  43 {}
  44
  45
  46 template<class Point, class PointRef>
  47 inline Foam::triangle<Point, PointRef>::triangle
  48 (
  49     const UList<Point>& points,
  50     const FixedList<label, 3>& indices
  51 )
  52 :
  53     a_(points[indices[0]]),
  54     b_(points[indices[1]]),
  55     c_(points[indices[2]])
  56 {}
  57
  58
  59
  60 template<class Point, class PointRef>
  61 inline Foam::triangle<Point, PointRef>::triangle(Istream& is)
  62 {
  63     is  >> *this;
  64 }
  65
  66
  67 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  68
  69 template<class Point, class PointRef>
  70 inline const Point& Foam::triangle<Point, PointRef>::a() const
  71 {
  72     return a_;
  73 }
  74
  75 template<class Point, class PointRef>
  76 inline const Point& Foam::triangle<Point, PointRef>::b() const
  77 {
  78     return b_;
  79 }
  80
  81 template<class Point, class PointRef>
  82 inline const Point& Foam::triangle<Point, PointRef>::c() const
  83 {
  84     return c_;
  85 }
  86
  87
  88 template<class Point, class PointRef>
  89 inline Point Foam::triangle<Point, PointRef>::centre() const
  90 {
  91     return (1.0/3.0)*(a_ + b_ + c_);
  92 }
  93
  94
  95 template<class Point, class PointRef>
  96 inline Foam::scalar Foam::triangle<Point, PointRef>::mag() const
  97 {
  98     return Foam::mag(normal());
  99 }
 100
 101
 102 template<class Point, class PointRef>
 103 inline Foam::vector Foam::triangle<Point, PointRef>::normal() const
 104 {
 105     return 0.5*((b_ - a_)^(c_ - a_));
 106 }
 107
 108
 109 template<class Point, class PointRef>
 110 inline Point Foam::triangle<Point, PointRef>::circumCentre() const
 111 {
 112     scalar d1 =  (c_ - a_) & (b_ - a_);
 113     scalar d2 = -(c_ - b_) & (b_ - a_);
 114     scalar d3 =  (c_ - a_) & (c_ - b_);
 115
 116     scalar c1 = d2*d3;
 117     scalar c2 = d3*d1;
 118     scalar c3 = d1*d2;
 119
 120     scalar c = c1 + c2 + c3;
 121
 122     if (Foam::mag(c) < ROOTVSMALL)
 123     {
 124         // Degenerate triangle, returning centre instead of circumCentre.
 125
 126         return centre();
 127     }
 128
 129     return
 130     (
 131         ((c2 + c3)*a_ + (c3 + c1)*b_ + (c1 + c2)*c_)/(2*c)
 132     );
 133 }
 134
 135
 136 template<class Point, class PointRef>
 137 inline Foam::scalar Foam::triangle<Point, PointRef>::circumRadius() const
 138 {
 139     scalar d1 =  (c_ - a_) & (b_ - a_);
 140     scalar d2 = -(c_ - b_) & (b_ - a_);
 141     scalar d3 =  (c_ - a_) & (c_ - b_);
 142
 143     scalar denom = d2*d3 + d3*d1 + d1*d2;
 144
 145     if (Foam::mag(denom) < VSMALL)
 146     {
 147         // Degenerate triangle, returning GREAT for circumRadius.
 148
 149         return GREAT;
 150     }
 151     else
 152     {
 153         scalar a = (d1 + d2)*(d2 + d3)*(d3 + d1) / denom;
 154
 155         return 0.5*Foam::sqrt(min(GREAT, max(0, a)));
 156     }
 157 }
 158
 159
 160 template<class Point, class PointRef>
 161 inline Foam::scalar Foam::triangle<Point, PointRef>::quality() const
 162 {
 163     return mag()/(Foam::sqr(circumRadius())*3.0*sqrt(3.0)/4.0 + VSMALL);
 164 }
 165
 166
 167 template<class Point, class PointRef>
 168 inline Foam::scalar Foam::triangle<Point, PointRef>::sweptVol
 169 (
 170     const triangle& t
 171 ) const
 172 {
 173     return (1.0/12.0)*
 174     (
 175         ((t.a_ - a_) & ((b_ - a_)^(c_ - a_)))
 176       + ((t.b_ - b_) & ((c_ - b_)^(t.a_ - b_)))
 177       + ((c_ - t.c_) & ((t.b_ - t.c_)^(t.a_ - t.c_)))
 178
 179       + ((t.a_ - a_) & ((b_ - a_)^(c_ - a_)))
 180       + ((b_ - t.b_) & ((t.a_ - t.b_)^(t.c_ - t.b_)))
 181       + ((c_ - t.c_) & ((b_ - t.c_)^(t.a_ - t.c_)))
 182     );
 183 }
 184
 185
 186 template<class Point, class PointRef>
 187 inline Foam::tensor Foam::triangle<Point, PointRef>::inertia
 188 (
 189     PointRef refPt,
 190     scalar density
 191 ) const
 192 {
 193     Point aRel = a_ - refPt;
 194     Point bRel = b_ - refPt;
 195     Point cRel = c_ - refPt;
 196
 197     tensor V
 198     (
 199         aRel.x(), aRel.y(), aRel.z(),
 200         bRel.x(), bRel.y(), bRel.z(),
 201         cRel.x(), cRel.y(), cRel.z()
 202     );
 203
 204     scalar a = Foam::mag((b_ - a_)^(c_ - a_));
 205
 206     tensor S = 1/24.0*(tensor::one + I);
 207
 208     return
 209     (
 210         a*I/24.0*
 211         (
 212             (aRel & aRel)
 213           + (bRel & bRel)
 214           + (cRel & cRel)
 215           + ((aRel + bRel + cRel) & (aRel + bRel + cRel))
 216         )
 217       - a*(V.T() & S & V)
 218     )
 219    *density;
 220 }
 221
 222
 223 template<class Point, class PointRef>
 224 inline Point Foam::triangle<Point, PointRef>::randomPoint(Random& rndGen) const
 225 {
 226     // Generating Random Points in Triangles
 227     // by Greg Turk
 228     // from "Graphics Gems", Academic Press, 1990
 229     // http://tog.acm.org/GraphicsGems/gems/TriPoints.c
 230
 231     scalar s = rndGen.scalar01();
 232
 233     scalar t = sqrt(rndGen.scalar01());
 234
 235     return (1 - t)*a_ + (1 - s)*t*b_ + s*t*c_;
 236 }
 237
 238
 239 template<class Point, class PointRef>
 240 Foam::scalar Foam::triangle<Point, PointRef>::barycentric
 241 (
 242     const point& pt,
 243     List<scalar>& bary
 244 ) const
 245 {
 246     // From:
 247     // Real-time collision detection, Christer Ericson, 2005, p47-48
 248
 249     vector v0 = b_ - a_;
 250     vector v1 = c_ - a_;
 251     vector v2 = pt - a_;
 252
 253     scalar d00 = v0 & v0;
 254     scalar d01 = v0 & v1;
 255     scalar d11 = v1 & v1;
 256     scalar d20 = v2 & v0;
 257     scalar d21 = v2 & v1;
 258
 259     scalar denom = d00*d11 - d01*d01;
 260
 261     if (Foam::mag(denom) < SMALL)
 262     {
 263         // Degenerate triangle, returning 1/3 barycentric coordinates.
 264
 265         bary = List<scalar>(3, 1.0/3.0);
 266
 267         return denom;
 268     }
 269
 270     bary.setSize(3);
 271
 272     bary[1] = (d11*d20 - d01*d21)/denom;
 273     bary[2] = (d00*d21 - d01*d20)/denom;
 274     bary[0] = 1.0 - bary[1] - bary[2];
 275
 276     return denom;
 277 }
 278
 279
 280 template<class Point, class PointRef>
 281 inline Foam::pointHit Foam::triangle<Point, PointRef>::ray
 282 (
 283     const point& p,
 284     const vector& q,
 285     const intersection::algorithm alg,
 286     const intersection::direction dir
 287 ) const
 288 {
 289     // Express triangle in terms of baseVertex (point a_) and
 290     // two edge vectors
 291     vector E0 = b_ - a_;
 292     vector E1 = c_ - a_;
 293
 294     // Initialize intersection to miss.
 295     pointHit inter(p);
 296
 297     vector n(0.5*(E0 ^ E1));
 298     scalar area = Foam::mag(n);
 299
 300     if (area < VSMALL)
 301     {
 302         // Ineligible miss.
 303         inter.setMiss(false);
 304
 305         // The miss point is the nearest point on the triangle. Take any one.
 306         inter.setPoint(a_);
 307
 308         // The distance to the miss is the distance between the
 309         // original point and plane of intersection. No normal so use
 310         // distance from p to a_
 311         inter.setDistance(Foam::mag(a_ - p));
 312
 313         return inter;
 314     }
 315
 316     vector q1 = q/Foam::mag(q);
 317
 318     if (dir == intersection::CONTACT_SPHERE)
 319     {
 320         n /= area;
 321
 322         return ray(p, q1 - n, alg, intersection::VECTOR);
 323     }
 324
 325     // Intersection point with triangle plane
 326     point pInter;
 327     // Is intersection point inside triangle
 328     bool hit;
 329     {
 330         // Reuse the fast ray intersection routine below in FULL_RAY
 331         // mode since the original intersection routine has rounding problems.
 332         pointHit fastInter = intersection(p, q1, intersection::FULL_RAY);
 333         hit = fastInter.hit();
 334
 335         if (hit)
 336         {
 337             pInter = fastInter.rawPoint();
 338         }
 339         else
 340         {
 341             // Calculate intersection of ray with triangle plane
 342             vector v = a_ - p;
 343             pInter = p + (q1&v)*q1;
 344         }
 345     }
 346
 347     // Distance to intersection point
 348     scalar dist = q1 & (pInter - p);
 349
 350     const scalar planarPointTol =
 351         Foam::min
 352         (
 353             Foam::min
 354             (
 355                 Foam::mag(E0),
 356                 Foam::mag(E1)
 357             ),
 358             Foam::mag(c_ - b_)
 359         )*intersection::planarTol();
 360
 361     bool eligible =
 362         alg == intersection::FULL_RAY
 363      || (alg == intersection::HALF_RAY && dist > -planarPointTol)
 364      || (
 365             alg == intersection::VISIBLE
 366          && ((q1 & normal()) < -VSMALL)
 367         );
 368
 369     if (hit && eligible)
 370     {
 371         // Hit. Set distance to intersection.
 372         inter.setHit();
 373         inter.setPoint(pInter);
 374         inter.setDistance(dist);
 375     }
 376     else
 377     {
 378         // Miss or ineligible hit.
 379         inter.setMiss(eligible);
 380
 381         // The miss point is the nearest point on the triangle
 382         inter.setPoint(nearestPoint(p).rawPoint());
 383
 384         // The distance to the miss is the distance between the
 385         // original point and plane of intersection
 386         inter.setDistance(Foam::mag(pInter - p));
 387     }
 388
 389
 390     return inter;
 391 }
 392
 393
 394 // From "Fast, Minimum Storage Ray/Triangle Intersection"
 395 // Moeller/Trumbore.
 396 template<class Point, class PointRef>
 397 inline Foam::pointHit Foam::triangle<Point, PointRef>::intersection
 398 (
 399     const point& orig,
 400     const vector& dir,
 401     const intersection::algorithm alg,
 402     const scalar tol
 403 ) const
 404 {
 405     const vector edge1 = b_ - a_;
 406     const vector edge2 = c_ - a_;
 407
 408     // begin calculating determinant - also used to calculate U parameter
 409     const vector pVec = dir ^ edge2;
 410
 411     // if determinant is near zero, ray lies in plane of triangle
 412     const scalar det = edge1 & pVec;
 413
 414     // Initialise to miss
 415     pointHit intersection(false, vector::zero, GREAT, false);
 416
 417     if (alg == intersection::VISIBLE)
 418     {
 419         // Culling branch
 420         if (det < ROOTVSMALL)
 421         {
 422             // Ray on wrong side of triangle. Return miss
 423             return intersection;
 424         }
 425     }
 426     else if (alg == intersection::HALF_RAY || alg == intersection::FULL_RAY)
 427     {
 428         // Non-culling branch
 429         if (det > -ROOTVSMALL && det < ROOTVSMALL)
 430         {
 431             // Ray parallel to triangle. Return miss
 432             return intersection;
 433         }
 434     }
 435
 436     const scalar inv_det = 1.0 / det;
 437
 438     /* calculate distance from a_ to ray origin */
 439     const vector tVec = orig-a_;
 440
 441     /* calculate U parameter and test bounds */
 442     const scalar u = (tVec & pVec)*inv_det;
 443
 444     if (u < -tol || u > 1.0+tol)
 445     {
 446         // return miss
 447         return intersection;
 448     }
 449
 450     /* prepare to test V parameter */
 451     const vector qVec = tVec ^ edge1;
 452
 453     /* calculate V parameter and test bounds */
 454     const scalar v = (dir & qVec) * inv_det;
 455
 456     if (v < -tol || u + v > 1.0+tol)
 457     {
 458         // return miss
 459         return intersection;
 460     }
 461
 462     /* calculate t, scale parameters, ray intersects triangle */
 463     const scalar t = (edge2 & qVec) * inv_det;
 464
 465     if (alg == intersection::HALF_RAY && t < -tol)
 466     {
 467         // Wrong side of orig. Return miss
 468         return intersection;
 469     }
 470
 471     intersection.setHit();
 472     intersection.setPoint(a_ + u*edge1 + v*edge2);
 473     intersection.setDistance(t);
 474
 475     return intersection;
 476 }
 477
 478
 479 template<class Point, class PointRef>
 480 Foam::pointHit Foam::triangle<Point, PointRef>::nearestPointClassify
 481 (
 482     const point& p,
 483     label& nearType,
 484     label& nearLabel
 485 ) const
 486 {
 487     // Adapted from:
 488     // Real-time collision detection, Christer Ericson, 2005, p136-142
 489
 490     // Check if P in vertex region outside A
 491     vector ab = b_ - a_;
 492     vector ac = c_ - a_;
 493     vector ap = p - a_;
 494
 495     scalar d1 = ab & ap;
 496     scalar d2 = ac & ap;
 497
 498     if (d1 <= 0.0 && d2 <= 0.0)
 499     {
 500         // barycentric coordinates (1, 0, 0)
 501
 502         nearType = POINT;
 503         nearLabel = 0;
 504         return pointHit(false, a_, Foam::mag(a_ - p), true);
 505     }
 506
 507     // Check if P in vertex region outside B
 508     vector bp = p - b_;
 509     scalar d3 = ab & bp;
 510     scalar d4 = ac & bp;
 511
 512     if (d3 >= 0.0 && d4 <= d3)
 513     {
 514         // barycentric coordinates (0, 1, 0)
 515
 516         nearType = POINT;
 517         nearLabel = 1;
 518         return pointHit(false, b_, Foam::mag(b_ - p), true);
 519     }
 520
 521     // Check if P in edge region of AB, if so return projection of P onto AB
 522     scalar vc = d1*d4 - d3*d2;
 523
 524     if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0)
 525     {
 526         if ((d1 - d3) < ROOTVSMALL)
 527         {
 528             // Degenerate triangle, for d1 = d3, a_ and b_ are likely coincident
 529             nearType = POINT;
 530             nearLabel = 0;
 531             return pointHit(false, a_, Foam::mag(a_ - p), true);
 532         }
 533
 534         // barycentric coordinates (1-v, v, 0)
 535         scalar v = d1/(d1 - d3);
 536
 537         point nearPt =  a_ + v*ab;
 538         nearType = EDGE;
 539         nearLabel = 0;
 540         return pointHit(false, nearPt, Foam::mag(nearPt - p), true);
 541     }
 542
 543     // Check if P in vertex region outside C
 544     vector cp = p - c_;
 545     scalar d5 = ab & cp;
 546     scalar d6 = ac & cp;
 547
 548     if (d6 >= 0.0 && d5 <= d6)
 549     {
 550         // barycentric coordinates (0, 0, 1)
 551
 552         nearType = POINT;
 553         nearLabel = 2;
 554         return pointHit(false, c_, Foam::mag(c_ - p), true);
 555     }
 556
 557     // Check if P in edge region of AC, if so return projection of P onto AC
 558     scalar vb = d5*d2 - d1*d6;
 559
 560     if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0)
 561     {
 562         if ((d2 - d6) < ROOTVSMALL)
 563         {
 564             // Degenerate triangle, for d2 = d6, a_ and c_ are likely coincident
 565             nearType = POINT;
 566             nearLabel = 0;
 567             return pointHit(false, a_, Foam::mag(a_ - p), true);
 568         }
 569
 570         // barycentric coordinates (1-w, 0, w)
 571         scalar w = d2/(d2 - d6);
 572
 573         point nearPt = a_ + w*ac;
 574         nearType = EDGE;
 575         nearLabel = 2;
 576         return pointHit(false, nearPt, Foam::mag(nearPt - p), true);
 577     }
 578
 579     // Check if P in edge region of BC, if so return projection of P onto BC
 580     scalar va = d3*d6 - d5*d4;
 581
 582     if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0)
 583     {
 584         if (((d4 - d3) + (d5 - d6)) < ROOTVSMALL)
 585         {
 586             // Degenerate triangle, for (d4 - d3) = (d6 - d5), b_ and c_ are
 587             // likely coincident
 588             nearType = POINT;
 589             nearLabel = 1;
 590             return pointHit(false, b_, Foam::mag(b_ - p), true);
 591         }
 592
 593         // barycentric coordinates (0, 1-w, w)
 594         scalar w = (d4 - d3)/((d4 - d3) + (d5 - d6));
 595
 596         point nearPt = b_ + w*(c_ - b_);
 597         nearType = EDGE;
 598         nearLabel = 1;
 599         return pointHit(false, nearPt, Foam::mag(nearPt - p), true);
 600     }
 601
 602     // P inside face region. Compute Q through its barycentric
 603     // coordinates (u, v, w)
 604
 605     if ((va + vb + vc) < ROOTVSMALL)
 606     {
 607         // Degenerate triangle, return the centre because no edge or points are
 608         // closest
 609         point nearPt = centre();
 610         nearType = NONE,
 611         nearLabel = -1;
 612         return pointHit(true, nearPt, Foam::mag(nearPt - p), false);
 613     }
 614
 615     scalar denom = 1.0/(va + vb + vc);
 616     scalar v = vb * denom;
 617     scalar w = vc * denom;
 618
 619     // = u*a + v*b + w*c, u = va*denom = 1.0 - v - w
 620
 621     point nearPt = a_ + ab*v + ac*w;
 622     nearType = NONE,
 623     nearLabel = -1;
 624     return pointHit(true, nearPt, Foam::mag(nearPt - p), false);
 625 }
 626
 627
 628 template<class Point, class PointRef>
 629 inline Foam::pointHit Foam::triangle<Point, PointRef>::nearestPoint
 630 (
 631     const point& p
 632 ) const
 633 {
 634     // Dummy labels
 635     label nearType = -1;
 636     label nearLabel = -1;
 637
 638     return nearestPointClassify(p, nearType, nearLabel);
 639 }
 640
 641
 642 template<class Point, class PointRef>
 643 inline bool Foam::triangle<Point, PointRef>::classify
 644 (
 645     const point& p,
 646     label& nearType,
 647     label& nearLabel
 648 ) const
 649 {
 650     return nearestPointClassify(p, nearType, nearLabel).hit();
 651 }
 652
 653
 654 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 655
 656 template<class Point, class PointRef>
 657 inline Foam::Istream& Foam::operator>>
 658 (
 659     Istream& is,
 660     triangle<Point, PointRef>& t
 661 )
 662 {
 663     is.readBegin("triangle");
 664     is  >> t.a_ >> t.b_ >> t.c_;
 665     is.readEnd("triangle");
 666
 667     is.check("Istream& operator>>(Istream&, triangle&)");
 668     return is;
 669 }
 670
 671
 672 template<class Point, class PointRef>
 673 inline Foam::Ostream& Foam::operator<<
 674 (
 675     Ostream& os,
 676     const triangle<Point, PointRef>& t
 677 )
 678 {
 679     os  << nl
 680         << token::BEGIN_LIST
 681         << t.a_ << token::SPACE
 682         << t.b_ << token::SPACE
 683         << t.c_
 684         << token::END_LIST;
 685
 686     return os;
 687 }
 688
 689
 690 // ************************************************************************* //