BoundaryGeometry.cpp 86.5 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
//
// Declaration of the BoundaryGeometry class
//
// Copyright (c) 200x - 2020, Achim Gsell,
//                            Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved.
//
// This file is part of OPAL.
//
// OPAL 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.
//
// You should have received a copy of the GNU General Public License
16
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
17 18 19
//

//#define   ENABLE_DEBUG
20

kraus's avatar
kraus committed
21 22
#include "Structure/BoundaryGeometry.h"

23
#include <cmath>
24 25
#include <fstream>
#include <string>
26
#include <algorithm>
gsell's avatar
gsell committed
27 28

#include "H5hut.h"
29
#include <cfloat>
30

frey_m's avatar
frey_m committed
31
#include "AbstractObjects/OpalData.h"
32
#include "Algorithms/PartBunchBase.h"
gsell's avatar
gsell committed
33 34
#include "Expressions/SRefExpr.h"
#include "Elements/OpalBeamline.h"
kraus's avatar
kraus committed
35
#include "Utilities/Options.h"
kraus's avatar
kraus committed
36
#include "Utilities/OpalException.h"
37
#include <gsl/gsl_sys.h>
gsell's avatar
gsell committed
38

gsell's avatar
gsell committed
39
extern Inform* gmsg;
40

gsell's avatar
gsell committed
41
#define SQR(x) ((x)*(x))
42 43
#define PointID(triangle_id, vertex_id) Triangles_m[triangle_id][vertex_id]
#define Point(triangle_id, vertex_id)   Points_m[Triangles_m[triangle_id][vertex_id]]
44

45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
/*
  In the following namespaces various approximately floating point
  comparisons are implemented. The used implementation is selected
  via

  namespaces cmp = IMPLEMENTATION;

*/

/*
  First we define some macros for function common in all namespaces.
*/
#define FUNC_EQ(x, y) inline bool eq(double x, double y) { \
        return almost_eq(x, y);                            \
    }

#define FUNC_EQ_ZERO(x) inline bool eq_zero(double x) { \
        return almost_eq_zero(x);                       \
    }

#define FUNC_LE(x, y) inline bool le(double x, double y) { \
        if (almost_eq(x, y)) {                             \
            return true;                                   \
        }                                                  \
        return x < y;                                      \
    }

#define FUNC_LE_ZERO(x) inline bool le_zero(double x) { \
        if (almost_eq_zero(x)) {                        \
            return true;                                \
        }                                               \
        return x < 0.0;                                 \
    }

#define FUNC_LT(x, y) inline bool lt(double x, double y) { \
        if (almost_eq(x, y)) {                             \
            return false;                                  \
        }                                                  \
        return x < y;                                      \
    }

#define FUNC_LT_ZERO(x) inline bool lt_zero(double x) {      \
        if (almost_eq_zero(x)) {                             \
            return false;                                    \
        }                                                    \
        return x < 0.0;                                      \
    }

#define FUNC_GE(x, y) inline bool ge(double x, double y) {      \
        if (almost_eq(x, y)) {                                  \
            return true;                                        \
        }                                                       \
        return x > y;                                           \
    }
99

100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
#define FUNC_GE_ZERO(x) inline bool ge_zero(double x) { \
        if (almost_eq_zero(x)) {                        \
            return true;                                \
        }                                               \
        return x > 0.0;                                 \
    }

#define FUNC_GT(x, y) inline bool gt(double x, double y) { \
        if (almost_eq(x, y)) {                             \
            return false;                                  \
        }                                                  \
        return x > y;                                      \
    }

#define FUNC_GT_ZERO(x) inline bool gt_zero(double x) { \
        if (almost_eq_zero(x)) {                        \
            return false;                               \
        }                                               \
        return x > 0.0;                                 \
    }
120 121 122 123 124 125 126

namespace cmp_diff {

    /*
      Link:
      https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
    */
127
    inline bool almost_eq(double A, double B, double maxDiff = 1e-15, double maxRelDiff = DBL_EPSILON) {
128 129
        // Check if the numbers are really close -- needed
        // when comparing numbers near zero.
130
        const double diff = std::abs(A - B);
131 132 133
        if (diff <= maxDiff)
            return true;

134 135
        A = std::abs(A);
        B = std::abs(B);
136 137 138 139 140 141 142
        const double largest = (B > A) ? B : A;

        if (diff <= largest * maxRelDiff)
            return true;
        return false;
    }

143
    inline bool almost_eq_zero(double A, double maxDiff = 1e-15) {
144
        const double diff = std::abs(A);
145
        return (diff <= maxDiff);
146 147
    }

148 149 150 151 152 153 154 155 156 157
    FUNC_EQ(x, y);
    FUNC_EQ_ZERO(x);
    FUNC_LE(x, y);
    FUNC_LE_ZERO(x);
    FUNC_LT(x, y);
    FUNC_LT_ZERO(x);
    FUNC_GE(x, y);
    FUNC_GE_ZERO(x);
    FUNC_GT(x, y);
    FUNC_GT_ZERO(x);
158 159
}

160 161 162 163 164 165
namespace cmp_ulp_obsolete {
    /*
      See:
      https://www.cygnus-software.com/papers/comparingfloats/comparing_floating_point_numbers_obsolete.htm
    */
    inline bool almost_eq(double A, double B, double maxDiff = 1e-20, int maxUlps = 1000) {
166 167 168 169 170 171 172 173 174 175 176 177 178 179
        // Make sure maxUlps is non-negative and small enough that the
        // default NAN won't compare as equal to anything.
        // assert(maxUlps > 0 && maxUlps < 4 * 1024 * 1024);

        // handle NaN's
        // Note: comparing something with a NaN is always false!
        if (A != A || B != B) {
            return false;
        }

        if (std::abs(A - B) <= maxDiff) {
            return true;
        }

180 181
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wstrict-aliasing"
182
        auto aInt = *(int64_t*)&A;
183
#pragma GCC diagnostic pop
184 185 186 187 188
        // Make aInt lexicographically ordered as a twos-complement int
        if (aInt < 0) {
            aInt = 0x8000000000000000 - aInt;
        }

189 190
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wstrict-aliasing"
191
        auto bInt = *(int64_t*)&B;
192
#pragma GCC diagnostic pop
193 194 195 196 197 198 199 200 201 202 203
        // Make bInt lexicographically ordered as a twos-complement int
        if (bInt < 0) {
            bInt = 0x8000000000000000 - bInt;
        }

        if (std::abs(aInt - bInt) <= maxUlps) {
            return true;
        }
        return false;
    }

204 205 206
    inline bool almost_eq_zero(double A, double maxDiff = 1e-15) {
        // no need to handle NaN's!
        return (std::abs(A) <= maxDiff);
207
    }
208 209 210 211 212 213 214 215 216 217 218
    FUNC_EQ(x, y);
    FUNC_EQ_ZERO(x);
    FUNC_LE(x, y);
    FUNC_LE_ZERO(x);
    FUNC_LT(x, y);
    FUNC_LT_ZERO(x);
    FUNC_GE(x, y);
    FUNC_GE_ZERO(x);
    FUNC_GT(x, y);
    FUNC_GT_ZERO(x);
}
219

220 221 222 223 224 225 226 227 228
namespace cmp_ulp {
    /*
      See:
      https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
    */
    inline bool almost_eq(double A, double B, double maxDiff = 1e-20, int maxUlps = 1000) {
        // handle NaN's
        // Note: comparing something with a NaN is always false!
        if (A != A || B != B) {
229 230 231
            return false;
        }

232 233
        // Check if the numbers are really close -- needed
        // when comparing numbers near zero.
234
        if (std::abs(A - B) <= maxDiff)
235 236
            return true;

237 238 239 240 241 242 243 244
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wstrict-aliasing"
        auto aInt = *(int64_t*)&A;
        auto bInt = *(int64_t*)&B;
#pragma GCC diagnostic pop

        // Different signs means they do not match.
        // Note: a negative floating point number is also negative as integer.
245
        if ( std::signbit(aInt) != std::signbit(bInt))
246
            return false;
247 248 249 250 251 252
 
        // Find the difference in ULPs.
        return (std::abs(aInt - bInt) <= maxUlps);
    }

    inline bool almost_eq_zero(double A, double maxDiff = 1e-15) {
253
        return (std::abs(A) <= maxDiff);
254
    }
255 256 257 258 259 260 261 262 263 264
    FUNC_EQ(x, y);
    FUNC_EQ_ZERO(x);
    FUNC_LE(x, y);
    FUNC_LE_ZERO(x);
    FUNC_LT(x, y);
    FUNC_LT_ZERO(x);
    FUNC_GE(x, y);
    FUNC_GE_ZERO(x);
    FUNC_GT(x, y);
    FUNC_GT_ZERO(x);
265 266 267
}

namespace cmp = cmp_ulp;
268 269
/*

gsell's avatar
gsell committed
270
  Some
271 272
   _   _      _
  | | | | ___| |_ __   ___ _ __
gsell's avatar
gsell committed
273
  | |_| |/ _ \ | '_ \ / _ \ '__|
274 275 276
  |  _  |  __/ | |_) |  __/ |
  |_| |_|\___|_| .__/ \___|_|
             |_|
277

gsell's avatar
gsell committed
278
  functions
279
 */
snuverink_j's avatar
snuverink_j committed
280
namespace {
gsell's avatar
gsell committed
281 282
struct VectorLessX {
    bool operator() (Vector_t x1, Vector_t x2) {
283
        return cmp::lt (x1(0), x2(0));
gsell's avatar
gsell committed
284 285
    }
};
286

gsell's avatar
gsell committed
287 288
struct VectorLessY {
    bool operator() (Vector_t x1, Vector_t x2) {
289
        return cmp::lt(x1(1), x2 (1));
gsell's avatar
gsell committed
290 291
    }
};
292

gsell's avatar
gsell committed
293 294
struct VectorLessZ {
    bool operator() (Vector_t x1, Vector_t x2) {
295
        return cmp::lt(x1(2), x2(2));
296
    }
gsell's avatar
gsell committed
297 298 299 300 301
};

/**
   Calculate the maximum of coordinates of geometry,i.e the maximum of X,Y,Z
 */
gsell's avatar
gsell committed
302
Vector_t get_max_extent (std::vector<Vector_t>& coords) {
gsell's avatar
gsell committed
303 304 305 306 307 308
    const Vector_t x = *max_element (
        coords.begin (), coords.end (), VectorLessX ());
    const Vector_t y = *max_element (
        coords.begin (), coords.end (), VectorLessY ());
    const Vector_t z = *max_element (
        coords.begin (), coords.end (), VectorLessZ ());
309
    return Vector_t (x(0), y(1), z(2));
310 311
}

snuverink_j's avatar
snuverink_j committed
312

gsell's avatar
gsell committed
313 314 315
/*
   Compute the minimum of coordinates of geometry, i.e the minimum of X,Y,Z
 */
gsell's avatar
gsell committed
316
Vector_t get_min_extent (std::vector<Vector_t>& coords) {
gsell's avatar
gsell committed
317 318 319 320 321 322
    const Vector_t x = *min_element (
        coords.begin (), coords.end (), VectorLessX ());
    const Vector_t y = *min_element (
        coords.begin (), coords.end (), VectorLessY ());
    const Vector_t z = *min_element (
        coords.begin (), coords.end (), VectorLessZ ());
323
    return Vector_t (x(0), y(1), z(2));
gsell's avatar
gsell committed
324
}
325

326 327 328 329
/*
  write legacy VTK file of voxel mesh
*/
static void write_voxel_mesh (
330 331 332 333
    const std::unordered_map< int, std::unordered_set<int> >& ids,
    const Vector_t& hr_m,
    const Vektor<int,3>& nr,
    const Vector_t& origin
334 335 336 337
    ) {
    /*----------------------------------------------------------------------*/
    const size_t numpoints = 8 * ids.size ();
    std::ofstream of;
338 339 340 341 342
    std::string fname = Util::combineFilePath({
        OpalData::getInstance()->getAuxiliaryOutputDirectory(),
        "testBBox.vtk"
    });
    of.open (fname);
343 344
    assert (of.is_open ());
    of.precision (6);
345

346
    of << "# vtk DataFile Version 2.0" << std::endl;
347
    of << "generated using BoundaryGeometry::computeMeshVoxelization"
348 349 350 351
       << std::endl;
    of << "ASCII" << std::endl << std::endl;
    of << "DATASET UNSTRUCTURED_GRID" << std::endl;
    of << "POINTS " << numpoints << " float" << std::endl;
352 353

    const auto nr0_times_nr1 = nr[0] * nr[1];
354 355
    for (auto& elem: ids) {
        auto id = elem.first;
356 357
        int k = (id - 1) / nr0_times_nr1;
        int rest = (id - 1) % nr0_times_nr1;
358
        int j = rest / nr[0];
359
        int i = rest % nr[0];
360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375

        Vector_t P;
        P[0] = i * hr_m[0] + origin[0];
        P[1] = j * hr_m[1] + origin[1];
        P[2] = k * hr_m[2] + origin[2];

        of << P[0]           << " " << P[1]           << " " << P[2]           << std::endl;
        of << P[0] + hr_m[0] << " " << P[1]           << " " << P[2]           << std::endl;
        of << P[0]           << " " << P[1] + hr_m[1] << " " << P[2]           << std::endl;
        of << P[0] + hr_m[0] << " " << P[1] + hr_m[1] << " " << P[2]           << std::endl;
        of << P[0]           << " " << P[1]           << " " << P[2] + hr_m[2] << std::endl;
        of << P[0] + hr_m[0] << " " << P[1]           << " " << P[2] + hr_m[2] << std::endl;
        of << P[0]           << " " << P[1] + hr_m[1] << " " << P[2] + hr_m[2] << std::endl;
        of << P[0] + hr_m[0] << " " << P[1] + hr_m[1] << " " << P[2] + hr_m[2] << std::endl;
    }
    of << std::endl;
376 377 378
    const auto num_cells = ids.size ();
    of << "CELLS " << num_cells << " " << 9 * num_cells << std::endl;
    for (size_t i = 0; i < num_cells; i++)
379 380 381
        of << "8 "
           << 8 * i << " " << 8 * i + 1 << " " << 8 * i + 2 << " " << 8 * i + 3 << " "
           << 8 * i + 4 << " " << 8 * i + 5 << " " << 8 * i + 6 << " " << 8 * i + 7 << std::endl;
382 383
    of << "CELL_TYPES " << num_cells << std::endl;
    for (size_t i = 0; i <  num_cells; i++)
384
        of << "11" << std::endl;
385
    of << "CELL_DATA " << num_cells << std::endl;
386 387
    of << "SCALARS " << "cell_attribute_data" << " float " << "1" << std::endl;
    of << "LOOKUP_TABLE " << "default" << std::endl;
388
    for (size_t i = 0; i <  num_cells; i++)
389 390 391
        of << (float)(i) << std::endl;
    of << std::endl;
    of << "COLOR_SCALARS " << "BBoxColor " << 4 << std::endl;
392
    for (size_t i = 0; i < num_cells; i++) {
393 394 395 396
        of << "1.0" << " 1.0 " << "0.0 " << "1.0" << std::endl;
    }
    of << std::endl;
}
snuverink_j's avatar
snuverink_j committed
397
}
398

399 400 401 402 403 404 405 406 407
/*___________________________________________________________________________

  Triangle-cube intersection test.

  See:
  http://tog.acm.org/resources/GraphicsGems/gemsiii/triangleCube.c

 */

408 409 410 411

#define INSIDE 0
#define OUTSIDE 1

412 413 414 415 416 417 418 419
class Triangle {
public:
    Triangle () { }
    Triangle (const Vector_t& v1, const Vector_t& v2, const Vector_t& v3) {
        pts[0] = v1;
        pts[1] = v2;
        pts[2] = v3;
    }
420

421 422 423
    inline const Vector_t& v1() const {
        return pts[0];
    }
gsell's avatar
gsell committed
424
    inline double v1(int i) const {
425 426 427 428 429
        return pts[0][i];
    }
    inline const Vector_t& v2() const {
        return pts[1];
    }
gsell's avatar
gsell committed
430
    inline double v2(int i) const {
431 432 433 434 435
        return pts[1][i];
    }
    inline const Vector_t& v3() const {
        return pts[2];
    }
gsell's avatar
gsell committed
436
    inline double v3(int i) const {
437 438
        return pts[2][i];
    }
439

440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461

    inline void scale (
        const Vector_t& scaleby,
        const Vector_t& shiftby
        ) {
        pts[0][0] *= scaleby[0];
        pts[0][1] *= scaleby[1];
        pts[0][2] *= scaleby[2];
        pts[1][0] *= scaleby[0];
        pts[1][1] *= scaleby[1];
        pts[1][2] *= scaleby[2];
        pts[2][0] *= scaleby[0];
        pts[2][1] *= scaleby[1];
        pts[2][2] *= scaleby[2];
        pts[0] -= shiftby;
        pts[1] -= shiftby;
        pts[2] -= shiftby;
    }


    Vector_t pts[3];
};
462 463 464 465 466 467 468

/*___________________________________________________________________________*/

/* Which of the six face-plane(s) is point P outside of? */

static inline int
face_plane (
469
    const Vector_t& p
470
    ) {
471
    int outcode_fcmp = 0;
472

473 474 475 476 477 478
    if (cmp::gt(p[0],  0.5)) outcode_fcmp |= 0x01;
    if (cmp::lt(p[0], -0.5)) outcode_fcmp |= 0x02;
    if (cmp::gt(p[1],  0.5)) outcode_fcmp |= 0x04;
    if (cmp::lt(p[1], -0.5)) outcode_fcmp |= 0x08;
    if (cmp::gt(p[2],  0.5)) outcode_fcmp |= 0x10;
    if (cmp::lt(p[2], -0.5)) outcode_fcmp |= 0x20;
479 480

    return(outcode_fcmp);
481 482 483 484 485 486 487 488
}

/*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . */

/* Which of the twelve edge plane(s) is point P outside of? */

static inline int
bevel_2d (
489
    const Vector_t& p
490
    ) {
491 492
    int outcode_fcmp = 0;

493 494 495 496 497 498 499 500 501 502 503 504
    if (cmp::gt( p[0] + p[1], 1.0)) outcode_fcmp |= 0x001;
    if (cmp::gt( p[0] - p[1], 1.0)) outcode_fcmp |= 0x002;
    if (cmp::gt(-p[0] + p[1], 1.0)) outcode_fcmp |= 0x004;
    if (cmp::gt(-p[0] - p[1], 1.0)) outcode_fcmp |= 0x008;
    if (cmp::gt( p[0] + p[2], 1.0)) outcode_fcmp |= 0x010;
    if (cmp::gt( p[0] - p[2], 1.0)) outcode_fcmp |= 0x020;
    if (cmp::gt(-p[0] + p[2], 1.0)) outcode_fcmp |= 0x040;
    if (cmp::gt(-p[0] - p[2], 1.0)) outcode_fcmp |= 0x080;
    if (cmp::gt( p[1] + p[2], 1.0)) outcode_fcmp |= 0x100;
    if (cmp::gt( p[1] - p[2], 1.0)) outcode_fcmp |= 0x200;
    if (cmp::gt(-p[1] + p[2], 1.0)) outcode_fcmp |= 0x400;
    if (cmp::gt(-p[1] - p[2], 1.0)) outcode_fcmp |= 0x800;
505

506
    return(outcode_fcmp);
507 508
}

509
/*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
510 511 512 513 514

  Which of the eight corner plane(s) is point P outside of?
*/
static inline int
bevel_3d (
515
    const Vector_t& p
516
    ) {
517 518
    int outcode_fcmp = 0;

519 520 521 522 523 524 525 526
    if (cmp::gt( p[0] + p[1] + p[2], 1.5)) outcode_fcmp |= 0x01;
    if (cmp::gt( p[0] + p[1] - p[2], 1.5)) outcode_fcmp |= 0x02;
    if (cmp::gt( p[0] - p[1] + p[2], 1.5)) outcode_fcmp |= 0x04;
    if (cmp::gt( p[0] - p[1] - p[2], 1.5)) outcode_fcmp |= 0x08;
    if (cmp::gt(-p[0] + p[1] + p[2], 1.5)) outcode_fcmp |= 0x10;
    if (cmp::gt(-p[0] + p[1] - p[2], 1.5)) outcode_fcmp |= 0x20;
    if (cmp::gt(-p[0] - p[1] + p[2], 1.5)) outcode_fcmp |= 0x40;
    if (cmp::gt(-p[0] - p[1] - p[2], 1.5)) outcode_fcmp |= 0x80;
527

528
    return(outcode_fcmp);
529 530 531 532 533 534 535 536 537 538 539
}

/*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Test the point "alpha" of the way from P1 to P2
  See if it is on a face of the cube
  Consider only faces in "mask"
*/

static inline int
check_point (
540 541
    const Vector_t& p1,
    const Vector_t& p2,
542 543 544 545 546
    const double alpha,
    const int mask
    ) {
    Vector_t plane_point;

547 548 549 550 551 552
#define LERP(a, b, t) (a + t * (b - a))
    // with C++20 we can use: std::lerp(a, b, t)
    plane_point[0] = LERP(p1[0], p2[0], alpha);
    plane_point[1] = LERP(p1[1], p2[1], alpha);
    plane_point[2] = LERP(p1[2], p2[2], alpha);
#undef LERP
553 554 555 556 557 558 559 560 561 562 563 564
    return(face_plane(plane_point) & mask);
}

/*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Compute intersection of P1 --> P2 line segment with face planes
  Then test intersection point to see if it is on cube face
  Consider only face planes in "outcode_diff"
  Note: Zero bits in "outcode_diff" means face line is outside of
*/
static inline int
check_line (
565 566
    const Vector_t& p1,
    const Vector_t& p2,
567 568 569 570 571 572
    const int outcode_diff
    ) {
    if ((0x01 & outcode_diff) != 0)
        if (check_point(p1,p2,( .5-p1[0])/(p2[0]-p1[0]),0x3e) == INSIDE) return(INSIDE);
    if ((0x02 & outcode_diff) != 0)
        if (check_point(p1,p2,(-.5-p1[0])/(p2[0]-p1[0]),0x3d) == INSIDE) return(INSIDE);
573
    if ((0x04 & outcode_diff) != 0)
574
        if (check_point(p1,p2,( .5-p1[1])/(p2[1]-p1[1]),0x3b) == INSIDE) return(INSIDE);
575
    if ((0x08 & outcode_diff) != 0)
576
        if (check_point(p1,p2,(-.5-p1[1])/(p2[1]-p1[1]),0x37) == INSIDE) return(INSIDE);
577
    if ((0x10 & outcode_diff) != 0)
578
        if (check_point(p1,p2,( .5-p1[2])/(p2[2]-p1[2]),0x2f) == INSIDE) return(INSIDE);
579
    if ((0x20 & outcode_diff) != 0)
580 581 582 583 584 585 586 587
        if (check_point(p1,p2,(-.5-p1[2])/(p2[2]-p1[2]),0x1f) == INSIDE) return(INSIDE);
    return(OUTSIDE);
}

/*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Test if 3D point is inside 3D triangle
*/
588
#define EPS 10e-10
589 590 591 592 593 594 595 596
static inline int
SIGN3 (
    Vector_t A
    ) {
    return ((A[0] < EPS) ? 4 : 0 | (A[0] > -EPS) ? 32 : 0 |
            (A[1] < EPS) ? 2 : 0 | (A[1] > -EPS) ? 16 : 0 |
            (A[2] < EPS) ? 1 : 0 | (A[2] > -EPS) ? 8 : 0);
}
597
#undef EPS
598

599 600
static int
point_triangle_intersection (
601 602
    const Vector_t& p,
    const Triangle& t
603 604 605 606 607
    ) {
    /*
      First, a quick bounding-box test:
      If P is outside triangle bbox, there cannot be an intersection.
    */
608 609 610 611 612 613
    if (cmp::gt(p[0], std::max({t.v1(0), t.v2(0), t.v3(0)}))) return(OUTSIDE);
    if (cmp::gt(p[1], std::max({t.v1(1), t.v2(1), t.v3(1)}))) return(OUTSIDE);
    if (cmp::gt(p[2], std::max({t.v1(2), t.v2(2), t.v3(2)}))) return(OUTSIDE);
    if (cmp::lt(p[0], std::min({t.v1(0), t.v2(0), t.v3(0)}))) return(OUTSIDE);
    if (cmp::lt(p[1], std::min({t.v1(1), t.v2(1), t.v3(1)}))) return(OUTSIDE);
    if (cmp::lt(p[2], std::min({t.v1(2), t.v2(2), t.v3(2)}))) return(OUTSIDE);
614

615 616 617 618 619
    /*
      For each triangle side, make a vector out of it by subtracting vertexes;
      make another vector from one vertex to point P.
      The crossproduct of these two vectors is orthogonal to both and the
      signs of its X,Y,Z components indicate whether P was to the inside or
620
      to the outside of this triangle side.
621
    */
622 623
    const Vector_t vect12 = t.v1() - t.v2();
    const Vector_t vect1h = t.v1() - p;
624
    const Vector_t cross12_1p = cross (vect12, vect1h);
625
    const int sign12 = SIGN3(cross12_1p);      /* Extract X,Y,Z signs as 0..7 or 0...63 integer */
626

627 628
    const Vector_t vect23 = t.v2() - t.v3();
    const Vector_t vect2h = t.v2() - p;
629
    const Vector_t cross23_2p = cross (vect23, vect2h);
630
    const int sign23 = SIGN3(cross23_2p);
631

632 633
    const Vector_t vect31 = t.v3() - t.v1();
    const Vector_t vect3h = t.v3() - p;
634
    const Vector_t cross31_3p = cross (vect31, vect3h);
635
    const int sign31 = SIGN3(cross31_3p);
636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654

    /*
      If all three crossproduct vectors agree in their component signs,
      then the point must be inside all three.
      P cannot be OUTSIDE all three sides simultaneously.
    */
    return ((sign12 & sign23 & sign31) == 0) ? OUTSIDE : INSIDE;
}


/*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  This is the main algorithm procedure.
  Triangle t is compared with a unit cube,
  centered on the origin.
  It returns INSIDE (0) or OUTSIDE(1) if t
  intersects or does not intersect the cube.
*/
static int
655 656
triangle_intersects_cube (
    const Triangle& t
657 658 659 660 661 662 663 664 665
    ) {
    int v1_test;
    int v2_test;
    int v3_test;

    /*
      First compare all three vertexes with all six face-planes
      If any vertex is inside the cube, return immediately!
    */
666 667 668
    if ((v1_test = face_plane(t.v1())) == INSIDE) return(INSIDE);
    if ((v2_test = face_plane(t.v2())) == INSIDE) return(INSIDE);
    if ((v3_test = face_plane(t.v3())) == INSIDE) return(INSIDE);
669 670 671 672 673 674 675 676 677 678

   /*
     If all three vertexes were outside of one or more face-planes,
     return immediately with a trivial rejection!
   */
   if ((v1_test & v2_test & v3_test) != 0) return(OUTSIDE);

   /*
     Now do the same trivial rejection test for the 12 edge planes
   */
679 680
   v1_test |= bevel_2d(t.v1()) << 8;
   v2_test |= bevel_2d(t.v2()) << 8;
681
   v3_test |= bevel_2d(t.v3()) << 8;
682
   if ((v1_test & v2_test & v3_test) != 0) return(OUTSIDE);
683 684 685 686

   /*
     Now do the same trivial rejection test for the 8 corner planes
   */
687 688 689 690
   v1_test |= bevel_3d(t.v1()) << 24;
   v2_test |= bevel_3d(t.v2()) << 24;
   v3_test |= bevel_3d(t.v3()) << 24;
   if ((v1_test & v2_test & v3_test) != 0) return(OUTSIDE);
691 692 693 694 695 696 697 698 699 700

   /*
     If vertex 1 and 2, as a pair, cannot be trivially rejected
     by the above tests, then see if the v1-->v2 triangle edge
     intersects the cube.  Do the same for v1-->v3 and v2-->v3./
     Pass to the intersection algorithm the "OR" of the outcode
     bits, so that only those cube faces which are spanned by
     each triangle edge need be tested.
   */
   if ((v1_test & v2_test) == 0)
701
       if (check_line (t.v1(), t.v2(), v1_test|v2_test) == INSIDE) return(INSIDE);
702
   if ((v1_test & v3_test) == 0)
703
       if (check_line (t.v1(), t.v3(), v1_test|v3_test) == INSIDE) return(INSIDE);
704
   if ((v2_test & v3_test) == 0)
705
       if (check_line (t.v2(), t.v3(), v2_test|v3_test) == INSIDE) return(INSIDE);
706 707 708 709 710 711 712 713 714 715 716

   /*
     By now, we know that the triangle is not off to any side,
     and that its sides do not penetrate the cube.  We must now
     test for the cube intersecting the interior of the triangle.
     We do this by looking for intersections between the cube
     diagonals and the triangle...first finding the intersection
     of the four diagonals with the plane of the triangle, and
     then if that intersection is inside the cube, pursuing
     whether the intersection point is inside the triangle itself.

717
     To find plane of the triangle, first perform crossproduct on
718
     two triangle side vectors to compute the normal vector.
719 720 721
   */
   Vector_t vect12 = t.v1() - t.v2();
   Vector_t vect13 = t.v1() - t.v3();
722 723 724 725 726 727 728 729 730 731 732 733
   Vector_t norm = cross (vect12, vect13);

   /*
     The normal vector "norm" X,Y,Z components are the coefficients
     of the triangles AX + BY + CZ + D = 0 plane equation.  If we
     solve the plane equation for X=Y=Z (a diagonal), we get
     -D/(A+B+C) as a metric of the distance from cube center to the
     diagonal/plane intersection.  If this is between -0.5 and 0.5,
     the intersection is inside the cube.  If so, we continue by
     doing a point/triangle intersection.
     Do this for all four diagonals.
   */
734
   double d = norm[0] * t.v1(0) + norm[1] * t.v1(1) + norm[2] * t.v1(2);
735 736 737 738 739

   /*
     if one of the diagonals is parallel to the plane, the other will
     intersect the plane
   */
740 741
   double denom = norm[0] + norm[1] + norm[2];
   if (cmp::eq_zero(std::abs(denom)) == false) {
742
       /* skip parallel diagonals to the plane; division by 0 can occure */
743
       Vector_t hitpp = d / denom;
744 745 746
       if (cmp::le(std::abs(hitpp[0]), 0.5))
           if (point_triangle_intersection(hitpp,t) == INSIDE)
               return(INSIDE);
747
   }
748 749
   denom = norm[0] + norm[1] - norm[2];
   if (cmp::eq_zero(std::abs(denom)) == false) {
750 751
       Vector_t hitpn;
       hitpn[2] = -(hitpn[0] = hitpn[1] = d / denom);
752 753 754
       if (cmp::le(std::abs(hitpn[0]), 0.5))
           if (point_triangle_intersection(hitpn,t) == INSIDE)
               return(INSIDE);
755
   }
756 757
   denom = norm[0] - norm[1] + norm[2];
   if (cmp::eq_zero(std::abs(denom)) == false) {
758 759
       Vector_t hitnp;
       hitnp[1] = -(hitnp[0] = hitnp[2] = d / denom);
760 761 762
       if (cmp::le(std::abs(hitnp[0]), 0.5))
           if (point_triangle_intersection(hitnp,t) == INSIDE)
               return(INSIDE);
763
   }
764 765
   denom = norm[0] - norm[1] - norm[2];
   if (cmp::eq_zero(std::abs(denom)) == false) {
766 767
       Vector_t hitnn;
       hitnn[1] = hitnn[2] = -(hitnn[0] = d / denom);
768 769 770
       if (cmp::le(std::abs(hitnn[0]), 0.5))
           if (point_triangle_intersection(hitnn,t) == INSIDE)
               return(INSIDE);
771
   }
772

773 774 775 776 777 778 779
   /*
     No edge touched the cube; no cube diagonal touched the triangle.
     We're done...there was no intersection.
   */
   return(OUTSIDE);
}

780 781 782 783 784 785 786
/*
 * Ray class, for use with the optimized ray-box intersection test
 * described in:
 *
 *      Amy Williams, Steve Barrus, R. Keith Morley, and Peter Shirley
 *      "An Efficient and Robust Ray-Box Intersection Algorithm"
 *      Journal of graphics tools, 10(1):49-54, 2005
787
 *
788 789 790 791
 */

class Ray {
public:
792 793 794 795 796 797 798 799 800 801 802 803 804 805 806
    Ray () { }
    Ray (Vector_t o, Vector_t d) {
        origin = o;
        direction = d;
        inv_direction = Vector_t (1/d[0], 1/d[1], 1/d[2]);
        sign[0] = (inv_direction[0] < 0);
        sign[1] = (inv_direction[1] < 0);
        sign[2] = (inv_direction[2] < 0);
    }
    Ray(const Ray &r) {
        origin = r.origin;
        direction = r.direction;
        inv_direction = r.inv_direction;
        sign[0] = r.sign[0]; sign[1] = r.sign[1]; sign[2] = r.sign[2];
    }
807
    const Ray &operator=(const Ray& a) = delete;
808

809 810 811 812
    Vector_t origin;
    Vector_t direction;
    Vector_t inv_direction;
    int sign[3];
813 814
};

815

816 817 818 819 820 821 822 823 824 825
/*
 * Axis-aligned bounding box class, for use with the optimized ray-box
 * intersection test described in:
 *
 *      Amy Williams, Steve Barrus, R. Keith Morley, and Peter Shirley
 *      "An Efficient and Robust Ray-Box Intersection Algorithm"
 *      Journal of graphics tools, 10(1):49-54, 2005
 *
 */

826
class Voxel {
827
public:
828 829 830 831 832 833 834 835 836 837 838 839 840 841
    Voxel () { }
    Voxel (const Vector_t &min, const Vector_t &max) {
        pts[0] = min;
        pts[1] = max;
    }
    inline void scale (
        const Vector_t& scale
        ) {
        pts[0][0] *= scale[0];
        pts[0][1] *= scale[1];
        pts[0][2] *= scale[2];
        pts[1][0] *= scale[0];
        pts[1][1] *= scale[1];
        pts[1][2] *= scale[2];
842
    }
843

844 845 846
    // (t0, t1) is the interval for valid hits
    bool intersect (
        const Ray& r,
847 848
        double& tmin,       // tmin and tmax are unchanged, if there is
        double& tmax        // no intersection
849
        ) const {
kraus's avatar
kraus committed
850 851 852 853
        double tmin_ = (pts[r.sign[0]][0]   - r.origin[0]) * r.inv_direction[0];
        double tmax_ = (pts[1-r.sign[0]][0] - r.origin[0]) * r.inv_direction[0];
        const double tymin = (pts[r.sign[1]][1]   - r.origin[1]) * r.inv_direction[1];
        const double tymax = (pts[1-r.sign[1]][1] - r.origin[1]) * r.inv_direction[1];
854
        if ( cmp::gt(tmin_, tymax) || cmp::gt(tymin, tmax_) )
855
            return 0;       // no intersection
856
        if (cmp::gt(tymin, tmin_))
857
            tmin_ = tymin;
858
        if (cmp::lt(tymax, tmax_))
859
            tmax_ = tymax;
kraus's avatar
kraus committed
860 861
        const double tzmin = (pts[r.sign[2]][2]   - r.origin[2]) * r.inv_direction[2];
        const double tzmax = (pts[1-r.sign[2]][2] - r.origin[2]) * r.inv_direction[2];
862
        if ( cmp::gt(tmin_, tzmax) || cmp::gt(tzmin, tmax_) )
863
            return 0;       // no intersection
864
        if (cmp::gt(tzmin, tmin_))
kraus's avatar
kraus committed
865 866
            tmin_ = tzmin;
        tmin = tmin_;
867
        if (cmp::lt(tzmax, tmax_))
kraus's avatar
kraus committed
868 869
            tmax_ = tzmax;
        tmax = tmax_;
870
        return cmp::ge_zero(tmax);
871 872
    }

873 874
    inline bool intersect (
        const Ray& r
875
        ) const {
876 877 878
        double tmin = 0.0;
        double tmax = 0.0;
        return intersect(r, tmin, tmax);
879
    }
880 881 882

    inline int intersect (
        const Triangle& t
883
        ) const {
884 885
        Voxel v_ = *this;
        Triangle t_ = t;
886
        const Vector_t scaleby = 1.0 / v_.extent();
887 888 889
        v_.scale (scaleby);
        t_.scale (scaleby , v_.pts[0] + 0.5);
        return triangle_intersects_cube (t_);
890
    }
891

gsell's avatar
gsell committed
892
    inline Vector_t extent () const {
893 894 895 896 897 898 899
        return (pts[1] - pts[0]);
    }

    inline bool isInside  (
        const Vector_t& P
        ) const {
        return (
900 901 902 903 904 905
                cmp::ge(P[0], pts[0][0])
                && cmp::ge(P[1], pts[0][1])
                && cmp::ge(P[2], pts[0][2])
                && cmp::le(P[0], pts[1][0])
                && cmp::le(P[1], pts[1][1])
                && cmp::le(P[2], pts[1][2]));
906 907
    }

908 909 910
    Vector_t pts[2];
};

911 912 913 914 915 916
static inline Vector_t normalVector (
    const Vector_t& A,
    const Vector_t& B,
    const Vector_t& C
    ) {
    const Vector_t N = cross (B - A, C - A);
917
    const double magnitude = std::sqrt (SQR (N (0)) + SQR (N (1)) + SQR (N (2)));
918
    PAssert (cmp::gt_zero(magnitude)); // in case we have degenerated triangles
919 920 921 922 923 924 925 926 927 928 929
    return N / magnitude;
}

// Calculate the area of triangle given by id.
static inline double computeArea (
    const Vector_t& A,
    const Vector_t& B,
    const Vector_t& C
    ) {
    const Vector_t AB = A - B;
    const Vector_t AC = C - A;
930
    return(0.5 * std::sqrt (dot (AB, AB) * dot (AC, AC) - dot (AB, AC) * dot (AB, AC)));
931 932
}

933

gsell's avatar
gsell committed
934
/*
935 936
  ____                           _
 / ___| ___  ___  _ __ ___   ___| |_ _ __ _   _
gsell's avatar
gsell committed
937 938 939 940 941
| |  _ / _ \/ _ \| '_ ` _ \ / _ \ __| '__| | | |
| |_| |  __/ (_) | | | | | |  __/ |_| |  | |_| |
 \____|\___|\___/|_| |_| |_|\___|\__|_|   \__, |
                                          |___/
*/
942

gsell's avatar
gsell committed
943
BoundaryGeometry::BoundaryGeometry() :
944
    Definition (
945
        SIZE, "GEOMETRY", "The \"GEOMETRY\" statement defines the beam pipe geometry.") {
946

gsell's avatar
gsell committed
947
    itsAttr[FGEOM] = Attributes::makeString
948
        ("FGEOM",
949
         "Specifies the geometry file [H5hut]",
950
         "");
gsell's avatar
gsell committed
951 952

    itsAttr[TOPO] = Attributes::makeString
953
        ("TOPO",
frey_m's avatar
frey_m committed
954
         "RECTANGULAR, BOXCORNER, ELLIPTIC if FGEOM is selected topo is over-written ",
955
         "ELLIPTIC");
gsell's avatar
gsell committed
956

957 958
    itsAttr[LENGTH] = Attributes::makeReal
        ("LENGTH",
959 960
         "Specifies the length of a tube shaped elliptic beam pipe [m]",
         1.0);
gsell's avatar
gsell committed
961 962

    itsAttr[S] = Attributes::makeReal
963 964 965
        ("S",
         "Specifies the start of a tube shaped elliptic beam pipe [m]",
         0.0);
gsell's avatar
gsell committed
966 967

    itsAttr[A] = Attributes::makeReal
968 969 970
        ("A",
         "Specifies the major semi-axis of a tube shaped elliptic beam pipe [m]",
         0.025);
gsell's avatar
gsell committed
971 972

    itsAttr[B] = Attributes::makeReal
973 974 975
        ("B",
         "Specifies the major semi-axis of a tube shaped elliptic beam pipe [m]",
         0.025);
gsell's avatar
gsell committed
976 977

    itsAttr[L1] = Attributes::makeReal
978
        ("L1",
979
         "In case of BOXCORNER Specifies first part with height == B [m]",
980
         0.5);
gsell's avatar
gsell committed
981 982

    itsAttr[L2] = Attributes::makeReal
983
        ("L2",
984
         "In case of BOXCORNER Specifies first second with height == B-C [m]",
985
         0.2);
gsell's avatar
gsell committed
986 987

    itsAttr[C] = Attributes::makeReal
988
        ("C",
989
         "In case of BOXCORNER Specifies height of corner C [m]",
990
         0.01);
gsell's avatar
gsell committed
991 992

    itsAttr[XYZSCALE] = Attributes::makeReal
993 994 995
        ("XYZSCALE",
         "Multiplicative scaling factor for coordinates ",
         1.0);
gsell's avatar
gsell committed
996

997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011
    itsAttr[XSCALE] = Attributes::makeReal
        ("XSCALE",
         "Multiplicative scaling factor for X coordinates ",
         1.0);

    itsAttr[YSCALE] = Attributes::makeReal
        ("YSCALE",
         "Multiplicative scaling factor for Y coordinates ",
         1.0);

    itsAttr[ZSCALE] = Attributes::makeReal
        ("ZSCALE",
         "Multiplicative scaling factor for Z coordinates ",
         1.0);

gsell's avatar
gsell committed
1012
    itsAttr[ZSHIFT] = Attributes::makeReal
1013 1014 1015
        ("ZSHIFT",
         "Shift in z direction",
         0.0);
gsell's avatar
gsell committed
1016

1017 1018 1019
    itsAttr[INSIDEPOINT] = Attributes::makeRealArray
        ("INSIDEPOINT", "A point inside the geometry");

1020 1021
    registerOwnership(AttributeHandler::STATEMENT);

1022
    BoundaryGeometry* defGeometry = clone ("UNNAMED_GEOMETRY");
gsell's avatar
gsell committed
1023 1024
    defGeometry->builtin = true;

gsell's avatar
gsell committed
1025 1026 1027 1028 1029 1030
    Tinitialize_m =   IpplTimings::getTimer ("Initialize geometry");
    TisInside_m =     IpplTimings::getTimer ("Inside test");
    TfastIsInside_m = IpplTimings::getTimer ("Fast inside test");
    TRayTrace_m =     IpplTimings::getTimer ("Ray tracing");
    TPartInside_m =   IpplTimings::getTimer ("Particle Inside");

1031
    h5FileName_m = Attributes::getString (itsAttr[FGEOM]);
gsell's avatar
gsell committed
1032
    try {
1033 1034 1035
        defGeometry->update ();
        OpalData::getInstance ()->define (defGeometry);
    } catch (...) {
gsell's avatar
gsell committed
1036 1037
        delete defGeometry;
    }
1038
    gsl_rng_env_setup();
1039
    randGen_m = gsl_rng_alloc(gsl_rng_default);
1040

1041 1042
    if (!h5FileName_m.empty ())
        initialize ();
1043

gsell's avatar
gsell committed
1044 1045 1046
}

BoundaryGeometry::BoundaryGeometry(
1047
    const std::string& name,
1048
    BoundaryGeometry* parent
1049
    ) : Definition (name, parent) {
1050
    gsl_rng_env_setup();
1051
    randGen_m = gsl_rng_alloc(gsl_rng_default);
1052

gsell's avatar
gsell committed
1053 1054 1055 1056 1057
    Tinitialize_m =   IpplTimings::getTimer ("Initialize geometry");
    TisInside_m =     IpplTimings::getTimer ("Inside test");
    TfastIsInside_m = IpplTimings::getTimer ("Fast inside test");
    TRayTrace_m =     IpplTimings::getTimer ("Ray tracing");
    TPartInside_m =   IpplTimings::getTimer ("Particle Inside");
1058 1059 1060 1061

    h5FileName_m = Attributes::getString (itsAttr[FGEOM]);
    if (!h5FileName_m.empty ())
        initialize ();
1062
 }
gsell's avatar
gsell committed
1063 1064

BoundaryGeometry::~BoundaryGeometry() {
1065
    gsl_rng_free(randGen_m);
gsell's avatar
gsell committed
1066 1067
}

1068
bool BoundaryGeometry::canReplaceBy (Object* object) {
gsell's avatar
gsell committed
1069
    // Can replace only by another GEOMETRY.
1070
    return dynamic_cast<BGeometryBase*>(object) != 0;
gsell's avatar
gsell committed
1071 1072
}

1073
BoundaryGeometry* BoundaryGeometry::clone (const std::string& name) {
1074
    return new BoundaryGeometry (name, this);
gsell's avatar
gsell committed
1075 1076
}

1077 1078
void BoundaryGeometry::update () {
    if (getOpalName ().empty ()) setOpalName ("UNNAMED_GEOMETRY");
gsell's avatar
gsell committed
1079 1080 1081
}


1082 1083
void BoundaryGeometry::execute () {
    update ();
gsell's avatar
gsell committed
1084 1085 1086 1087 1088
    Tinitialize_m =   IpplTimings::getTimer ("Initialize geometry");
    TisInside_m =     IpplTimings::getTimer ("Inside test");
    TfastIsInside_m = IpplTimings::getTimer ("Fast inside test");
    TRayTrace_m =     IpplTimings::getTimer ("Ray tracing");
    TPartInside_m =   IpplTimings::getTimer ("Particle Inside");
gsell's avatar
gsell committed
1089 1090
}

1091
BoundaryGeometry* BoundaryGeometry::find (const std::string& name) {
1092 1093
    BoundaryGeometry* geom = dynamic_cast<BoundaryGeometry*>(
        OpalData::getInstance ()->find (name));
gsell's avatar
gsell committed
1094

1095
    if (geom == 0)
1096
        throw OpalException ("BoundaryGeometry::find()", "Geometry \""
1097 1098
                             + name + "\" not found.");
    return geom;
gsell's avatar
gsell committed
1099 1100
}

1101
void BoundaryGeometry::updateElement (ElementBase* /*element*/) {
gsell's avatar
gsell committed
1102 1103
}

1104 1105 1106
int
BoundaryGeometry::intersectTriangleVoxel (
    const int triangle_id,
1107 1108 1109
    const int i,
    const int j,
    const int k
1110
    ) {
gsell's avatar
gsell committed
1111 1112 1113 1114 1115
    const Triangle t(
        getPoint (triangle_id, 1),
        getPoint (triangle_id, 2),
        getPoint (triangle_id, 3)
        );
1116

1117
    const Vector_t P(
gsell's avatar
gsell committed
1118 1119 1120
        i * voxelMesh_m.sizeOfVoxel [0] + voxelMesh_m.minExtent[0],
        j * voxelMesh_m.sizeOfVoxel [1] + voxelMesh_m.minExtent[1],
        k * voxelMesh_m.sizeOfVoxel [2] + voxelMesh_m.minExtent[2]
1121
        );
1122

gsell's avatar
gsell committed
1123
    Voxel v(P, P+voxelMesh_m.sizeOfVoxel);
1124

1125
    return v.intersect (t);
1126 1127
}

gsell's avatar
gsell committed
1128 1129 1130 1131 1132
/*
  Find the 3D intersection of a line segment, ray or line with a triangle.

  Input:
        kind: type of test: SEGMENT, RAY or LINE
1133
        P0, P0: defining
gsell's avatar
gsell committed
1134 1135 1136 1137
            a line segment from P0 to P1 or
            a ray starting at P0 with directional vector P1-P0 or
            a line through P0 and P1
        V0, V1, V2: the triangle vertices
1138

gsell's avatar
gsell committed
1139 1140 1141
  Output:
        I: intersection point (when it exists)

1142 1143
  Return values for line segment and ray test :
        -1 = triangle is degenerated (a segment or point)
gsell's avatar
gsell committed
1144
        0 =  disjoint (no intersect)
1145 1146 1147 1148 1149 1150 1151 1152 1153 1154
        1 =  are in the same plane
        2 =  intersect in unique point I1

  Return values for line intersection test :
        -1: triangle is degenerated (a segment or point)
        0:  disjoint (no intersect)
        1:  are in the same plane
        2:  intersect in unique point I1, with r < 0.0
        3:  intersect in unique point I1, with 0.0 <= r <= 1.0
        4:  intersect in unique point I1, with 1.0 < r
gsell's avatar
gsell committed
1155 1156 1157 1158 1159 1160 1161 1162 1163 1164

  For algorithm and implementation see:
  http://geomalgorithms.com/a06-_intersect-2.html

  Copyright 2001 softSurfer, 2012 Dan Sunday
  This code may be freely used and modified for any purpose
  providing that this copyright notice is included with it.
  SoftSurfer makes no warranty for this code, and cannot be held
  liable for any real or imagined damage resulting from its use.
  Users of this code must verify correctness for their application.
gsell's avatar
gsell committed
1165
 */
gsell's avatar
gsell committed
1166 1167

int
1168
BoundaryGeometry::intersectLineTriangle (
1169
    const enum INTERSECTION_TESTS kind,
gsell's avatar
gsell committed
1170 1171 1172 1173
    const Vector_t& P0,
    const Vector_t& P1,
    const int triangle_id,
    Vector_t& I
1174
    ) {
1175 1176 1177
    const Vector_t V0 = getPoint (triangle_id, 1);
    const Vector_t V1 = getPoint (triangle_id, 2);
    const Vector_t V2 = getPoint (triangle_id, 3);
gsell's avatar
gsell committed
1178 1179 1180 1181 1182 1183 1184

    // get triangle edge vectors and plane normal
    const Vector_t u = V1 - V0;         // triangle vectors
    const Vector_t v = V2 - V0;
    const Vector_t n = cross (u, v);
    if (n == (Vector_t)0)               // triangle is degenerate
        return -1;                      // do not deal with this case
1185

gsell's avatar
gsell committed
1186 1187 1188 1189
    const Vector_t dir = P1 - P0;       // ray direction vector
    const Vector_t w0 = P0 - V0;
    const double a = -dot(n,w0);
    const double b = dot(n,dir);
1190 1191
    if (cmp::eq_zero(b)) {              // ray is  parallel to triangle plane
        if (cmp::eq_zero(a)) {          // ray lies in triangle plane
1192
            return 1;
1193
        } else {                        // ray disjoint from plane
1194
            return 0;
1195
        }
gsell's avatar
gsell committed
1196
    }
1197

gsell's avatar
gsell committed
1198 1199 1200 1201
    // get intersect point of ray with triangle plane