BoundaryGeometry.cpp 97.5 KB
Newer Older
gsell's avatar
gsell committed
1
/*
gsell's avatar
gsell committed
2
  Implementation of the class BoundaryGeometry.
gsell's avatar
gsell committed
3

gsell's avatar
gsell committed
4
  Copyright & License: See Copyright.readme in src directory
gsell's avatar
gsell committed
5 6
 */

7
#include <fstream>
gsell's avatar
gsell committed
8 9

#include "H5hut.h"
10

gsell's avatar
gsell committed
11 12 13 14
#include "Structure/BoundaryGeometry.h"
#include "Structure/PriEmissionPhysics.h"
#include "Expressions/SRefExpr.h"
#include "Elements/OpalBeamline.h"
gsell's avatar
gsell committed
15

gsell's avatar
gsell committed
16
extern Inform* gmsg;
17

gsell's avatar
gsell committed
18 19 20
#define SQR(x) ((x)*(x))
#define PointID(triangle_id, vertex_id) allbfaces_m[4 * (triangle_id) + (vertex_id)]
#define Point(triangle_id, vertex_id)   geo3Dcoords_m[allbfaces_m[4 * (triangle_id) + (vertex_id)]]
21 22 23

/*

gsell's avatar
gsell committed
24 25 26 27 28 29 30
  Some
   _   _      _                 
  | | | | ___| |_ __   ___ _ __ 
  | |_| |/ _ \ | '_ \ / _ \ '__|
  |  _  |  __/ | |_) |  __/ |   
  |_| |_|\___|_| .__/ \___|_|   
             |_|  
31

gsell's avatar
gsell committed
32
  functions
33
 */
gsell's avatar
gsell committed
34 35 36 37 38
struct VectorLessX {
    bool operator() (Vector_t x1, Vector_t x2) {
        return x1 (0) < x2 (0);
    }
};
39

gsell's avatar
gsell committed
40 41 42 43 44
struct VectorLessY {
    bool operator() (Vector_t x1, Vector_t x2) {
        return x1 (1) < x2 (1);
    }
};
45

gsell's avatar
gsell committed
46 47 48
struct VectorLessZ {
    bool operator() (Vector_t x1, Vector_t x2) {
        return x1 (2) < x2 (2);
49
    }
gsell's avatar
gsell committed
50 51 52 53 54 55 56 57 58 59 60 61 62
};

/**
   Calculate the maximum of coordinates of geometry,i.e the maximum of X,Y,Z
 */
Vector_t get_max_extend (std::vector<Vector_t>& coords) {
    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 ());
    return Vector_t (x (0), y (1), z (2));
63 64
}

gsell's avatar
gsell committed
65 66 67 68 69 70 71 72 73 74 75 76
/*
   Compute the minimum of coordinates of geometry, i.e the minimum of X,Y,Z
 */
Vector_t get_min_extend (std::vector<Vector_t>& coords) {
    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 ());
    return Vector_t (x (0), y (1), z (2));
}
77

gsell's avatar
gsell committed
78 79 80 81
/*
  "ULP compare" for double precision floating point numbers.
  See:
    http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
82

gsell's avatar
gsell committed
83 84 85
  Note:
    An updated version of this document with improved code is here:
    http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
86

gsell's avatar
gsell committed
87 88
 */
static int64_t fcmp (
89 90 91
    const double A,
    const double B,
    const int maxUlps ) {
gsell's avatar
gsell committed
92 93 94 95 96 97 98 99
                    
    // 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);
    assert (sizeof (long long) == sizeof (int64_t) );
    assert (sizeof (long long) == sizeof (double) );
                    
    // Make [ab]Int lexicographically ordered as a twos-complement int
100
    const double* pa = &A;
gsell's avatar
gsell committed
101 102 103 104
    int64_t aInt = *(int64_t*)pa;
    if (aInt < 0)
        aInt = 0x8000000000000000LL - aInt;
                    
105
    const double* pb = &B;
gsell's avatar
gsell committed
106 107 108 109
    int64_t bInt = *(int64_t*)pb;
    if (bInt < 0)
        bInt = 0x8000000000000000LL - bInt;
                    
110
    const int64_t intDiff = aInt - bInt;
gsell's avatar
gsell committed
111 112 113 114
    if (llabs(intDiff) <= maxUlps)
        return 0;
    return intDiff;
}
gsell's avatar
gsell committed
115

116 117 118 119 120 121 122 123 124 125 126 127 128 129
static inline bool is_in_voxel (Vector_t& point, Vector_t& min, Vector_t& max) {
    if (fcmp (point [0], min [0], 10) < 0) return false;
    if (fcmp (point [1], min [1], 10) < 0) return false;
    if (fcmp (point [2], min [2], 10) < 0) return false;
    if (fcmp (point [0], max [0], 10) > 0) return false;
    if (fcmp (point [1], max [1], 10) > 0) return false;
    if (fcmp (point [2], max [2], 10) > 0) return false;
    return true;
}

/*
  write legacy VTK file of voxel mesh
*/
static void write_voxel_mesh (
130
    const std::unordered_set<int> ids,
131 132 133 134 135 136 137 138 139 140 141 142
    const Vector_t hr_m,
    const Vektor<int,3> nr,
    const Vector_t origin
    ) {
    /*----------------------------------------------------------------------*/
    const size_t numpoints = 8 * ids.size ();
    std::ofstream of;
    of.open (string ("data/testBBox.vtk").c_str ());
    assert (of.is_open ());
    of.precision (6);
    
    of << "# vtk DataFile Version 2.0" << std::endl;
143
    of << "generated using BoundaryGeometry::computeMeshVoxelization"
144 145 146 147
       << std::endl;
    of << "ASCII" << std::endl << std::endl;
    of << "DATASET UNSTRUCTURED_GRID" << std::endl;
    of << "POINTS " << numpoints << " float" << std::endl;
148 149 150 151 152 153 154 155

    const auto end_it = ids.end();
    const auto nr0_times_nr1 = nr[0] * nr[1];
    for (auto id = ids.begin (); id != end_it; id++) {
        int k = (*id - 1) / nr0_times_nr1;
        int rest = (*id - 1) % nr0_times_nr1;
        int j = rest / nr[0];
        int i = rest % nr[0]; 
156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171

        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;
172 173 174
    const auto num_cells = ids.size ();
    of << "CELLS " << num_cells << " " << 9 * num_cells << std::endl;
    for (size_t i = 0; i < num_cells; i++)
175 176 177
        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;
178 179
    of << "CELL_TYPES " << num_cells << std::endl;
    for (size_t i = 0; i <  num_cells; i++)
180
        of << "11" << std::endl;
181
    of << "CELL_DATA " << num_cells << std::endl;
182 183
    of << "SCALARS " << "cell_attribute_data" << " float " << "1" << std::endl;
    of << "LOOKUP_TABLE " << "default" << std::endl;
184
    for (size_t i = 0; i <  num_cells; i++)
185 186 187
        of << (float)(i) << std::endl;
    of << std::endl;
    of << "COLOR_SCALARS " << "BBoxColor " << 4 << std::endl;
188
    for (size_t i = 0; i < num_cells; i++) {
189 190 191 192 193
        of << "1.0" << " 1.0 " << "0.0 " << "1.0" << std::endl;
    }
    of << std::endl;
}

194 195 196 197 198 199 200 201 202
/*___________________________________________________________________________

  Triangle-cube intersection test.

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

 */

203 204 205 206
#include <math.h>


#define LERP( A, B, C) ((B)+(A)*((C)-(B)))
207 208
#define MIN2(a,b) (((a) < (b)) ? (a) : (b))
#define MAX2(a,b) (((a) > (b)) ? (a) : (b))
209 210 211 212 213
#define MIN3(a,b,c) ((((a)<(b))&&((a)<(c))) ? (a) : (((b)<(c)) ? (b) : (c)))
#define MAX3(a,b,c) ((((a)>(b))&&((a)>(c))) ? (a) : (((b)>(c)) ? (b) : (c)))
#define INSIDE 0
#define OUTSIDE 1

214 215 216 217 218
typedef struct {
    Vector_t v1;
    Vector_t v2;
} LineSegment;

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 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 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 347 348 349 350 351
typedef struct {
    Vector_t v1;                 /* Vertex1 */
    Vector_t v2;                 /* Vertex2 */
    Vector_t v3;                 /* Vertex3 */
} Triangle; 

typedef struct {
    Vector_t v1;
    Vector_t v2;
} Cube;

/*___________________________________________________________________________*/

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

static inline int
face_plane (
    const Vector_t p
    ) {
    int outcode;

    outcode = 0;
    if (p[0] >  .5) outcode |= 0x01;
    if (p[0] < -.5) outcode |= 0x02;
    if (p[1] >  .5) outcode |= 0x04;
    if (p[1] < -.5) outcode |= 0x08;
    if (p[2] >  .5) outcode |= 0x10;
    if (p[2] < -.5) outcode |= 0x20;
    return(outcode);
}

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

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

static inline int
bevel_2d (
    const Vector_t p
    ) {
    int outcode;

    outcode = 0;
    if ( p[0] + p[1] > 1.0) outcode |= 0x001;
    if ( p[0] - p[1] > 1.0) outcode |= 0x002;
    if (-p[0] + p[1] > 1.0) outcode |= 0x004;
    if (-p[0] - p[1] > 1.0) outcode |= 0x008;
    if ( p[0] + p[2] > 1.0) outcode |= 0x010;
    if ( p[0] - p[2] > 1.0) outcode |= 0x020;
    if (-p[0] + p[2] > 1.0) outcode |= 0x040;
    if (-p[0] - p[2] > 1.0) outcode |= 0x080;
    if ( p[1] + p[2] > 1.0) outcode |= 0x100;
    if ( p[1] - p[2] > 1.0) outcode |= 0x200;
    if (-p[1] + p[2] > 1.0) outcode |= 0x400;
    if (-p[1] - p[2] > 1.0) outcode |= 0x800;
    return(outcode);
}

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

  Which of the eight corner plane(s) is point P outside of?
*/
static inline int
bevel_3d (
    const Vector_t p
    ) {
    int outcode;

    outcode = 0;
    if (( p[0] + p[1] + p[2]) > 1.5) outcode |= 0x01;
    if (( p[0] + p[1] - p[2]) > 1.5) outcode |= 0x02;
    if (( p[0] - p[1] + p[2]) > 1.5) outcode |= 0x04;
    if (( p[0] - p[1] - p[2]) > 1.5) outcode |= 0x08;
    if ((-p[0] + p[1] + p[2]) > 1.5) outcode |= 0x10;
    if ((-p[0] + p[1] - p[2]) > 1.5) outcode |= 0x20;
    if ((-p[0] - p[1] + p[2]) > 1.5) outcode |= 0x40;
    if ((-p[0] - p[1] - p[2]) > 1.5) outcode |= 0x80;
    return(outcode);
}

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

  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 (
    const Vector_t p1,
    const Vector_t p2,
    const double alpha,
    const int mask
    ) {
    Vector_t plane_point;

    plane_point[0] = LERP(alpha, p1[0], p2[0]);
    plane_point[1] = LERP(alpha, p1[1], p2[1]);
    plane_point[2] = LERP(alpha, p1[2], p2[2]);
    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 (
    const Vector_t p1,
    const Vector_t p2,
    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);
    if ((0x04 & outcode_diff) != 0) 
        if (check_point(p1,p2,( .5-p1[1])/(p2[1]-p1[1]),0x3b) == INSIDE) return(INSIDE);
    if ((0x08 & outcode_diff) != 0) 
        if (check_point(p1,p2,(-.5-p1[1])/(p2[1]-p1[1]),0x37) == INSIDE) return(INSIDE);
    if ((0x10 & outcode_diff) != 0) 
        if (check_point(p1,p2,( .5-p1[2])/(p2[2]-p1[2]),0x2f) == INSIDE) return(INSIDE);
    if ((0x20 & outcode_diff) != 0) 
        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
*/
352

353 354 355 356 357 358 359 360 361 362 363
#define EPS 10e-5

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

364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389
static int
point_triangle_intersection (
    const Vector_t p,
    const Triangle t
    ) {
    /*
      First, a quick bounding-box test:
      If P is outside triangle bbox, there cannot be an intersection.
    */
    if (p[0] > MAX3(t.v1[0], t.v2[0], t.v3[0])) return(OUTSIDE);  
    if (p[1] > MAX3(t.v1[1], t.v2[1], t.v3[1])) return(OUTSIDE);
    if (p[2] > MAX3(t.v1[2], t.v2[2], t.v3[2])) return(OUTSIDE);
    if (p[0] < MIN3(t.v1[0], t.v2[0], t.v3[0])) return(OUTSIDE);
    if (p[1] < MIN3(t.v1[1], t.v2[1], t.v3[1])) return(OUTSIDE);
    if (p[2] < MIN3(t.v1[2], t.v2[2], t.v3[2])) return(OUTSIDE);
    
    /*
      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
      to the outside of this triangle side.                                
    */
    const Vector_t vect12 = t.v1 - t.v2;
    const Vector_t vect1h = t.v1 - p;
    const Vector_t cross12_1p = cross (vect12, vect1h);
390
    const int sign12 = SIGN3(cross12_1p);      /* Extract X,Y,Z signs as 0..7 or 0...63 integer */
391 392 393 394

    const Vector_t vect23 = t.v2 - t.v3;
    const Vector_t vect2h = t.v2 - p;
    const Vector_t cross23_2p = cross (vect23, vect2h);
395
    const int sign23 = SIGN3(cross23_2p);
396 397 398 399

    const Vector_t vect31 = t.v3 - t.v1;
    const Vector_t vect3h = t.v3 - p;
    const Vector_t cross31_3p = cross (vect31, vect3h);
400
    const int sign31 = SIGN3(cross31_3p);
401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419

    /*
      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
420 421
triangle_intersects_cube (
    const Triangle& t
422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506
    ) {
    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!
    */
   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);

   /*
     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
   */
   v1_test |= bevel_2d(t.v1) << 8; 
   v2_test |= bevel_2d(t.v2) << 8; 
   v3_test |= bevel_2d(t.v3) << 8;
   if ((v1_test & v2_test & v3_test) != 0) return(OUTSIDE);  

   /*
     Now do the same trivial rejection test for the 8 corner planes
   */
   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);   

   /*
     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)
      if (check_line(t.v1,t.v2,v1_test|v2_test) == INSIDE) return(INSIDE);
   if ((v1_test & v3_test) == 0)
      if (check_line(t.v1,t.v3,v1_test|v3_test) == INSIDE) return(INSIDE);
   if ((v2_test & v3_test) == 0)
      if (check_line(t.v2,t.v3,v2_test|v3_test) == INSIDE) return(INSIDE);

   /*
     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.

     To find plane of the triangle, first perform crossproduct on 
     two triangle side vectors to compute the normal vector.
   */  
   Vector_t vect12 = t.v1 - t.v2;
   Vector_t vect13 = t.v1 - t.v3;
   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.
   */
   double d = norm[0] * t.v1[0] + norm[1] * t.v1[1] + norm[2] * t.v1[2];

   /*
     if one of the diagonals is parallel to the plane, the other will
     intersect the plane
   */
   double denom;
   if(fabs(denom=(norm[0] + norm[1] + norm[2]))>EPS) {
507
       /* skip parallel diagonals to the plane; division by 0 can occure */
508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553
       Vector_t hitpp = d / denom;
       if (fabs(hitpp[0]) <= 0.5)
           if (point_triangle_intersection(hitpp,t) == INSIDE) return(INSIDE);
   }
   if(fabs(denom=(norm[0] + norm[1] - norm[2]))>EPS) {
       Vector_t hitpn;
       hitpn[2] = -(hitpn[0] = hitpn[1] = d / denom);
       if (fabs(hitpn[0]) <= 0.5)
           if (point_triangle_intersection(hitpn,t) == INSIDE) return(INSIDE);
   }       
   if(fabs(denom=(norm[0] - norm[1] + norm[2]))>EPS) {       
       Vector_t hitnp;
       hitnp[1] = -(hitnp[0] = hitnp[2] = d / denom);
       if (fabs(hitnp[0]) <= 0.5)
           if (point_triangle_intersection(hitnp,t) == INSIDE) return(INSIDE);
   }
   if(fabs(denom=(norm[0] - norm[1] - norm[2]))>EPS) {
       Vector_t hitnn;
       hitnn[1] = hitnn[2] = -(hitnn[0] = d / denom);
       if (fabs(hitnn[0]) <= 0.5)
           if (point_triangle_intersection(hitnn,t) == INSIDE) return(INSIDE);
   }
   
   /*
     No edge touched the cube; no cube diagonal touched the triangle.
     We're done...there was no intersection.
   */
   return(OUTSIDE);
}

static inline void
scaleCube (
    Cube& c,
    const Vector_t& scale
    ) {
    c.v1[0] *= scale[0];
    c.v1[1] *= scale[1];
    c.v1[2] *= scale[2];
    c.v2[0] *= scale[0];
    c.v2[1] *= scale[1];
    c.v2[2] *= scale[2];
}

static inline void
scaleTriangle (
    Triangle& t,
554 555
    const Vector_t& scale,
    const Vector_t& shift
556 557 558 559 560 561 562 563 564 565 566 567 568 569 570
    ) {
    t.v1[0] *= scale[0];
    t.v1[1] *= scale[1];
    t.v1[2] *= scale[2];
    t.v2[0] *= scale[0];
    t.v2[1] *= scale[1];
    t.v2[2] *= scale[2];
    t.v3[0] *= scale[0];
    t.v3[1] *= scale[1];
    t.v3[2] *= scale[2];
    t.v1 -= shift;
    t.v2 -= shift;
    t.v3 -= shift;
}

571
static inline int
572 573 574 575 576 577 578 579 580
intersect3dTriangleCube (
    const Triangle& t,
    const Cube& c
    ) {
    Cube c_ = c;
    Triangle t_ = t;
    Vector_t extend = c_.v2 - c_.v1;
    const Vector_t scale = 1.0 / extend; 
    scaleCube (c_, scale);
581 582 583
    scaleTriangle (t_, scale , c_.v1 + 0.5);
    return triangle_intersects_cube (t_);
}
584

585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699

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

  Ray-cube intersection
 */
/*
 * 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
 * 
 */

class Ray {
public:
	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];
	}

	Vector_t origin;
	Vector_t direction;
	Vector_t inv_direction;
	int sign[3];
};

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

class Box {
public:
    Box() { }
    Box(const Vector_t &min, const Vector_t &max) {
        parameters[0] = min;
        parameters[1] = max;
    }
    // (t0, t1) is the interval for valid hits
    bool intersect (
        const Ray& r,
        double& tmin,
        double& tmax
        ) const {
	double tmin_ = (parameters[r.sign[0]][0]   - r.origin[0]) * r.inv_direction[0];
	double tmax_ = (parameters[1-r.sign[0]][0] - r.origin[0]) * r.inv_direction[0];
	double tymin = (parameters[r.sign[1]][1]   - r.origin[1]) * r.inv_direction[1];
	double tymax = (parameters[1-r.sign[1]][1] - r.origin[1]) * r.inv_direction[1];
	if ( (tmin_ > tymax) || (tymin > tmax_) ) 
		return false;
	if (tymin > tmin_)
		tmin_ = tymin;
	if (tymax < tmax_)
		tmax_ = tymax;
	double tzmin = (parameters[r.sign[2]][2]   - r.origin[2]) * r.inv_direction[2];
	double tzmax = (parameters[1-r.sign[2]][2] - r.origin[2]) * r.inv_direction[2];
	if ( (tmin_ > tzmax) || (tzmin > tmax_) ) 
		return false;
	if (tzmin > tmin_)
		tmin_ = tzmin;
	tmin = tmin_;
	if (tzmax < tmax_)
		tmax_ = tzmax;
	tmax = tmax_;
	return true;
    }

    // corners
    Vector_t parameters[2];
};

/*
  Result:
    0   no intersection
    1   ???
    2   intersection, origin of ray outside cube
    3   intersection, origin of ray inside cube
    4   intersection on extended line (in opposite direction)
 */
static inline int
intersect3dRayCube (
    const Ray& r,
    const Cube& c,
    double& tmin,
    double& tmax
    ) {
    Box b = Box(c.v1, c.v2);
    if (!b.intersect(r, tmin, tmax)) {
        return 0;
    }
    if (tmin > 0) {
        return 2;
    }
    if (tmin <= 0 && tmax > 0) {
        return 3;
    }
    return 4;
700 701
}

gsell's avatar
gsell committed
702 703 704 705 706 707 708 709
/*
  ____                           _              
 / ___| ___  ___  _ __ ___   ___| |_ _ __ _   _ 
| |  _ / _ \/ _ \| '_ ` _ \ / _ \ __| '__| | | |
| |_| |  __/ (_) | | | | | |  __/ |_| |  | |_| |
 \____|\___|\___/|_| |_| |_|\___|\__|_|   \__, |
                                          |___/
*/
710

gsell's avatar
gsell committed
711
BoundaryGeometry::BoundaryGeometry() :
712 713 714 715 716 717
    Definition (
        SIZE, "GEOMETRY", "The \"GEOMETRY\" statement defines the beam pipe geometry."),
    Tribarycent_m (NULL),
    TriPrPartloss_m (NULL),
    TriSePartloss_m (NULL),
    TriFEPartloss_m (NULL),
718 719
    allbfaces_m (NULL) {

gsell's avatar
gsell committed
720
    itsAttr[FGEOM] = Attributes::makeString
721 722 723
        ("FGEOM",
         "Specifies the geometry file [h5fed]",
         "");
gsell's avatar
gsell committed
724 725

    itsAttr[TOPO] = Attributes::makeString
726
        ("TOPO",
727
         "BOX, BOXCORNER, ELLIPTIC if FGEOM is selected topo is over-written ",
728
         "ELLIPTIC");
gsell's avatar
gsell committed
729

730 731
    itsAttr[LENGTH] = Attributes::makeReal
        ("LENGTH",
732 733
         "Specifies the length of a tube shaped elliptic beam pipe [m]",
         1.0);
gsell's avatar
gsell committed
734 735

    itsAttr[S] = Attributes::makeReal
736 737 738
        ("S",
         "Specifies the start of a tube shaped elliptic beam pipe [m]",
         0.0);
gsell's avatar
gsell committed
739 740

    itsAttr[A] = Attributes::makeReal
741 742 743
        ("A",
         "Specifies the major semi-axis of a tube shaped elliptic beam pipe [m]",
         0.025);
gsell's avatar
gsell committed
744 745

    itsAttr[B] = Attributes::makeReal
746 747 748
        ("B",
         "Specifies the major semi-axis of a tube shaped elliptic beam pipe [m]",
         0.025);
gsell's avatar
gsell committed
749 750

    itsAttr[L1] = Attributes::makeReal
751 752 753
        ("L1",
         "In case of BOXCORNER Specifies first part with hight == B [m]",
         0.5);
gsell's avatar
gsell committed
754 755

    itsAttr[L2] = Attributes::makeReal
756 757 758
        ("L2",
         "In case of BOXCORNER Specifies first second with hight == B-C [m]",
         0.2);
gsell's avatar
gsell committed
759 760

    itsAttr[C] = Attributes::makeReal
761 762 763
        ("C",
         "In case of BOXCORNER Specifies hight of corner C [m]",
         0.01);
gsell's avatar
gsell committed
764 765

    itsAttr[DISTR] = Attributes::makeString
766 767 768
        ("DISTR",
         "Distribution to be generated on the surface",
         "");
gsell's avatar
gsell committed
769
    itsAttr[DISTRS] = Attributes::makeStringArray
770 771
        ("DISTRS",
         "Distribution array to be generated on the surface");
gsell's avatar
gsell committed
772 773

    itsAttr[XYZSCALE] = Attributes::makeReal
774 775 776
        ("XYZSCALE",
         "Multiplicative scaling factor for coordinates ",
         1.0);
gsell's avatar
gsell committed
777

778 779 780 781 782 783 784 785 786 787 788 789 790 791 792
    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
793
    itsAttr[ZSHIFT] = Attributes::makeReal
794 795 796
        ("ZSHIFT",
         "Shift in z direction",
         0.0);
gsell's avatar
gsell committed
797

adelmann's avatar
adelmann committed
798 799 800
    itsAttr[APERTURE]  = Attributes::makeRealArray
        ("APERTURE", "The element aperture");

801
    BoundaryGeometry* defGeometry = clone ("UNNAMED_GEOMETRY");
gsell's avatar
gsell committed
802 803
    defGeometry->builtin = true;

804 805 806 807 808
    TPInside_m = IpplTimings::getTimer ("Particle Inside");
    TPreProc_m = IpplTimings::getTimer ("Pre Processing");
    TRayTrace_m = IpplTimings::getTimer ("Ray tracing");
    h5FileName_m = Attributes::getString (itsAttr[FGEOM]);
    Tinward_m = IpplTimings::getTimer ("Check inward");
gsell's avatar
gsell committed
809
    try {
810 811 812
        defGeometry->update ();
        OpalData::getInstance ()->define (defGeometry);
    } catch (...) {
gsell's avatar
gsell committed
813 814
        delete defGeometry;
    }
815 816 817
    gsl_rng_env_setup();
    randGen_m = gsl_rng_alloc(gsl_rng_default);    

818 819
    if (!h5FileName_m.empty ())
        initialize ();
820

gsell's avatar
gsell committed
821 822 823
}

BoundaryGeometry::BoundaryGeometry(
824 825 826 827 828 829 830 831
    const string& name,
    BoundaryGeometry* parent
    ) :
    Definition (name, parent),
    Tribarycent_m (NULL),
    TriPrPartloss_m (NULL),
    TriSePartloss_m (NULL),
    TriFEPartloss_m (NULL),
832
    allbfaces_m (NULL) {
833 834 835
    gsl_rng_env_setup();
    randGen_m = gsl_rng_alloc(gsl_rng_default);    

836 837 838 839 840 841 842
    h5FileName_m = Attributes::getString (itsAttr[FGEOM]);
    if (!h5FileName_m.empty ())
        initialize ();
    TPInside_m = IpplTimings::getTimer ("Particle Inside");
    TPreProc_m = IpplTimings::getTimer ("zshiftPre Processing");
    TRayTrace_m = IpplTimings::getTimer ("Ray tracing");
    Tinward_m = IpplTimings::getTimer ("Check inward");
gsell's avatar
gsell committed
843

844
 }
gsell's avatar
gsell committed
845 846 847

BoundaryGeometry::~BoundaryGeometry() {
    /*
848
       if (allbfaces_m)
gsell's avatar
gsell committed
849
          delete allbfaces_m;
850
       if (Tribarycent_m)
gsell's avatar
gsell committed
851
          delete Tribarycent_m;
852 853 854 855 856 857 858
       if (TriPrPartloss_m )
       delete TriPrPartloss_m ;
       if (TriFEPartloss_m)
       delete TriFEPartloss_m;
       if (TriSePartloss_m)
       delete TriSePartloss_m;
     */
859 860

    gsl_rng_free(randGen_m);
gsell's avatar
gsell committed
861 862
}

863
bool BoundaryGeometry::canReplaceBy (Object* object) {
gsell's avatar
gsell committed
864
    // Can replace only by another GEOMETRY.
865
    return dynamic_cast<BGeometryBase*>(object) != 0;
gsell's avatar
gsell committed
866 867
}

868 869
BoundaryGeometry* BoundaryGeometry::clone (const string& name) {
    return new BoundaryGeometry (name, this);
gsell's avatar
gsell committed
870 871
}

872 873
void BoundaryGeometry::update () {
    if (getOpalName ().empty ()) setOpalName ("UNNAMED_GEOMETRY");
gsell's avatar
gsell committed
874 875 876
}


877 878 879 880 881 882
void BoundaryGeometry::execute () {
    update ();
    TPInside_m = IpplTimings::getTimer ("Particle Inside");
    TPreProc_m = IpplTimings::getTimer ("Pre Processing");
    TRayTrace_m = IpplTimings::getTimer ("Ray tracing");
    Tinward_m = IpplTimings::getTimer ("Check inward");
gsell's avatar
gsell committed
883 884
}

885 886 887
BoundaryGeometry* BoundaryGeometry::find (const string& name) {
    BoundaryGeometry* geom = dynamic_cast<BoundaryGeometry*>(
        OpalData::getInstance ()->find (name));
gsell's avatar
gsell committed
888

889
    if (geom == 0)
890
        throw OpalException ("BoundaryGeometry::find()", "Geometry \""
891 892
                             + name + "\" not found.");
    return geom;
gsell's avatar
gsell committed
893 894
}

gsell's avatar
gsell committed
895
void BoundaryGeometry::updateElement (ElementBase* element) {
gsell's avatar
gsell committed
896 897
}

898 899 900
int
BoundaryGeometry::intersectTriangleVoxel (
    const int triangle_id,
901 902 903
    const int i,
    const int j,
    const int k
904 905 906 907 908
    ) {
    const Triangle t = {getPoint (triangle_id, 1),
                        getPoint (triangle_id, 2),
                        getPoint (triangle_id, 3)};

909 910 911 912
    //const int k = (voxel_id - 1) / (nr_m[0] * nr_m[1]);
    //const int rest = (voxel_id - 1) % (nr_m[0] * nr_m[1]);
    //const int j = rest / nr_m[0];
    //const int i = rest % nr_m[0]; 
913 914 915 916 917 918 919 920
    
    Cube c;
    c.v1 = {
        i * hr_m[0] + voxelMesh_m.minExtend[0],
        j * hr_m[1] + voxelMesh_m.minExtend[1],
        k * hr_m[2] + voxelMesh_m.minExtend[2]
    };

921
    c.v2 = c.v1 + hr_m;
922 923 924 925

    return intersect3dTriangleCube (t, c);
}

gsell's avatar
gsell committed
926 927 928 929 930 931 932 933 934 935 936 937 938 939
/*
  Find the 3D intersection of a line segment, ray or line with a triangle.

  Input:
        kind: type of test: SEGMENT, RAY or LINE
        P0, P0: defining 
            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
        
  Output:
        I: intersection point (when it exists)

940 941
  Return values for line segment and ray test :
        -1 = triangle is degenerated (a segment or point)
gsell's avatar
gsell committed
942
        0 =  disjoint (no intersect)
943 944 945 946 947 948 949 950 951 952
        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
953 954 955 956 957 958 959 960 961 962

  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
963
 */
gsell's avatar
gsell committed
964 965

int
966
BoundaryGeometry::intersectLineTriangle (
967
    const enum INTERSECTION_TESTS kind,
gsell's avatar
gsell committed
968 969 970 971
    const Vector_t& P0,
    const Vector_t& P1,
    const int triangle_id,
    Vector_t& I
972
    ) {
973 974 975
    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
976 977 978 979 980 981 982 983 984 985 986 987

    // 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
    
    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);
988
    if (fcmp (b, 0.0, 10) == 0) {       // ray is  parallel to triangle plane
989
        if (a == 0) {                   // ray lies in triangle plane
990
            return 1;
991
        } else {                        // ray disjoint from plane
992
            return 0;
993
        }
gsell's avatar
gsell committed
994
    }
gsell's avatar
gsell committed
995 996 997 998 999
    
    // get intersect point of ray with triangle plane
    const double r = a / b;
    switch (kind) {
    case RAY:
1000
        if (r < 0.0) {                  // ray goes away from triangle
gsell's avatar
gsell committed
1001
            return 0;                   // => no intersect
1002
        }
gsell's avatar
gsell committed
1003
    case SEGMENT:
1004
        if (r < 0 || 1.0 < r) {         // intersection on extended
gsell's avatar
gsell committed
1005
            return 0;                   // segment
1006
        }
gsell's avatar
gsell committed
1007 1008 1009
    case LINE:
        break;
    };
1010
    I = P0 + r * dir;                   // intersect point of ray and plane
gsell's avatar
gsell committed
1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022
    
    // is I inside T?
    const double uu = dot(u,u);
    const double uv = dot(u,v);
    const double vv = dot(v,v);
    const Vector_t w = I - V0;
    const double wu = dot(w,u);
    const double wv = dot(w,v);
    const double D = uv * uv - uu * vv;
    
    // get and test parametric coords
    const double s = (uv * wv - vv * wu) / D;
1023
    if (s < 0.0 || s > 1.0) {           // I is outside T
gsell's avatar
gsell committed
1024
        return 0;
1025
    }
gsell's avatar
gsell committed
1026
    const double t = (uv * wu - uu * wv) / D;
1027
    if (t < 0.0 || (s + t) > 1.0) {     // I is outside T
gsell's avatar
gsell committed
1028
        return 0;
1029
    }
1030
    // intersection point is in triangle
1031 1032 1033
    if (r < 0.0) {                      // in extended segment in opposite
        return 2;                       // direction of ray
    } else if ((0.0 <= r) && (r <= 1.0)) { // in segment
1034
        return 3;
1035 1036
    } else {                            // in extended segment in
        return 4;                       // direction of ray 
1037
    }
gsell's avatar
gsell committed
1038 1039
}

1040 1041 1042 1043 1044 1045
static inline double magnitude (
    const Vector_t& v
    ) {
    return sqrt (dot (v,v));
}

1046 1047 1048 1049 1050 1051 1052 1053 1054
/*
  Game plan:
  Count number of intersection of the line segment defined by P and a reference
  pt with the boundary. If the reference pt is inside the boundary and the number
  of intersections is even, then P is inside the geometry. Otherwise P is outside.
  To count the number of intersection, we divide the line segment in N segments
  and run the line-segment boundary intersection test for all these segments.
  N must be choosen carefully. It shouldn't be to large to avoid needless test.
 */
1055 1056 1057 1058
int
BoundaryGeometry::fastIsInside (
    const Vector_t reference_pt,        // [in] reference pt inside the boundary
    const Vector_t P                    // [in] pt to test
1059
    ) {
1060
    const Vector_t v = reference_pt - P; // Bugfix: Switched the direction here -DW
1061 1062 1063 1064 1065 1066 1067 1068
    const int N = ceil (magnitude (v) / MIN3 (hr_m[0], hr_m[1], hr_m[2]));
    const Vector_t v_ = v / N;
    Vector_t P0 = P;
    Vector_t P1 = P + v_;
    Vector_t I;
    int triangle_id = -1;
    int result = 0;
    for (int i = 0; i < N; i++) {
1069 1070
        //result += intersectLineSegmentBoundary4PartInside (P0, P1, I, triangle_id) == 3 ? 1 : 0;
        result += intersectLineSegmentBoundary (P0, P1, I, triangle_id) == 3 ? 1 : 0;
1071 1072 1073 1074 1075 1076
        P0 = P1;
        P1 += v_;
    }
    return result;
}

1077 1078 1079
/*
  P must be *inside* the boundary geometry!
 */
1080 1081
int
BoundaryGeometry::intersectRayBoundary (
1082 1083
    const Vector_t& P,
    const Vector_t& v,
1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094
    Vector_t& I
    ) {
    /*
      set P1 to intersection of ray with bbox of voxel mesh 
      run line segment boundary intersection test with P and P1
     */
    Ray r = Ray (P, v);
    Box c = Box (mincoords_m, maxcoords_m);
    double tmin = 0.0;
    double tmax = 0.0;
    c.intersect (r, tmin, tmax);
1095
    int triangle_id = -1;
1096
    return (3 == intersectLineSegmentBoundary (P, P + tmax*v, I, triangle_id)) ? 1 : 0;
1097
}
1098

gsell's avatar
gsell committed
1099 1100 1101 1102 1103 1104 1105
/*
  Map point to unique voxel ID.

  Remember:
  * hr_m:  is the  mesh size
  * nr_m:  number of mesh points
  */
1106
inline int
1107
BoundaryGeometry::mapVoxelIndices2ID (
1108 1109 1110 1111 1112 1113 1114
    const int i,
    const int j,
    const int k
    ) {
    if (i < 0 || i >= nr_m[0] ||
        j < 0 || j >= nr_m[1] ||
        k < 0 || k >= nr_m[2]) {
gsell's avatar
gsell committed
1115
        return 0;
gsell's avatar
gsell committed
1116
    }
1117 1118
    return 1 + k * nr_m[0] * nr_m[1] + j * nr_m[0] + i;
}
1119

1120
inline bool
1121 1122 1123 1124 1125 1126 1127 1128 1129
BoundaryGeometry::mapPoint2VoxelIndices(
    const Vector_t pt,
    int& i,
    int& j,
    int& k
    ) {
    i = floor ((pt[0] - voxelMesh_m.minExtend [0]) / hr_m[0]);
    j = floor ((pt[1] - voxelMesh_m.minExtend [1]) / hr_m[1]);
    k = floor ((pt[2] - voxelMesh_m.minExtend [2]) / hr_m[2]);
1130 1131 1132 1133 1134 1135
    if (0 <= i && i < nr_m[0] &&
        0 <= j && j < nr_m[1] &&
        0 <= k && k < nr_m[2]) {
        return true;
    }
    return false;
1136
}
1137 1138 1139 1140 1141
    
inline int
BoundaryGeometry::mapPoint2VoxelID (
    const Vector_t x
    ) {
1142
    int i, j, k;
1143 1144 1145 1146
    if (mapPoint2VoxelIndices (x, i, j, k)) {
        return mapVoxelIndices2ID (i, j, k);
    }
    return -1;
1147 1148
}

1149
inline Vector_t
1150 1151 1152 1153
BoundaryGeometry::mapIndices2Voxel (
    const int i,
    const int j,
    const int k
1154
    ) {
1155
    return Vector_t (
1156 1157
        i * hr_m[0] + voxelMesh_m.minExtend[0],
        j * hr_m[1] + voxelMesh_m.minExtend[1],
1158
        k * hr_m[2] + voxelMesh_m.minExtend[2]);
gsell's avatar
gsell committed
1159 1160
}

1161
inline Vector_t
1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172
BoundaryGeometry::mapPoint2Voxel (
    const Vector_t& pt
    ) {
    const int i = floor ((pt[0] - voxelMesh_m.minExtend [0]) / hr_m[0]);
    const int j = floor ((pt[1] - voxelMesh_m.minExtend [1]) / hr_m[1]);
    const int k = floor ((pt[2] - voxelMesh_m.minExtend [2]) / hr_m[2]);

    return mapIndices2Voxel (i, j, k);
}


gsell's avatar
gsell committed
1173
void BoundaryGeometry::initialize () {
gsell's avatar
gsell committed
1174

gsell's avatar
gsell committed
1175
    class Local {
gsell's avatar
gsell committed
1176

gsell's avatar
gsell committed
1177
    public:
gsell's avatar
gsell committed
1178

gsell's avatar
gsell committed
1179
        static void computeGeometryInterval (BoundaryGeometry* bg) {
gsell's avatar
gsell committed
1180

gsell's avatar
gsell committed
1181 1182 1183
            bg->mincoords_m = get_min_extend (bg->geo3Dcoords_m);
            bg->maxcoords_m = get_max_extend (bg->geo3Dcoords_m);
            bg->len_m = bg->maxcoords_m - bg->mincoords_m;
gsell's avatar
gsell committed
1184

gsell's avatar
gsell committed
1185 1186 1187 1188
            /*
              Calculate the maximum dimension of triangles. This value will be used to
              define the cubic box size
            */
gsell's avatar
gsell committed
1189 1190 1191
            bg->longest_side_max_m = 0.0;
            bg->longest_side_min_m = 0.01;
            for (int i = 0; i < bg->num_triangles_m; i++) {
1192
                // compute length of longest edge
1193 1194 1195 1196
                const Vector_t x1 = bg->getPoint (i, 1);
                const Vector_t x2 = bg->getPoint (i, 2);
                const Vector_t x3 = bg->getPoint (i, 3);
                const double length_edge1 = sqrt (
1197
                    SQR (x1[0] - x2[0]) + SQR (x1[1] - x2[1]) + SQR (x1[2] - x2[2]));
1198
                const double length_edge2 = sqrt (
1199
                    SQR (x3[0] - x2[0]) + SQR (x3[1] - x2[1]) + SQR (x3[2] - x2[2]));
1200
                const double length_edge3 = sqrt (
1201
                    SQR (x3[0] - x1[0]) + SQR (x3[1] - x1[1]) + SQR (x3[2] - x1[2]));
1202
        
1203 1204 1205
                double max = length_edge1;
                if (length_edge2 > max) max = length_edge2;
                if (length_edge3 > max) max = length_edge3;
1206
        
1207
                // save min and max of length of longest edge
gsell's avatar
gsell committed
1208 1209
                if (bg->longest_side_max_m < max) bg->longest_side_max_m = max;
                if (bg->longest_side_min_m > max) bg->longest_side_min_m = max;