triangleI.H 21.8 KB
Newer Older
1 2 3 4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
OpenFOAM bot's avatar
OpenFOAM bot committed
5
    \\  /    A nd           | www.openfoam.com
OpenFOAM bot's avatar
OpenFOAM bot committed
6 7
     \\/     M anipulation  |
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8 9
    Copyright (C) 2011-2017 OpenFOAM Foundation
    Copyright (C) 2018 OpenCFD Ltd.
10 11 12 13
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

14 15 16 17
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
18 19 20 21 22 23 24

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
25
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
26 27 28 29 30

\*---------------------------------------------------------------------------*/

#include "IOstreams.H"
#include "pointHit.H"
31
#include "triPoints.H"
32
#include "mathematicalConstants.H"
33 34 35 36

// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

template<class Point, class PointRef>
37
inline Foam::triangle<Point, PointRef>::triangle
38 39 40 41 42 43 44 45 46 47 48 49
(
    const Point& a,
    const Point& b,
    const Point& c
)
:
    a_(a),
    b_(b),
    c_(c)
{}


50 51 52 53 54 55 56 57 58 59 60 61
template<class Point, class PointRef>
inline Foam::triangle<Point, PointRef>::triangle
(
    const FixedList<Point, 3>& tri
)
:
    a_(tri[0]),
    b_(tri[1]),
    c_(tri[2])
{}


62
template<class Point, class PointRef>
63 64 65 66 67 68 69 70 71 72
inline Foam::triangle<Point, PointRef>::triangle
(
    const UList<Point>& points,
    const FixedList<label, 3>& indices
)
:
    a_(points[indices[0]]),
    b_(points[indices[1]]),
    c_(points[indices[2]])
{}
73 74 75



76 77 78 79
template<class Point, class PointRef>
inline Foam::triangle<Point, PointRef>::triangle(Istream& is)
{
    is  >> *this;
80 81 82 83 84 85
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class Point, class PointRef>
86
inline const Point& Foam::triangle<Point, PointRef>::a() const
87 88 89 90 91
{
    return a_;
}

template<class Point, class PointRef>
92
inline const Point& Foam::triangle<Point, PointRef>::b() const
93 94 95 96 97
{
    return b_;
}

template<class Point, class PointRef>
98
inline const Point& Foam::triangle<Point, PointRef>::c() const
99 100 101 102 103 104
{
    return c_;
}


template<class Point, class PointRef>
105
inline Point Foam::triangle<Point, PointRef>::centre() const
106 107 108 109 110 111
{
    return (1.0/3.0)*(a_ + b_ + c_);
}


template<class Point, class PointRef>
112
inline Foam::vector Foam::triangle<Point, PointRef>::areaNormal() const
113
{
114
    return 0.5*((b_ - a_)^(c_ - a_));
115 116 117 118
}


template<class Point, class PointRef>
119
inline Foam::vector Foam::triangle<Point, PointRef>::unitNormal() const
120
{
121 122 123 124 125 126 127 128 129 130
    const vector n(areaNormal());
    const scalar s(Foam::mag(n));
    return s < ROOTVSMALL ? Zero : n/s;
}


template<class Point, class PointRef>
inline Foam::scalar Foam::triangle<Point, PointRef>::mag() const
{
    return ::Foam::mag(areaNormal());
131 132 133 134
}


template<class Point, class PointRef>
135
inline Point Foam::triangle<Point, PointRef>::circumCentre() const
136
{
137 138 139
    scalar d1 =  (c_ - a_) & (b_ - a_);
    scalar d2 = -(c_ - b_) & (b_ - a_);
    scalar d3 =  (c_ - a_) & (c_ - b_);
140 141 142 143 144 145 146

    scalar c1 = d2*d3;
    scalar c2 = d3*d1;
    scalar c3 = d1*d2;

    scalar c = c1 + c2 + c3;

147 148
    if (Foam::mag(c) < ROOTVSMALL)
    {
149
        // Degenerate triangle, returning centre instead of circumCentre.
150 151 152 153

        return centre();
    }

154
    return
155 156 157 158 159 160 161
    (
        ((c2 + c3)*a_ + (c3 + c1)*b_ + (c1 + c2)*c_)/(2*c)
    );
}


template<class Point, class PointRef>
162
inline Foam::scalar Foam::triangle<Point, PointRef>::circumRadius() const
163
{
164 165 166
    const scalar d1 =  (c_ - a_) & (b_ - a_);
    const scalar d2 = -(c_ - b_) & (b_ - a_);
    const scalar d3 =  (c_ - a_) & (c_ - b_);
167

168
    const scalar denom = d2*d3 + d3*d1 + d1*d2;
169 170 171

    if (Foam::mag(denom) < VSMALL)
    {
172
        // Degenerate triangle, returning GREAT for circumRadius.
173

174 175 176
        return GREAT;
    }

177 178
    const scalar a = (d1 + d2)*(d2 + d3)*(d3 + d1) / denom;
    return 0.5*Foam::sqrt(min(GREAT, max(0, a)));
179 180 181 182
}


template<class Point, class PointRef>
183
inline Foam::scalar Foam::triangle<Point, PointRef>::quality() const
184
{
185 186 187 188 189 190 191 192 193
    scalar c = circumRadius();

    if (c < ROOTVSMALL)
    {
        // zero circumRadius, something has gone wrong.
        return SMALL;
    }

    return mag()/(Foam::sqr(c)*3.0*sqrt(3.0)/4.0);
194 195 196 197
}


template<class Point, class PointRef>
198 199 200 201
inline Foam::scalar Foam::triangle<Point, PointRef>::sweptVol
(
    const triangle& t
) const
202
{
203
    return (1.0/12.0)*
204 205 206 207
    (
        ((t.a_ - a_) & ((b_ - a_)^(c_ - a_)))
      + ((t.b_ - b_) & ((c_ - b_)^(t.a_ - b_)))
      + ((c_ - t.c_) & ((t.b_ - t.c_)^(t.a_ - t.c_)))
208 209 210 211

      + ((t.a_ - a_) & ((b_ - a_)^(c_ - a_)))
      + ((b_ - t.b_) & ((t.a_ - t.b_)^(t.c_ - t.b_)))
      + ((c_ - t.c_) & ((b_ - t.c_)^(t.a_ - t.c_)))
212 213 214 215
    );
}


216
template<class Point, class PointRef>
217
inline Foam::tensor Foam::triangle<Point, PointRef>::inertia
218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251
(
    PointRef refPt,
    scalar density
) const
{
    Point aRel = a_ - refPt;
    Point bRel = b_ - refPt;
    Point cRel = c_ - refPt;

    tensor V
    (
        aRel.x(), aRel.y(), aRel.z(),
        bRel.x(), bRel.y(), bRel.z(),
        cRel.x(), cRel.y(), cRel.z()
    );

    scalar a = Foam::mag((b_ - a_)^(c_ - a_));

    tensor S = 1/24.0*(tensor::one + I);

    return
    (
        a*I/24.0*
        (
            (aRel & aRel)
          + (bRel & bRel)
          + (cRel & cRel)
          + ((aRel + bRel + cRel) & (aRel + bRel + cRel))
        )
      - a*(V.T() & S & V)
    )
   *density;
}

252 253

template<class Point, class PointRef>
254
inline Point Foam::triangle<Point, PointRef>::randomPoint(Random& rndGen) const
255
{
256 257
    return barycentricToPoint(barycentric2D01(rndGen));
}
258

259

260 261 262 263 264 265 266 267
template<class Point, class PointRef>
inline Point Foam::triangle<Point, PointRef>::barycentricToPoint
(
    const barycentric2D& bary
) const
{
    return bary[0]*a_ + bary[1]*b_ + bary[2]*c_;
}
268

269 270 271 272 273 274 275 276 277 278

template<class Point, class PointRef>
inline Foam::barycentric2D Foam::triangle<Point, PointRef>::pointToBarycentric
(
    const point& pt
) const
{
    barycentric2D bary;
    pointToBarycentric(pt, bary);
    return bary;
279 280 281
}


282
template<class Point, class PointRef>
283
inline Foam::scalar Foam::triangle<Point, PointRef>::pointToBarycentric
284 285
(
    const point& pt,
286
    barycentric2D& bary
287 288
) const
{
Henry Weller's avatar
Henry Weller committed
289
    // Reference:
290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305
    // Real-time collision detection, Christer Ericson, 2005, p47-48

    vector v0 = b_ - a_;
    vector v1 = c_ - a_;
    vector v2 = pt - a_;

    scalar d00 = v0 & v0;
    scalar d01 = v0 & v1;
    scalar d11 = v1 & v1;
    scalar d20 = v2 & v0;
    scalar d21 = v2 & v1;

    scalar denom = d00*d11 - d01*d01;

    if (Foam::mag(denom) < SMALL)
    {
306
        // Degenerate triangle, returning 1/3 barycentric coordinates.
307

308
        bary = barycentric2D(1.0/3.0, 1.0/3.0, 1.0/3.0);
309 310 311 312

        return denom;
    }

313 314 315
    bary[1] = (d11*d20 - d01*d21)/denom;
    bary[2] = (d00*d21 - d01*d20)/denom;
    bary[0] = 1.0 - bary[1] - bary[2];
316 317 318 319 320

    return denom;
}


321
template<class Point, class PointRef>
322
inline Foam::pointHit Foam::triangle<Point, PointRef>::ray
323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360
(
    const point& p,
    const vector& q,
    const intersection::algorithm alg,
    const intersection::direction dir
) const
{
    // Express triangle in terms of baseVertex (point a_) and
    // two edge vectors
    vector E0 = b_ - a_;
    vector E1 = c_ - a_;

    // Initialize intersection to miss.
    pointHit inter(p);

    vector n(0.5*(E0 ^ E1));
    scalar area = Foam::mag(n);

    if (area < VSMALL)
    {
        // Ineligible miss.
        inter.setMiss(false);

        // The miss point is the nearest point on the triangle. Take any one.
        inter.setPoint(a_);

        // The distance to the miss is the distance between the
        // original point and plane of intersection. No normal so use
        // distance from p to a_
        inter.setDistance(Foam::mag(a_ - p));

        return inter;
    }

    vector q1 = q/Foam::mag(q);

    if (dir == intersection::CONTACT_SPHERE)
    {
mattijs's avatar
mattijs committed
361 362
        n /= area;

363 364 365
        return ray(p, q1 - n, alg, intersection::VECTOR);
    }

366
    // Intersection point with triangle plane
367
    point pInter;
368
    // Is intersection point inside triangle
mattijs's avatar
mattijs committed
369 370 371 372 373 374
    bool hit;
    {
        // Reuse the fast ray intersection routine below in FULL_RAY
        // mode since the original intersection routine has rounding problems.
        pointHit fastInter = intersection(p, q1, intersection::FULL_RAY);
        hit = fastInter.hit();
375 376 377 378 379 380 381 382 383 384 385

        if (hit)
        {
            pInter = fastInter.rawPoint();
        }
        else
        {
            // Calculate intersection of ray with triangle plane
            vector v = a_ - p;
            pInter = p + (q1&v)*q1;
        }
mattijs's avatar
mattijs committed
386
    }
387

388
    // Distance to intersection point
389 390 391 392 393 394 395 396 397 398 399 400 401 402 403
    scalar dist = q1 & (pInter - p);

    const scalar planarPointTol =
        Foam::min
        (
            Foam::min
            (
                Foam::mag(E0),
                Foam::mag(E1)
            ),
            Foam::mag(c_ - b_)
        )*intersection::planarTol();

    bool eligible =
        alg == intersection::FULL_RAY
404
     || (alg == intersection::HALF_RAY && dist > -planarPointTol)
405 406
     || (
            alg == intersection::VISIBLE
407
         && ((q1 & areaNormal()) < -VSMALL)
408 409 410 411 412 413 414 415 416 417 418 419 420 421 422
        );

    if (hit && eligible)
    {
        // Hit. Set distance to intersection.
        inter.setHit();
        inter.setPoint(pInter);
        inter.setDistance(dist);
    }
    else
    {
        // Miss or ineligible hit.
        inter.setMiss(eligible);

        // The miss point is the nearest point on the triangle
423
        inter.setPoint(nearestPoint(p).rawPoint());
424 425 426 427 428 429

        // The distance to the miss is the distance between the
        // original point and plane of intersection
        inter.setDistance(Foam::mag(pInter - p));
    }

mattijs's avatar
mattijs committed
430

431 432 433 434
    return inter;
}


mattijs's avatar
mattijs committed
435 436 437
// From "Fast, Minimum Storage Ray/Triangle Intersection"
// Moeller/Trumbore.
template<class Point, class PointRef>
438
inline Foam::pointHit Foam::triangle<Point, PointRef>::intersection
mattijs's avatar
mattijs committed
439 440
(
    const point& orig,
441
    const vector& dir,
mattijs's avatar
mattijs committed
442 443
    const intersection::algorithm alg,
    const scalar tol
mattijs's avatar
mattijs committed
444 445 446 447 448 449 450 451 452 453 454 455
) const
{
    const vector edge1 = b_ - a_;
    const vector edge2 = c_ - a_;

    // begin calculating determinant - also used to calculate U parameter
    const vector pVec = dir ^ edge2;

    // if determinant is near zero, ray lies in plane of triangle
    const scalar det = edge1 & pVec;

    // Initialise to miss
Henry Weller's avatar
Henry Weller committed
456
    pointHit intersection(false, Zero, GREAT, false);
mattijs's avatar
mattijs committed
457 458 459 460

    if (alg == intersection::VISIBLE)
    {
        // Culling branch
mattijs's avatar
mattijs committed
461
        if (det < ROOTVSMALL)
mattijs's avatar
mattijs committed
462
        {
mattijs's avatar
mattijs committed
463
            // Ray on wrong side of triangle. Return miss
mattijs's avatar
mattijs committed
464 465 466 467 468 469
            return intersection;
        }
    }
    else if (alg == intersection::HALF_RAY || alg == intersection::FULL_RAY)
    {
        // Non-culling branch
mattijs's avatar
mattijs committed
470
        if (det > -ROOTVSMALL && det < ROOTVSMALL)
mattijs's avatar
mattijs committed
471
        {
mattijs's avatar
mattijs committed
472
            // Ray parallel to triangle. Return miss
mattijs's avatar
mattijs committed
473 474
            return intersection;
        }
mattijs's avatar
mattijs committed
475
    }
mattijs's avatar
mattijs committed
476

mattijs's avatar
mattijs committed
477
    const scalar inv_det = 1.0 / det;
mattijs's avatar
mattijs committed
478

mattijs's avatar
mattijs committed
479 480
    /* calculate distance from a_ to ray origin */
    const vector tVec = orig-a_;
mattijs's avatar
mattijs committed
481

mattijs's avatar
mattijs committed
482 483
    /* calculate U parameter and test bounds */
    const scalar u = (tVec & pVec)*inv_det;
mattijs's avatar
mattijs committed
484

mattijs's avatar
mattijs committed
485 486 487 488 489
    if (u < -tol || u > 1.0+tol)
    {
        // return miss
        return intersection;
    }
mattijs's avatar
mattijs committed
490

mattijs's avatar
mattijs committed
491 492 493 494 495 496 497 498 499 500
    /* prepare to test V parameter */
    const vector qVec = tVec ^ edge1;

    /* calculate V parameter and test bounds */
    const scalar v = (dir & qVec) * inv_det;

    if (v < -tol || u + v > 1.0+tol)
    {
        // return miss
        return intersection;
mattijs's avatar
mattijs committed
501
    }
mattijs's avatar
mattijs committed
502 503 504 505 506

    /* calculate t, scale parameters, ray intersects triangle */
    const scalar t = (edge2 & qVec) * inv_det;

    if (alg == intersection::HALF_RAY && t < -tol)
mattijs's avatar
mattijs committed
507
    {
mattijs's avatar
mattijs committed
508 509
        // Wrong side of orig. Return miss
        return intersection;
mattijs's avatar
mattijs committed
510 511
    }

mattijs's avatar
mattijs committed
512 513 514 515
    intersection.setHit();
    intersection.setPoint(a_ + u*edge1 + v*edge2);
    intersection.setDistance(t);

mattijs's avatar
mattijs committed
516 517 518
    return intersection;
}

519 520

template<class Point, class PointRef>
521
Foam::pointHit Foam::triangle<Point, PointRef>::nearestPointClassify
522 523 524 525 526 527
(
    const point& p,
    label& nearType,
    label& nearLabel
) const
{
528
    // Adapted from:
529
    // Real-time collision detection, Christer Ericson, 2005, p136-142
530

531 532 533 534
    // Check if P in vertex region outside A
    vector ab = b_ - a_;
    vector ac = c_ - a_;
    vector ap = p - a_;
535

536 537
    scalar d1 = ab & ap;
    scalar d2 = ac & ap;
538

539 540 541
    if (d1 <= 0.0 && d2 <= 0.0)
    {
        // barycentric coordinates (1, 0, 0)
542

543 544 545 546
        nearType = POINT;
        nearLabel = 0;
        return pointHit(false, a_, Foam::mag(a_ - p), true);
    }
547

548 549 550 551
    // Check if P in vertex region outside B
    vector bp = p - b_;
    scalar d3 = ab & bp;
    scalar d4 = ac & bp;
552

553 554 555
    if (d3 >= 0.0 && d4 <= d3)
    {
        // barycentric coordinates (0, 1, 0)
556

557 558 559 560
        nearType = POINT;
        nearLabel = 1;
        return pointHit(false, b_, Foam::mag(b_ - p), true);
    }
561

562 563
    // Check if P in edge region of AB, if so return projection of P onto AB
    scalar vc = d1*d4 - d3*d2;
564

565
    if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0)
566
    {
567 568
        if ((d1 - d3) < ROOTVSMALL)
        {
569
            // Degenerate triangle, for d1 = d3, a_ and b_ are likely coincident
570 571 572 573 574
            nearType = POINT;
            nearLabel = 0;
            return pointHit(false, a_, Foam::mag(a_ - p), true);
        }

575 576
        // barycentric coordinates (1-v, v, 0)
        scalar v = d1/(d1 - d3);
577

578 579 580 581
        point nearPt =  a_ + v*ab;
        nearType = EDGE;
        nearLabel = 0;
        return pointHit(false, nearPt, Foam::mag(nearPt - p), true);
582 583
    }

584 585 586 587
    // Check if P in vertex region outside C
    vector cp = p - c_;
    scalar d5 = ab & cp;
    scalar d6 = ac & cp;
588

589 590 591
    if (d6 >= 0.0 && d5 <= d6)
    {
        // barycentric coordinates (0, 0, 1)
592

593 594 595 596
        nearType = POINT;
        nearLabel = 2;
        return pointHit(false, c_, Foam::mag(c_ - p), true);
    }
597

598 599
    // Check if P in edge region of AC, if so return projection of P onto AC
    scalar vb = d5*d2 - d1*d6;
mattijs's avatar
mattijs committed
600

601
    if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0)
mattijs's avatar
mattijs committed
602
    {
603 604
        if ((d2 - d6) < ROOTVSMALL)
        {
605
            // Degenerate triangle, for d2 = d6, a_ and c_ are likely coincident
606 607 608 609 610
            nearType = POINT;
            nearLabel = 0;
            return pointHit(false, a_, Foam::mag(a_ - p), true);
        }

611 612 613 614 615 616 617
        // barycentric coordinates (1-w, 0, w)
        scalar w = d2/(d2 - d6);

        point nearPt = a_ + w*ac;
        nearType = EDGE;
        nearLabel = 2;
        return pointHit(false, nearPt, Foam::mag(nearPt - p), true);
mattijs's avatar
mattijs committed
618 619
    }

620 621
    // Check if P in edge region of BC, if so return projection of P onto BC
    scalar va = d3*d6 - d5*d4;
622

623
    if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0)
624
    {
625 626
        if (((d4 - d3) + (d5 - d6)) < ROOTVSMALL)
        {
627 628
            // Degenerate triangle, for (d4 - d3) = (d6 - d5), b_ and c_ are
            // likely coincident
629 630 631 632 633
            nearType = POINT;
            nearLabel = 1;
            return pointHit(false, b_, Foam::mag(b_ - p), true);
        }

634 635
        // barycentric coordinates (0, 1-w, w)
        scalar w = (d4 - d3)/((d4 - d3) + (d5 - d6));
636

637 638 639 640
        point nearPt = b_ + w*(c_ - b_);
        nearType = EDGE;
        nearLabel = 1;
        return pointHit(false, nearPt, Foam::mag(nearPt - p), true);
641 642
    }

643 644
    // P inside face region. Compute Q through its barycentric
    // coordinates (u, v, w)
645 646 647

    if ((va + vb + vc) < ROOTVSMALL)
    {
648 649
        // Degenerate triangle, return the centre because no edge or points are
        // closest
650 651 652 653 654 655
        point nearPt = centre();
        nearType = NONE,
        nearLabel = -1;
        return pointHit(true, nearPt, Foam::mag(nearPt - p), false);
    }

656 657 658
    scalar denom = 1.0/(va + vb + vc);
    scalar v = vb * denom;
    scalar w = vc * denom;
659

660
    // = u*a + v*b + w*c, u = va*denom = 1.0 - v - w
661

662 663 664 665
    point nearPt = a_ + ab*v + ac*w;
    nearType = NONE,
    nearLabel = -1;
    return pointHit(true, nearPt, Foam::mag(nearPt - p), false);
666 667 668
}


669
template<class Point, class PointRef>
670
inline Foam::pointHit Foam::triangle<Point, PointRef>::nearestPoint
671 672 673 674 675 676 677 678 679 680
(
    const point& p
) const
{
    // Dummy labels
    label nearType = -1;
    label nearLabel = -1;

    return nearestPointClassify(p, nearType, nearLabel);
}
681 682


683
template<class Point, class PointRef>
684
inline bool Foam::triangle<Point, PointRef>::classify
685 686 687 688 689 690 691 692 693
(
    const point& p,
    label& nearType,
    label& nearLabel
) const
{
    return nearestPointClassify(p, nearType, nearLabel).hit();
}

694

mattijs's avatar
mattijs committed
695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750
template<class Point, class PointRef>
inline Foam::pointHit Foam::triangle<Point, PointRef>::nearestPoint
(
    const linePointRef& ln,
    pointHit& lnInfo
) const
{
    vector q = ln.vec();
    pointHit triInfo
    (
        triangle<Point, PointRef>::intersection
        (
            ln.start(),
            q,
            intersection::FULL_RAY
        )
    );

    if (triInfo.hit())
    {
        // Line hits triangle. Find point on line.
        if (triInfo.distance() > 1)
        {
            // Hit beyond endpoint
            lnInfo.setMiss(true);
            lnInfo.setPoint(ln.end());
            scalar dist = Foam::mag(triInfo.hitPoint()-lnInfo.missPoint());
            lnInfo.setDistance(dist);
            triInfo.setMiss(true);
            triInfo.setDistance(dist);
        }
        else if (triInfo.distance() < 0)
        {
            // Hit beyond startpoint
            lnInfo.setMiss(true);
            lnInfo.setPoint(ln.start());
            scalar dist = Foam::mag(triInfo.hitPoint()-lnInfo.missPoint());
            lnInfo.setDistance(dist);
            triInfo.setMiss(true);
            triInfo.setDistance(dist);
        }
        else
        {
            // Hit on line
            lnInfo.setHit();
            lnInfo.setPoint(triInfo.hitPoint());
            lnInfo.setDistance(0.0);
            triInfo.setDistance(0.0);
        }
    }
    else
    {
        // Line skips triangle. See which triangle edge it gets closest to

        point nearestEdgePoint;
        point nearestLinePoint;
mattijs's avatar
mattijs committed
751
        //label minEdgeIndex = 0;
mattijs's avatar
mattijs committed
752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772
        scalar minDist = ln.nearestDist
        (
            linePointRef(a_, b_),
            nearestLinePoint,
            nearestEdgePoint
        );

        {
            point linePoint;
            point triEdgePoint;
            scalar dist = ln.nearestDist
            (
                linePointRef(b_, c_),
                linePoint,
                triEdgePoint
            );
            if (dist < minDist)
            {
                minDist = dist;
                nearestEdgePoint = triEdgePoint;
                nearestLinePoint = linePoint;
mattijs's avatar
mattijs committed
773
                //minEdgeIndex = 1;
mattijs's avatar
mattijs committed
774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790
            }
        }

        {
            point linePoint;
            point triEdgePoint;
            scalar dist = ln.nearestDist
            (
                linePointRef(c_, a_),
                linePoint,
                triEdgePoint
            );
            if (dist < minDist)
            {
                minDist = dist;
                nearestEdgePoint = triEdgePoint;
                nearestLinePoint = linePoint;
mattijs's avatar
mattijs committed
791
                //minEdgeIndex = 2;
mattijs's avatar
mattijs committed
792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820
            }
        }

        lnInfo.setDistance(minDist);
        triInfo.setDistance(minDist);
        triInfo.setMiss(false);
        triInfo.setPoint(nearestEdgePoint);

        // Convert point on line to pointHit
        if (Foam::mag(nearestLinePoint-ln.start()) < SMALL)
        {
            lnInfo.setMiss(true);
            lnInfo.setPoint(ln.start());
        }
        else if (Foam::mag(nearestLinePoint-ln.end()) < SMALL)
        {
            lnInfo.setMiss(true);
            lnInfo.setPoint(ln.end());
        }
        else
        {
            lnInfo.setHit();
            lnInfo.setPoint(nearestLinePoint);
        }
    }
    return triInfo;
}


821 822 823 824 825 826 827 828 829 830 831 832 833
template<class Point, class PointRef>
inline int Foam::triangle<Point, PointRef>::sign
(
    const point& p,
    const scalar tol
) const
{
    const scalar dist = ((p - a_) & unitNormal());

    return ((dist < -tol) ? -1 : (dist > tol) ? +1 : 0);
}


834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893
template<class Point, class PointRef>
inline void Foam::triangle<Point, PointRef>::dummyOp::operator()
(
    const triPoints&
)
{}


template<class Point, class PointRef>
inline Foam::triangle<Point, PointRef>::sumAreaOp::sumAreaOp()
:
    area_(0.0)
{}


template<class Point, class PointRef>
inline void Foam::triangle<Point, PointRef>::sumAreaOp::operator()
(
    const triPoints& tri
)
{
    area_ += triangle<Point, const Point&>(tri).mag();
}


template<class Point, class PointRef>
inline Foam::triangle<Point, PointRef>::storeOp::storeOp
(
    triIntersectionList& tris,
    label& nTris
)
:
    tris_(tris),
    nTris_(nTris)
{}


template<class Point, class PointRef>
inline void Foam::triangle<Point, PointRef>::storeOp::operator()
(
    const triPoints& tri
)
{
    tris_[nTris_++] = tri;
}


template<class Point, class PointRef>
inline Foam::point Foam::triangle<Point, PointRef>::planeIntersection
(
    const FixedList<scalar, 3>& d,
    const triPoints& t,
    const label negI,
    const label posI
)
{
    return (d[posI]*t[negI] - d[negI]*t[posI])/(-d[negI] + d[posI]);
}


894 895
// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //

896
template<class Point, class PointRef>
897 898
inline Foam::Istream& Foam::operator>>
(
899 900
    Istream& is,
    triangle<Point, PointRef>& t
901
)
902 903
{
    is.readBegin("triangle");
904
    is  >> t.a_ >> t.b_ >> t.c_;
905 906
    is.readEnd("triangle");

907
    is.check(FUNCTION_NAME);
908 909 910 911 912
    return is;
}


template<class Point, class PointRef>
913 914 915 916 917
inline Foam::Ostream& Foam::operator<<
(
    Ostream& os,
    const triangle<Point, PointRef>& t
)
918 919 920
{
    os  << nl
        << token::BEGIN_LIST
921 922 923
        << t.a_ << token::SPACE
        << t.b_ << token::SPACE
        << t.c_
924 925 926 927 928 929 930
        << token::END_LIST;

    return os;
}


// ************************************************************************* //