triangleI.H 20.2 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
mattijs's avatar
mattijs committed
5
    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
6
7
8
9
10
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

11
12
13
14
    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.
15
16
17
18
19
20
21

    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
22
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
23
24
25
26
27

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

#include "IOstreams.H"
#include "pointHit.H"
28
#include "mathematicalConstants.H"
29
30
31
32

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

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


template<class Point, class PointRef>
47
48
49
50
51
52
53
54
55
56
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]])
{}
57
58
59



60
61
62
63
template<class Point, class PointRef>
inline Foam::triangle<Point, PointRef>::triangle(Istream& is)
{
    is  >> *this;
64
65
66
67
68
69
}


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

template<class Point, class PointRef>
70
inline const Point& Foam::triangle<Point, PointRef>::a() const
71
72
73
74
75
{
    return a_;
}

template<class Point, class PointRef>
76
inline const Point& Foam::triangle<Point, PointRef>::b() const
77
78
79
80
81
{
    return b_;
}

template<class Point, class PointRef>
82
inline const Point& Foam::triangle<Point, PointRef>::c() const
83
84
85
86
87
88
{
    return c_;
}


template<class Point, class PointRef>
89
inline Point Foam::triangle<Point, PointRef>::centre() const
90
91
92
93
94
95
{
    return (1.0/3.0)*(a_ + b_ + c_);
}


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
106
107
108
109
{
    return 0.5*((b_ - a_)^(c_ - a_));
}


template<class Point, class PointRef>
110
inline Point Foam::triangle<Point, PointRef>::circumCentre() const
111
{
112
113
114
    scalar d1 =  (c_ - a_) & (b_ - a_);
    scalar d2 = -(c_ - b_) & (b_ - a_);
    scalar d3 =  (c_ - a_) & (c_ - b_);
115
116
117
118
119
120
121

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

    scalar c = c1 + c2 + c3;

122
123
    if (Foam::mag(c) < ROOTVSMALL)
    {
124
        // Degenerate triangle, returning centre instead of circumCentre.
125
126
127
128

        return centre();
    }

129
    return
130
131
132
133
134
135
136
    (
        ((c2 + c3)*a_ + (c3 + c1)*b_ + (c1 + c2)*c_)/(2*c)
    );
}


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

    scalar denom = d2*d3 + d3*d1 + d1*d2;

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

149
150
151
152
153
154
155
156
157
158
159
160
        return GREAT;
    }
    else
    {
        scalar a = (d1 + d2)*(d2 + d3)*(d3 + d1) / denom;

        return 0.5*Foam::sqrt(min(GREAT, max(0, a)));
    }
}


template<class Point, class PointRef>
161
inline Foam::scalar Foam::triangle<Point, PointRef>::quality() const
162
{
163
164
165
166
167
168
169
170
171
    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);
172
173
174
175
}


template<class Point, class PointRef>
176
177
178
179
inline Foam::scalar Foam::triangle<Point, PointRef>::sweptVol
(
    const triangle& t
) const
180
{
181
    return (1.0/12.0)*
182
183
184
185
    (
        ((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_)))
186
187
188
189

      + ((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_)))
190
191
192
193
    );
}


194
template<class Point, class PointRef>
195
inline Foam::tensor Foam::triangle<Point, PointRef>::inertia
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
(
    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;
}

230
231

template<class Point, class PointRef>
232
inline Point Foam::triangle<Point, PointRef>::randomPoint(Random& rndGen) const
233
234
235
236
237
238
239
240
241
242
243
244
245
246
{
    // Generating Random Points in Triangles
    // by Greg Turk
    // from "Graphics Gems", Academic Press, 1990
    // http://tog.acm.org/GraphicsGems/gems/TriPoints.c

    scalar s = rndGen.scalar01();

    scalar t = sqrt(rndGen.scalar01());

    return (1 - t)*a_ + (1 - s)*t*b_ + s*t*c_;
}


247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
template<class Point, class PointRef>
inline Point Foam::triangle<Point, PointRef>::randomPoint
(
    cachedRandom& rndGen
) const
{
    // Generating Random Points in Triangles
    // by Greg Turk
    // from "Graphics Gems", Academic Press, 1990
    // http://tog.acm.org/GraphicsGems/gems/TriPoints.c

    scalar s = rndGen.sample01<scalar>();

    scalar t = sqrt(rndGen.sample01<scalar>());

    return (1 - t)*a_ + (1 - s)*t*b_ + s*t*c_;
}


266
template<class Point, class PointRef>
267
Foam::scalar Foam::triangle<Point, PointRef>::barycentric
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
(
    const point& pt,
    List<scalar>& bary
) const
{
    // From:
    // 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)
    {
290
        // Degenerate triangle, returning 1/3 barycentric coordinates.
291
292
293
294
295
296
297
298

        bary = List<scalar>(3, 1.0/3.0);

        return denom;
    }

    bary.setSize(3);

299
300
301
    bary[1] = (d11*d20 - d01*d21)/denom;
    bary[2] = (d00*d21 - d01*d20)/denom;
    bary[0] = 1.0 - bary[1] - bary[2];
302
303
304
305
306

    return denom;
}


307
template<class Point, class PointRef>
308
inline Foam::pointHit Foam::triangle<Point, PointRef>::ray
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
(
    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
347
348
        n /= area;

349
350
351
        return ray(p, q1 - n, alg, intersection::VECTOR);
    }

352
    // Intersection point with triangle plane
353
    point pInter;
354
    // Is intersection point inside triangle
mattijs's avatar
mattijs committed
355
356
357
358
359
360
    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();
361
362
363
364
365
366
367
368
369
370
371

        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
372
    }
373

374
    // Distance to intersection point
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
    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
390
     || (alg == intersection::HALF_RAY && dist > -planarPointTol)
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
     || (
            alg == intersection::VISIBLE
         && ((q1 & normal()) < -VSMALL)
        );

    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
409
        inter.setPoint(nearestPoint(p).rawPoint());
410
411
412
413
414
415

        // 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
416

417
418
419
420
    return inter;
}


mattijs's avatar
mattijs committed
421
422
423
// From "Fast, Minimum Storage Ray/Triangle Intersection"
// Moeller/Trumbore.
template<class Point, class PointRef>
424
inline Foam::pointHit Foam::triangle<Point, PointRef>::intersection
mattijs's avatar
mattijs committed
425
426
(
    const point& orig,
427
    const vector& dir,
mattijs's avatar
mattijs committed
428
429
    const intersection::algorithm alg,
    const scalar tol
mattijs's avatar
mattijs committed
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
) 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
    pointHit intersection(false, vector::zero, GREAT, false);

    if (alg == intersection::VISIBLE)
    {
        // Culling branch
mattijs's avatar
mattijs committed
447
        if (det < ROOTVSMALL)
mattijs's avatar
mattijs committed
448
        {
mattijs's avatar
mattijs committed
449
            // Ray on wrong side of triangle. Return miss
mattijs's avatar
mattijs committed
450
451
452
453
454
455
            return intersection;
        }
    }
    else if (alg == intersection::HALF_RAY || alg == intersection::FULL_RAY)
    {
        // Non-culling branch
mattijs's avatar
mattijs committed
456
        if (det > -ROOTVSMALL && det < ROOTVSMALL)
mattijs's avatar
mattijs committed
457
        {
mattijs's avatar
mattijs committed
458
            // Ray parallel to triangle. Return miss
mattijs's avatar
mattijs committed
459
460
            return intersection;
        }
mattijs's avatar
mattijs committed
461
    }
mattijs's avatar
mattijs committed
462

mattijs's avatar
mattijs committed
463
    const scalar inv_det = 1.0 / det;
mattijs's avatar
mattijs committed
464

mattijs's avatar
mattijs committed
465
466
    /* calculate distance from a_ to ray origin */
    const vector tVec = orig-a_;
mattijs's avatar
mattijs committed
467

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

mattijs's avatar
mattijs committed
471
472
473
474
475
    if (u < -tol || u > 1.0+tol)
    {
        // return miss
        return intersection;
    }
mattijs's avatar
mattijs committed
476

mattijs's avatar
mattijs committed
477
478
479
480
481
482
483
484
485
486
    /* 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
487
    }
mattijs's avatar
mattijs committed
488
489
490
491
492

    /* 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
493
    {
mattijs's avatar
mattijs committed
494
495
        // Wrong side of orig. Return miss
        return intersection;
mattijs's avatar
mattijs committed
496
497
    }

mattijs's avatar
mattijs committed
498
499
500
501
    intersection.setHit();
    intersection.setPoint(a_ + u*edge1 + v*edge2);
    intersection.setDistance(t);

mattijs's avatar
mattijs committed
502
503
504
    return intersection;
}

505
506

template<class Point, class PointRef>
507
Foam::pointHit Foam::triangle<Point, PointRef>::nearestPointClassify
508
509
510
511
512
513
(
    const point& p,
    label& nearType,
    label& nearLabel
) const
{
514
    // Adapted from:
515
    // Real-time collision detection, Christer Ericson, 2005, p136-142
516

517
518
519
520
    // Check if P in vertex region outside A
    vector ab = b_ - a_;
    vector ac = c_ - a_;
    vector ap = p - a_;
521

522
523
    scalar d1 = ab & ap;
    scalar d2 = ac & ap;
524

525
526
527
    if (d1 <= 0.0 && d2 <= 0.0)
    {
        // barycentric coordinates (1, 0, 0)
528

529
530
531
532
        nearType = POINT;
        nearLabel = 0;
        return pointHit(false, a_, Foam::mag(a_ - p), true);
    }
533

534
535
536
537
    // Check if P in vertex region outside B
    vector bp = p - b_;
    scalar d3 = ab & bp;
    scalar d4 = ac & bp;
538

539
540
541
    if (d3 >= 0.0 && d4 <= d3)
    {
        // barycentric coordinates (0, 1, 0)
542

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

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

551
    if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0)
552
    {
553
554
        if ((d1 - d3) < ROOTVSMALL)
        {
555
            // Degenerate triangle, for d1 = d3, a_ and b_ are likely coincident
556
557
558
559
560
            nearType = POINT;
            nearLabel = 0;
            return pointHit(false, a_, Foam::mag(a_ - p), true);
        }

561
562
        // barycentric coordinates (1-v, v, 0)
        scalar v = d1/(d1 - d3);
563

564
565
566
567
        point nearPt =  a_ + v*ab;
        nearType = EDGE;
        nearLabel = 0;
        return pointHit(false, nearPt, Foam::mag(nearPt - p), true);
568
569
    }

570
571
572
573
    // Check if P in vertex region outside C
    vector cp = p - c_;
    scalar d5 = ab & cp;
    scalar d6 = ac & cp;
574

575
576
577
    if (d6 >= 0.0 && d5 <= d6)
    {
        // barycentric coordinates (0, 0, 1)
578

579
580
581
582
        nearType = POINT;
        nearLabel = 2;
        return pointHit(false, c_, Foam::mag(c_ - p), true);
    }
583

584
585
    // 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
586

587
    if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0)
mattijs's avatar
mattijs committed
588
    {
589
590
        if ((d2 - d6) < ROOTVSMALL)
        {
591
            // Degenerate triangle, for d2 = d6, a_ and c_ are likely coincident
592
593
594
595
596
            nearType = POINT;
            nearLabel = 0;
            return pointHit(false, a_, Foam::mag(a_ - p), true);
        }

597
598
599
600
601
602
603
        // 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
604
605
    }

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

609
    if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0)
610
    {
611
612
        if (((d4 - d3) + (d5 - d6)) < ROOTVSMALL)
        {
613
614
            // Degenerate triangle, for (d4 - d3) = (d6 - d5), b_ and c_ are
            // likely coincident
615
616
617
618
619
            nearType = POINT;
            nearLabel = 1;
            return pointHit(false, b_, Foam::mag(b_ - p), true);
        }

620
621
        // barycentric coordinates (0, 1-w, w)
        scalar w = (d4 - d3)/((d4 - d3) + (d5 - d6));
622

623
624
625
626
        point nearPt = b_ + w*(c_ - b_);
        nearType = EDGE;
        nearLabel = 1;
        return pointHit(false, nearPt, Foam::mag(nearPt - p), true);
627
628
    }

629
630
    // P inside face region. Compute Q through its barycentric
    // coordinates (u, v, w)
631
632
633

    if ((va + vb + vc) < ROOTVSMALL)
    {
634
635
        // Degenerate triangle, return the centre because no edge or points are
        // closest
636
637
638
639
640
641
        point nearPt = centre();
        nearType = NONE,
        nearLabel = -1;
        return pointHit(true, nearPt, Foam::mag(nearPt - p), false);
    }

642
643
644
    scalar denom = 1.0/(va + vb + vc);
    scalar v = vb * denom;
    scalar w = vc * denom;
645

646
    // = u*a + v*b + w*c, u = va*denom = 1.0 - v - w
647

648
649
650
651
    point nearPt = a_ + ab*v + ac*w;
    nearType = NONE,
    nearLabel = -1;
    return pointHit(true, nearPt, Foam::mag(nearPt - p), false);
652
653
654
}


655
template<class Point, class PointRef>
656
inline Foam::pointHit Foam::triangle<Point, PointRef>::nearestPoint
657
658
659
660
661
662
663
664
665
666
(
    const point& p
) const
{
    // Dummy labels
    label nearType = -1;
    label nearLabel = -1;

    return nearestPointClassify(p, nearType, nearLabel);
}
667
668


669
template<class Point, class PointRef>
670
inline bool Foam::triangle<Point, PointRef>::classify
671
672
673
674
675
676
677
678
679
(
    const point& p,
    label& nearType,
    label& nearLabel
) const
{
    return nearestPointClassify(p, nearType, nearLabel).hit();
}

680

mattijs's avatar
mattijs committed
681
682
683
684
685
686
687
688
689
690
691
692
693
694
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
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
737
        //label minEdgeIndex = 0;
mattijs's avatar
mattijs committed
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
        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
759
                //minEdgeIndex = 1;
mattijs's avatar
mattijs committed
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
            }
        }

        {
            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
777
                //minEdgeIndex = 2;
mattijs's avatar
mattijs committed
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
            }
        }

        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;
}


807
808
// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //

809
template<class Point, class PointRef>
810
811
inline Foam::Istream& Foam::operator>>
(
812
813
    Istream& is,
    triangle<Point, PointRef>& t
814
)
815
816
{
    is.readBegin("triangle");
817
    is  >> t.a_ >> t.b_ >> t.c_;
818
819
820
821
822
823
824
825
    is.readEnd("triangle");

    is.check("Istream& operator>>(Istream&, triangle&)");
    return is;
}


template<class Point, class PointRef>
826
827
828
829
830
inline Foam::Ostream& Foam::operator<<
(
    Ostream& os,
    const triangle<Point, PointRef>& t
)
831
832
833
{
    os  << nl
        << token::BEGIN_LIST
834
835
836
        << t.a_ << token::SPACE
        << t.b_ << token::SPACE
        << t.c_
837
838
839
840
841
842
843
        << token::END_LIST;

    return os;
}


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