1 /* libctl: flexible Guile-based control files for scientific software
2 * Copyright (C) 1998-2020 Massachusetts Institute of Technology and Steven G. Johnson
3 *
4 * This library is free software; you can redistribute it and/or
5 * modify it under the terms of the GNU Lesser General Public
6 * License as published by the Free Software Foundation; either
7 * version 2 of the License, or (at your option) any later version.
8 *
9 * This library is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 * Lesser General Public License for more details.
13 *
14 * You should have received a copy of the GNU Lesser General Public
15 * License along with this library; if not, write to the
16 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
17 * Boston, MA 02111-1307, USA.
18 *
19 * Steven G. Johnson can be contacted at stevenj@alum.mit.edu.
20 */
21
22 #define _GNU_SOURCE
23 #include <stdlib.h>
24 #include <stdio.h>
25 #include <string.h>
26 #include <math.h>
27 #include <stdarg.h>
28
29 #ifndef LIBCTLGEOM
30 #include "ctl-io.h"
31 #else
32 #define material_type void *
material_type_copy(void ** src,void ** dest)33 static void material_type_copy(void **src, void **dest) { *dest = *src; }
34 #endif
35 #include "ctlgeom.h"
36
37 #ifdef CXX_CTL_IO
38 using namespace ctlio;
39 #define CTLIO ctlio::
40 #define GEOM geometric_object::
41 #define BLK block::
42 #define CYL cylinder::
43 #define MAT material_type::
44 #else
45 #define CTLIO
46 #define GEOM
47 #define BLK
48 #define CYL
49 #define MAT
50 #endif
51
52 #ifdef __cplusplus
53 #define MALLOC(type, num) (new type[num])
54 #define MALLOC1(type) (new type)
55 #define FREE(p) delete[](p)
56 #define FREE1(p) delete (p)
57 #else
58 #define MALLOC(type, num) ((type *)malloc(sizeof(type) * (num)))
59 #define MALLOC1(type) MALLOC(type, 1)
60 #define FREE(p) free(p)
61 #define FREE1(p) free(p)
62 #endif
63
64 #define K_PI 3.14159265358979323846
65 #define CHECK(cond, s) \
66 if (!(cond)) { \
67 fprintf(stderr, s "\n"); \
68 exit(EXIT_FAILURE); \
69 }
70
71 #define MAX(a, b) ((a) > (b) ? (a) : (b))
72 #define MIN(a, b) ((a) < (b) ? (a) : (b))
73
74 // forward declarations of prism-related routines, at the bottom of this file
75 static boolean node_in_polygon(double qx, double qy, vector3 *nodes, int num_nodes);
76 static boolean point_in_prism(prism *prsm, vector3 pc);
77 static vector3 normal_to_prism(prism *prsm, vector3 pc);
78 static double intersect_line_segment_with_prism(prism *prsm, vector3 pc, vector3 dc, double a,
79 double b);
80 static double get_prism_volume(prism *prsm);
81 static void get_prism_bounding_box(prism *prsm, geom_box *box);
82 static void display_prism_info(int indentby, geometric_object *o);
83 static void init_prism(geometric_object *o);
84 /**************************************************************************/
85
86 /* Allows writing to Python's stdout when running from Meep's Python interface */
87 void (*ctl_printf_callback)(const char *s) = NULL;
88
ctl_printf(const char * fmt,...)89 void ctl_printf(const char *fmt, ...) {
90 va_list ap;
91 va_start(ap, fmt);
92 if (ctl_printf_callback) {
93 char *s;
94 CHECK(vasprintf(&s, fmt, ap) >= 0, "vasprintf failed");
95 ctl_printf_callback(s);
96 free(s);
97 }
98 else {
99 vprintf(fmt, ap);
100 fflush(stdout);
101 }
102 va_end(ap);
103 }
104
105 /* If v is a vector in the lattice basis, normalize v so that
106 its cartesian length is unity. */
lattice_normalize(vector3 * v)107 static void lattice_normalize(vector3 *v) {
108 *v = vector3_scale(
109 1.0 / sqrt(vector3_dot(*v, matrix3x3_vector3_mult(geometry_lattice.metric, *v))), *v);
110 }
111
lattice_to_cartesian(vector3 v)112 static vector3 lattice_to_cartesian(vector3 v) {
113 return matrix3x3_vector3_mult(geometry_lattice.basis, v);
114 }
115
cartesian_to_lattice(vector3 v)116 static vector3 cartesian_to_lattice(vector3 v) {
117 return matrix3x3_vector3_mult(matrix3x3_inverse(geometry_lattice.basis), v);
118 }
119
120 /* geom_fix_object_ptr is called after an object's externally-configurable parameters
121 have been initialized, but before any actual geometry calculations are done;
122 it is an opportunity to (re)compute internal data fields (such as cached
123 rotation matrices) that depend on externally-configurable parameters.
124
125 One example: "Fix" the parameters of the given object to account for the
126 geometry_lattice basis, which may be non-orthogonal. In particular,
127 this means that the normalization of several unit vectors, such
128 as the cylinder or block axes, needs to be changed.
129
130 Unfortunately, we can't do this stuff at object-creation time
131 in Guile, because the geometry_lattice variable may not have
132 been assigned to its final value at that point. */
geom_fix_object_ptr(geometric_object * o)133 void geom_fix_object_ptr(geometric_object *o) {
134 switch (o->which_subclass) {
135 case GEOM CYLINDER:
136 lattice_normalize(&o->subclass.cylinder_data->axis);
137 if (o->subclass.cylinder_data->which_subclass == CYL WEDGE) {
138 vector3 a = o->subclass.cylinder_data->axis;
139 vector3 s = o->subclass.cylinder_data->subclass.wedge_data->wedge_start;
140 double p = vector3_dot(s, matrix3x3_vector3_mult(geometry_lattice.metric, a));
141 o->subclass.cylinder_data->subclass.wedge_data->e1 = vector3_minus(s, vector3_scale(p, a));
142 lattice_normalize(&o->subclass.cylinder_data->subclass.wedge_data->e1);
143 o->subclass.cylinder_data->subclass.wedge_data->e2 = cartesian_to_lattice(vector3_cross(
144 lattice_to_cartesian(o->subclass.cylinder_data->axis),
145 lattice_to_cartesian(o->subclass.cylinder_data->subclass.wedge_data->e1)));
146 }
147 break;
148 case GEOM BLOCK: {
149 matrix3x3 m;
150 lattice_normalize(&o->subclass.block_data->e1);
151 lattice_normalize(&o->subclass.block_data->e2);
152 lattice_normalize(&o->subclass.block_data->e3);
153 m.c0 = o->subclass.block_data->e1;
154 m.c1 = o->subclass.block_data->e2;
155 m.c2 = o->subclass.block_data->e3;
156 o->subclass.block_data->projection_matrix = matrix3x3_inverse(m);
157 break;
158 }
159 case GEOM PRISM: {
160 init_prism(o);
161 break;
162 }
163 case GEOM COMPOUND_GEOMETRIC_OBJECT: {
164 int i;
165 int n = o->subclass.compound_geometric_object_data->component_objects.num_items;
166 geometric_object *os = o->subclass.compound_geometric_object_data->component_objects.items;
167 for (i = 0; i < n; ++i) {
168 #if MATERIAL_TYPE_ABSTRACT
169 if (os[i].material.which_subclass == MAT MATERIAL_TYPE_SELF)
170 material_type_copy(&o->material, &os[i].material);
171 #endif
172 geom_fix_object_ptr(os + i);
173 }
174 break;
175 }
176 case GEOM GEOMETRIC_OBJECT_SELF:
177 case GEOM SPHERE: break; /* these objects are fine */
178 }
179 }
180
181 // deprecated API — doesn't work for prisms
geom_fix_object(geometric_object o)182 void geom_fix_object(geometric_object o) { geom_fix_object_ptr(&o); }
183
184 /* fix all objects in the geometry list as described in
185 geom_fix_object, above */
geom_fix_object_list(geometric_object_list geometry)186 void geom_fix_object_list(geometric_object_list geometry) {
187 int index;
188
189 for (index = 0; index < geometry.num_items; ++index)
190 geom_fix_object_ptr(geometry.items + index);
191 }
192
geom_fix_objects0(geometric_object_list geometry)193 void geom_fix_objects0(geometric_object_list geometry) { geom_fix_object_list(geometry); }
194
geom_fix_objects(void)195 void geom_fix_objects(void) { geom_fix_object_list(geometry); }
196
geom_fix_lattice0(lattice * L)197 void geom_fix_lattice0(lattice *L) {
198 L->basis1 = unit_vector3(L->basis1);
199 L->basis2 = unit_vector3(L->basis2);
200 L->basis3 = unit_vector3(L->basis3);
201 L->b1 = vector3_scale(L->basis_size.x, L->basis1);
202 L->b2 = vector3_scale(L->basis_size.y, L->basis2);
203 L->b3 = vector3_scale(L->basis_size.z, L->basis3);
204 L->basis.c0 = L->b1;
205 L->basis.c1 = L->b2;
206 L->basis.c2 = L->b3;
207 L->metric = matrix3x3_mult(matrix3x3_transpose(L->basis), L->basis);
208 }
209
geom_fix_lattice(void)210 void geom_fix_lattice(void) { geom_fix_lattice0(&geometry_lattice); }
211
geom_cartesian_lattice0(lattice * L)212 void geom_cartesian_lattice0(lattice *L) {
213 L->basis1.x = 1;
214 L->basis1.y = 0;
215 L->basis1.z = 0;
216 L->basis2.x = 0;
217 L->basis2.y = 1;
218 L->basis2.z = 0;
219 L->basis3.x = 0;
220 L->basis3.y = 0;
221 L->basis3.z = 1;
222 L->basis_size.x = L->basis_size.y = L->basis_size.z = 1;
223 geom_fix_lattice0(L);
224 }
225
geom_cartesian_lattice(void)226 void geom_cartesian_lattice(void) { geom_cartesian_lattice0(&geometry_lattice); }
227
geom_initialize(void)228 void geom_initialize(void) {
229 /* initialize many of the input variables that are normally
230 initialized from Scheme, except for default_material and
231 geometry_lattice.size. */
232 geom_cartesian_lattice();
233 geometry_center.x = geometry_center.y = geometry_center.z = 0;
234 dimensions = 3;
235 ensure_periodicity = 1;
236 geometry.num_items = 0;
237 geometry.items = 0;
238 }
239
240 /**************************************************************************/
241
242 /* Return whether or not the point p (in the lattice basis) is inside
243 the object o.
244
245 Requires that the global input var geometry_lattice already be
246 initialized.
247
248 point_in_fixed_objectp additionally requires that geom_fix_object
249 has been called on o (if the lattice basis is non-orthogonal). */
250
point_in_objectp(vector3 p,geometric_object o)251 boolean CTLIO point_in_objectp(vector3 p, geometric_object o) {
252 geom_fix_object_ptr(&o);
253 return point_in_fixed_objectp(p, o);
254 }
255
point_in_fixed_objectp(vector3 p,geometric_object o)256 boolean point_in_fixed_objectp(vector3 p, geometric_object o) {
257 return point_in_fixed_pobjectp(p, &o);
258 }
259
260 /* as point_in_fixed_objectp, but sets o to the object in question (if true)
261 (which may be different from the input o if o is a compound object) */
point_in_fixed_pobjectp(vector3 p,geometric_object * o)262 boolean point_in_fixed_pobjectp(vector3 p, geometric_object *o) {
263 vector3 r = vector3_minus(p, o->center);
264
265 switch (o->which_subclass) {
266 case GEOM GEOMETRIC_OBJECT_SELF: return 0;
267 case GEOM SPHERE: {
268 number radius = o->subclass.sphere_data->radius;
269 return (radius > 0.0 && vector3_dot(r, matrix3x3_vector3_mult(geometry_lattice.metric, r)) <=
270 radius * radius);
271 }
272 case GEOM CYLINDER: {
273 vector3 rm = matrix3x3_vector3_mult(geometry_lattice.metric, r);
274 number proj = vector3_dot(o->subclass.cylinder_data->axis, rm);
275 number height = o->subclass.cylinder_data->height;
276 if (fabs(proj) <= 0.5 * height) {
277 number radius = o->subclass.cylinder_data->radius;
278 if (o->subclass.cylinder_data->which_subclass == CYL CONE)
279 radius += (proj / height + 0.5) *
280 (o->subclass.cylinder_data->subclass.cone_data->radius2 - radius);
281 else if (o->subclass.cylinder_data->which_subclass == CYL WEDGE) {
282 number x = vector3_dot(rm, o->subclass.cylinder_data->subclass.wedge_data->e1);
283 number y = vector3_dot(rm, o->subclass.cylinder_data->subclass.wedge_data->e2);
284 number theta = atan2(y, x);
285 number wedge_angle = o->subclass.cylinder_data->subclass.wedge_data->wedge_angle;
286 if (wedge_angle > 0) {
287 if (theta < 0) theta = theta + 2 * K_PI;
288 if (theta > wedge_angle) return 0;
289 }
290 else {
291 if (theta > 0) theta = theta - 2 * K_PI;
292 if (theta < wedge_angle) return 0;
293 }
294 }
295 return (radius != 0.0 && vector3_dot(r, rm) - proj * proj <= radius * radius);
296 }
297 else
298 return 0;
299 }
300 case GEOM BLOCK: {
301 vector3 proj = matrix3x3_vector3_mult(o->subclass.block_data->projection_matrix, r);
302 switch (o->subclass.block_data->which_subclass) {
303 case BLK BLOCK_SELF: {
304 vector3 size = o->subclass.block_data->size;
305 return (fabs(proj.x) <= 0.5 * size.x && fabs(proj.y) <= 0.5 * size.y &&
306 fabs(proj.z) <= 0.5 * size.z);
307 }
308 case BLK ELLIPSOID: {
309 vector3 isa = o->subclass.block_data->subclass.ellipsoid_data->inverse_semi_axes;
310 double a = proj.x * isa.x, b = proj.y * isa.y, c = proj.z * isa.z;
311 return (a * a + b * b + c * c <= 1.0);
312 }
313 }
314 break; // never get here but silence compiler warning
315 }
316 case GEOM PRISM: {
317 return point_in_prism(o->subclass.prism_data, p);
318 }
319 case GEOM COMPOUND_GEOMETRIC_OBJECT: {
320 int i;
321 int n = o->subclass.compound_geometric_object_data->component_objects.num_items;
322 geometric_object *os = o->subclass.compound_geometric_object_data->component_objects.items;
323 vector3 shiftby = o->center;
324 for (i = 0; i < n; ++i) {
325 *o = os[i];
326 o->center = vector3_plus(o->center, shiftby);
327 if (point_in_fixed_pobjectp(p, o)) return 1;
328 }
329 break;
330 }
331 }
332 return 0;
333 }
334
335 /**************************************************************************/
336
337 /* convert a point p inside o to a coordinate in [0,1]^3 that
338 is some "natural" coordinate for the object */
to_geom_object_coords(vector3 p,geometric_object o)339 vector3 to_geom_object_coords(vector3 p, geometric_object o) {
340 const vector3 half = {0.5, 0.5, 0.5};
341 vector3 r = vector3_minus(p, o.center);
342
343 switch (o.which_subclass) {
344 default: {
345 vector3 po = {0, 0, 0};
346 return po;
347 }
348 case GEOM SPHERE: {
349 number radius = o.subclass.sphere_data->radius;
350 return vector3_plus(half, vector3_scale(0.5 / radius, r));
351 }
352 /* case GEOM CYLINDER:
353 NOT YET IMPLEMENTED */
354 case GEOM BLOCK: {
355 vector3 proj = matrix3x3_vector3_mult(o.subclass.block_data->projection_matrix, r);
356 vector3 size = o.subclass.block_data->size;
357 if (size.x != 0.0) proj.x /= size.x;
358 if (size.y != 0.0) proj.y /= size.y;
359 if (size.z != 0.0) proj.z /= size.z;
360 return vector3_plus(half, proj);
361 }
362 /* case GEOM PRISM:
363 NOT YET IMPLEMENTED */
364 }
365 }
366
367 /* inverse of to_geom_object_coords */
from_geom_object_coords(vector3 p,geometric_object o)368 vector3 from_geom_object_coords(vector3 p, geometric_object o) {
369 const vector3 half = {0.5, 0.5, 0.5};
370 p = vector3_minus(p, half);
371 switch (o.which_subclass) {
372 default: return o.center;
373 case GEOM SPHERE: {
374 number radius = o.subclass.sphere_data->radius;
375 return vector3_plus(o.center, vector3_scale(radius / 0.5, p));
376 }
377 /* case GEOM CYLINDER:
378 NOT YET IMPLEMENTED */
379 case GEOM BLOCK: {
380 vector3 size = o.subclass.block_data->size;
381 return vector3_plus(
382 o.center,
383 vector3_plus(vector3_scale(size.x * p.x, o.subclass.block_data->e1),
384 vector3_plus(vector3_scale(size.y * p.y, o.subclass.block_data->e2),
385 vector3_scale(size.z * p.z, o.subclass.block_data->e3))));
386 }
387 /* case GEOM PRISM:
388 NOT YET IMPLEMENTED */
389 }
390 }
391
392 /**************************************************************************/
393 /* Return the normal vector from the given object to the given point,
394 in lattice coordinates, using the surface of the object that the
395 point is "closest" to for some definition of "closest" that is
396 reasonable (at least for points near to the object). The length and
397 sign of the normal vector are arbitrary. */
398
normal_to_object(vector3 p,geometric_object o)399 vector3 CTLIO normal_to_object(vector3 p, geometric_object o) {
400 geom_fix_object_ptr(&o);
401 return normal_to_fixed_object(p, o);
402 }
403
normal_to_fixed_object(vector3 p,geometric_object o)404 vector3 normal_to_fixed_object(vector3 p, geometric_object o) {
405 vector3 r = vector3_minus(p, o.center);
406
407 switch (o.which_subclass) {
408
409 case GEOM CYLINDER: {
410 vector3 rm = matrix3x3_vector3_mult(geometry_lattice.metric, r);
411 double proj = vector3_dot(o.subclass.cylinder_data->axis, rm),
412 height = o.subclass.cylinder_data->height, radius, prad;
413 if (fabs(proj) > height * 0.5) return o.subclass.cylinder_data->axis;
414 radius = o.subclass.cylinder_data->radius;
415 prad = sqrt(fabs(vector3_dot(r, rm) - proj * proj));
416 if (o.subclass.cylinder_data->which_subclass == CYL CONE)
417 radius += (proj / height + 0.5) *
418 (o.subclass.cylinder_data->subclass.cone_data->radius2 - radius);
419 if (fabs(fabs(proj) - height * 0.5) < fabs(prad - radius))
420 return o.subclass.cylinder_data->axis;
421 if (o.subclass.cylinder_data->which_subclass == CYL CONE)
422 return vector3_minus(
423 r, vector3_scale(
424 proj + prad * (o.subclass.cylinder_data->subclass.cone_data->radius2 - radius) /
425 height,
426 o.subclass.cylinder_data->axis));
427 else
428 return vector3_minus(r, vector3_scale(proj, o.subclass.cylinder_data->axis));
429 } // case GEOM CYLINDER
430
431 case GEOM BLOCK: {
432 vector3 proj = matrix3x3_vector3_mult(o.subclass.block_data->projection_matrix, r);
433 switch (o.subclass.block_data->which_subclass) {
434 case BLK BLOCK_SELF: {
435 vector3 size = o.subclass.block_data->size;
436 double d1 = fabs(fabs(proj.x) - 0.5 * size.x);
437 double d2 = fabs(fabs(proj.y) - 0.5 * size.y);
438 double d3 = fabs(fabs(proj.z) - 0.5 * size.z);
439 if (d1 < d2 && d1 < d3)
440 return matrix3x3_row1(o.subclass.block_data->projection_matrix);
441 else if (d2 < d3)
442 return matrix3x3_row2(o.subclass.block_data->projection_matrix);
443 else
444 return matrix3x3_row3(o.subclass.block_data->projection_matrix);
445 } // case BLK BLOCK_SELF
446
447 case BLK ELLIPSOID:
448 default: {
449 vector3 isa = o.subclass.block_data->subclass.ellipsoid_data->inverse_semi_axes;
450 proj.x *= isa.x * isa.x;
451 proj.y *= isa.y * isa.y;
452 proj.z *= isa.z * isa.z;
453 return matrix3x3_transpose_vector3_mult(o.subclass.block_data->projection_matrix, proj);
454 } // case BLK ELLIPSOID
455
456 } // switch (o.subclass.block_data->which_subclass)
457
458 } // case GEOM BLOCK
459
460 case GEOM PRISM: return normal_to_prism(o.subclass.prism_data, p);
461
462 default: return r;
463 } // switch (o.which_subclass)
464
465 return r; // never get here
466 }
467
468 /**************************************************************************/
469
470 /* Here is a useful macro to loop over different possible shifts of
471 the lattice vectors. body is executed for each possible shift,
472 where the shift is given by the value of shiftby (which should
473 be a vector3 variable). I would much rather make this a function,
474 but C's lack of lambda-like function construction or closures makes
475 this easier to do as a macro. (One could at least wish for
476 an easier way to make multi-line macros.) */
477
478 #define LOOP_PERIODIC(shiftby, body) \
479 { \
480 switch (dimensions) { \
481 case 1: { \
482 int iii; \
483 shiftby.y = shiftby.z = 0; \
484 for (iii = -1; iii <= 1; ++iii) { \
485 shiftby.x = iii * geometry_lattice.size.x; \
486 body; \
487 } \
488 break; \
489 } \
490 case 2: { \
491 int iii, jjj; \
492 shiftby.z = 0; \
493 for (iii = -1; iii <= 1; ++iii) { \
494 shiftby.x = iii * geometry_lattice.size.x; \
495 for (jjj = -1; jjj <= 1; ++jjj) { \
496 shiftby.y = jjj * geometry_lattice.size.y; \
497 body; \
498 } \
499 } \
500 break; \
501 } \
502 case 3: { \
503 int iii, jjj, kkk; \
504 for (iii = -1; iii <= 1; ++iii) { \
505 shiftby.x = iii * geometry_lattice.size.x; \
506 for (jjj = -1; jjj <= 1; ++jjj) { \
507 shiftby.y = jjj * geometry_lattice.size.y; \
508 for (kkk = -1; kkk <= 1; ++kkk) { \
509 shiftby.z = kkk * geometry_lattice.size.z; \
510 body; \
511 if (geometry_lattice.size.z == 0) break; \
512 } \
513 if (geometry_lattice.size.y == 0) break; \
514 } \
515 if (geometry_lattice.size.x == 0) break; \
516 } \
517 break; \
518 } \
519 } \
520 }
521
522 /**************************************************************************/
523
524 /* Like point_in_objectp, but also checks the object shifted
525 by the lattice vectors: */
526
point_in_periodic_objectp(vector3 p,geometric_object o)527 boolean CTLIO point_in_periodic_objectp(vector3 p, geometric_object o) {
528 geom_fix_object_ptr(&o);
529 return point_in_periodic_fixed_objectp(p, o);
530 }
531
point_in_periodic_fixed_objectp(vector3 p,geometric_object o)532 boolean point_in_periodic_fixed_objectp(vector3 p, geometric_object o) {
533 vector3 shiftby;
534 LOOP_PERIODIC(shiftby, if (point_in_fixed_objectp(vector3_minus(p, shiftby), o)) return 1);
535 return 0;
536 }
537
point_shift_in_periodic_fixed_pobjectp(vector3 p,geometric_object * o,vector3 * shiftby)538 boolean point_shift_in_periodic_fixed_pobjectp(vector3 p, geometric_object *o, vector3 *shiftby) {
539 geometric_object o0 = *o;
540 LOOP_PERIODIC((*shiftby), {
541 *o = o0;
542 if (point_in_fixed_pobjectp(vector3_minus(p, *shiftby), o)) return 1;
543 });
544 return 0;
545 }
546
547 /**************************************************************************/
548
549 /* Functions to return the object or material type corresponding to
550 the point p (in the lattice basis). Returns default_material if p
551 is not in any object.
552
553 Requires that the global input vars geometry_lattice, geometry,
554 dimensions, default_material and ensure_periodicity already be
555 initialized.
556
557 Also requires that geom_fix_objects() has been called!
558
559 material_of_point_inobject is a variant that also returns whether
560 or not the point was in any object. */
561
object_of_point0(geometric_object_list geometry,vector3 p,vector3 * shiftby)562 geometric_object object_of_point0(geometric_object_list geometry, vector3 p, vector3 *shiftby) {
563 geometric_object o;
564 int index;
565 shiftby->x = shiftby->y = shiftby->z = 0;
566 /* loop in reverse order so that later items are given precedence: */
567 for (index = geometry.num_items - 1; index >= 0; --index) {
568 o = geometry.items[index];
569 if ((ensure_periodicity && point_shift_in_periodic_fixed_pobjectp(p, &o, shiftby)) ||
570 point_in_fixed_pobjectp(p, &o))
571 return o;
572 }
573 o.which_subclass = GEOM GEOMETRIC_OBJECT_SELF; /* no object found */
574 return o;
575 }
576
object_of_point(vector3 p,vector3 * shiftby)577 geometric_object object_of_point(vector3 p, vector3 *shiftby) {
578 return object_of_point0(geometry, p, shiftby);
579 }
580
material_of_point_inobject0(geometric_object_list geometry,vector3 p,boolean * inobject)581 material_type material_of_point_inobject0(geometric_object_list geometry, vector3 p,
582 boolean *inobject) {
583 vector3 shiftby;
584 geometric_object o = object_of_point0(geometry, p, &shiftby);
585 *inobject = o.which_subclass != GEOM GEOMETRIC_OBJECT_SELF;
586 ;
587 return (*inobject ? o.material : default_material);
588 }
589
material_of_point_inobject(vector3 p,boolean * inobject)590 material_type material_of_point_inobject(vector3 p, boolean *inobject) {
591 return material_of_point_inobject0(geometry, p, inobject);
592 }
593
material_of_point0(geometric_object_list geometry,vector3 p)594 material_type material_of_point0(geometric_object_list geometry, vector3 p) {
595 boolean inobject;
596 return material_of_point_inobject0(geometry, p, &inobject);
597 }
598
material_of_point(vector3 p)599 material_type material_of_point(vector3 p) { return material_of_point0(geometry, p); }
600
601 /**************************************************************************/
602
603 /* Given a geometric object o, display some information about it,
604 indented by indentby spaces. */
605
display_geometric_object_info(int indentby,geometric_object o)606 void CTLIO display_geometric_object_info(int indentby, geometric_object o) {
607 geom_fix_object_ptr(&o);
608 ctl_printf("%*s", indentby, "");
609 switch (o.which_subclass) {
610 case GEOM CYLINDER:
611 switch (o.subclass.cylinder_data->which_subclass) {
612 case CYL WEDGE: ctl_printf("wedge"); break;
613 case CYL CONE: ctl_printf("cone"); break;
614 case CYL CYLINDER_SELF: ctl_printf("cylinder"); break;
615 }
616 break;
617 case GEOM SPHERE: ctl_printf("sphere"); break;
618 case GEOM BLOCK:
619 switch (o.subclass.block_data->which_subclass) {
620 case BLK ELLIPSOID: ctl_printf("ellipsoid"); break;
621 case BLK BLOCK_SELF: ctl_printf("block"); break;
622 }
623 break;
624 case GEOM PRISM: ctl_printf("prism"); break;
625 case GEOM COMPOUND_GEOMETRIC_OBJECT: ctl_printf("compound object"); break;
626 default: ctl_printf("geometric object"); break;
627 }
628 ctl_printf(", center = (%g,%g,%g)\n", o.center.x, o.center.y, o.center.z);
629 switch (o.which_subclass) {
630 case GEOM CYLINDER:
631 ctl_printf("%*s radius %g, height %g, axis (%g, %g, %g)\n", indentby, "",
632 o.subclass.cylinder_data->radius, o.subclass.cylinder_data->height,
633 o.subclass.cylinder_data->axis.x, o.subclass.cylinder_data->axis.y,
634 o.subclass.cylinder_data->axis.z);
635 if (o.subclass.cylinder_data->which_subclass == CYL CONE)
636 ctl_printf("%*s radius2 %g\n", indentby, "",
637 o.subclass.cylinder_data->subclass.cone_data->radius2);
638 else if (o.subclass.cylinder_data->which_subclass == CYL WEDGE)
639 ctl_printf("%*s wedge-theta %g\n", indentby, "",
640 o.subclass.cylinder_data->subclass.wedge_data->wedge_angle);
641 break;
642 case GEOM SPHERE:
643 ctl_printf("%*s radius %g\n", indentby, "", o.subclass.sphere_data->radius);
644 break;
645 case GEOM BLOCK:
646 ctl_printf("%*s size (%g,%g,%g)\n", indentby, "", o.subclass.block_data->size.x,
647 o.subclass.block_data->size.y, o.subclass.block_data->size.z);
648 ctl_printf(
649 "%*s axes (%g,%g,%g), (%g,%g,%g), (%g,%g,%g)\n", indentby, "",
650 o.subclass.block_data->e1.x, o.subclass.block_data->e1.y, o.subclass.block_data->e1.z,
651 o.subclass.block_data->e2.x, o.subclass.block_data->e2.y, o.subclass.block_data->e2.z,
652 o.subclass.block_data->e3.x, o.subclass.block_data->e3.y, o.subclass.block_data->e3.z);
653 break;
654 case GEOM PRISM: display_prism_info(indentby, &o); break;
655 case GEOM COMPOUND_GEOMETRIC_OBJECT: {
656 int i;
657 int n = o.subclass.compound_geometric_object_data->component_objects.num_items;
658 geometric_object *os = o.subclass.compound_geometric_object_data->component_objects.items;
659 ctl_printf("%*s %d components:\n", indentby, "", n);
660 for (i = 0; i < n; ++i)
661 display_geometric_object_info(indentby + 5, os[i]);
662 break;
663 }
664 default: break;
665 }
666 }
667
668 /**************************************************************************/
669
670 /* Compute the intersections with o of a line along p+s*d, returning
671 the number of intersections (at most 2) and the two intersection "s"
672 values in s[0] and s[1]. (Note: o must not be a compound object.) */
intersect_line_with_object(vector3 p,vector3 d,geometric_object o,double s[2])673 int intersect_line_with_object(vector3 p, vector3 d, geometric_object o, double s[2]) {
674 p = vector3_minus(p, o.center);
675 s[0] = s[1] = 0;
676 switch (o.which_subclass) {
677 case GEOM SPHERE: {
678 number radius = o.subclass.sphere_data->radius;
679 vector3 dm = matrix3x3_vector3_mult(geometry_lattice.metric, d);
680 double a = vector3_dot(d, dm);
681 double b2 = -vector3_dot(dm, p);
682 double c =
683 vector3_dot(p, matrix3x3_vector3_mult(geometry_lattice.metric, p)) - radius * radius;
684 double discrim = b2 * b2 - a * c;
685 if (discrim < 0)
686 return 0;
687 else if (discrim == 0) {
688 s[0] = b2 / a;
689 return 1;
690 }
691 else {
692 discrim = sqrt(discrim);
693 s[0] = (b2 + discrim) / a;
694 s[1] = (b2 - discrim) / a;
695 return 2;
696 }
697 } // case GEOM SPHERE
698 case GEOM CYLINDER: {
699 vector3 dm = matrix3x3_vector3_mult(geometry_lattice.metric, d);
700 vector3 pm = matrix3x3_vector3_mult(geometry_lattice.metric, p);
701 number height = o.subclass.cylinder_data->height;
702 number radius = o.subclass.cylinder_data->radius;
703 number radius2 = o.subclass.cylinder_data->which_subclass == CYL CONE
704 ? o.subclass.cylinder_data->subclass.cone_data->radius2
705 : radius;
706 double dproj = vector3_dot(o.subclass.cylinder_data->axis, dm);
707 double pproj = vector3_dot(o.subclass.cylinder_data->axis, pm);
708 double D = (radius2 - radius) / height;
709 double L = radius + (radius2 - radius) * 0.5 + pproj * D;
710 double a = vector3_dot(d, dm) - dproj * dproj * (1 + D * D);
711 double b2 = dproj * (pproj + D * L) - vector3_dot(p, dm);
712 double c = vector3_dot(p, pm) - pproj * pproj - L * L;
713 double discrim = b2 * b2 - a * c;
714 int ret;
715 if (a == 0) { /* linear equation */
716 if (b2 == 0) {
717 if (c == 0) { /* infinite intersections */
718 s[0] = ((height * 0.5) - pproj) / dproj;
719 s[1] = -((height * 0.5) + pproj) / dproj;
720 return 2;
721 }
722 else
723 ret = 0;
724 }
725 else {
726 s[0] = 0.5 * c / b2;
727 ret = 1;
728 }
729 }
730 else if (discrim < 0)
731 ret = 0;
732 else if (discrim == 0) {
733 s[0] = b2 / a;
734 ret = 1;
735 }
736 else {
737 discrim = sqrt(discrim);
738 s[0] = (b2 + discrim) / a;
739 s[1] = (b2 - discrim) / a;
740 ret = 2;
741 }
742 if (ret == 2 && fabs(pproj + s[1] * dproj) > height * 0.5) ret = 1;
743 if (ret >= 1 && fabs(pproj + s[0] * dproj) > height * 0.5) {
744 --ret;
745 s[0] = s[1];
746 }
747 if (ret == 2 || dproj == 0) return ret;
748 /* find intersections with endcaps */
749 s[ret] = (height * 0.5 - pproj) / dproj;
750 if (a * s[ret] * s[ret] - 2 * b2 * s[ret] + c <= 0) ++ret;
751 if (ret < 2) {
752 s[ret] = -(height * 0.5 + pproj) / dproj;
753 if (a * s[ret] * s[ret] - 2 * b2 * s[ret] + c <= 0) ++ret;
754 }
755 if (ret == 2 && s[0] == s[1]) ret = 1;
756 return ret;
757 } // case GEOM CYLINDER
758 case GEOM BLOCK: {
759 vector3 dproj = matrix3x3_vector3_mult(o.subclass.block_data->projection_matrix, d);
760 vector3 pproj = matrix3x3_vector3_mult(o.subclass.block_data->projection_matrix, p);
761 switch (o.subclass.block_data->which_subclass) {
762 case BLK BLOCK_SELF: {
763 vector3 size = o.subclass.block_data->size;
764 int ret = 0;
765 size.x *= 0.5;
766 size.y *= 0.5;
767 size.z *= 0.5;
768 if (dproj.x != 0) {
769 s[ret] = (size.x - pproj.x) / dproj.x;
770 if (fabs(pproj.y + s[ret] * dproj.y) <= size.y &&
771 fabs(pproj.z + s[ret] * dproj.z) <= size.z)
772 ++ret;
773 s[ret] = (-size.x - pproj.x) / dproj.x;
774 if (fabs(pproj.y + s[ret] * dproj.y) <= size.y &&
775 fabs(pproj.z + s[ret] * dproj.z) <= size.z)
776 ++ret;
777 if (ret == 2) return 2;
778 }
779 if (dproj.y != 0) {
780 s[ret] = (size.y - pproj.y) / dproj.y;
781 if (fabs(pproj.x + s[ret] * dproj.x) <= size.x &&
782 fabs(pproj.z + s[ret] * dproj.z) <= size.z)
783 ++ret;
784 if (ret == 2) return 2;
785 s[ret] = (-size.y - pproj.y) / dproj.y;
786 if (fabs(pproj.x + s[ret] * dproj.x) <= size.x &&
787 fabs(pproj.z + s[ret] * dproj.z) <= size.z)
788 ++ret;
789 if (ret == 2) return 2;
790 }
791 if (dproj.z != 0) {
792 s[ret] = (size.z - pproj.z) / dproj.z;
793 if (fabs(pproj.x + s[ret] * dproj.x) <= size.x &&
794 fabs(pproj.y + s[ret] * dproj.y) <= size.y)
795 ++ret;
796 if (ret == 2) return 2;
797 s[ret] = (-size.z - pproj.z) / dproj.z;
798 if (fabs(pproj.x + s[ret] * dproj.x) <= size.x &&
799 fabs(pproj.y + s[ret] * dproj.y) <= size.y)
800 ++ret;
801 }
802 return ret;
803 } // case BLK BLOCK_SELF:
804
805 case BLK ELLIPSOID:
806 default: {
807 vector3 isa = o.subclass.block_data->subclass.ellipsoid_data->inverse_semi_axes;
808 double a, b2, c, discrim;
809 dproj.x *= isa.x;
810 dproj.y *= isa.y;
811 dproj.z *= isa.z;
812 pproj.x *= isa.x;
813 pproj.y *= isa.y;
814 pproj.z *= isa.z;
815 a = vector3_dot(dproj, dproj);
816 b2 = -vector3_dot(dproj, pproj);
817 c = vector3_dot(pproj, pproj) - 1;
818 discrim = b2 * b2 - a * c;
819 if (discrim < 0)
820 return 0;
821 else if (discrim == 0) {
822 s[0] = b2 / a;
823 return 1;
824 }
825 else {
826 discrim = sqrt(discrim);
827 s[0] = (b2 + discrim) / a;
828 s[1] = (b2 - discrim) / a;
829 return 2;
830 }
831 } // case BLK BLOCK_SELF, default
832
833 } // switch (o.subclass.block_data->which_subclass)
834
835 } // case GEOM BLOCK
836 default: return 0;
837 }
838 }
839
840 /* Compute the intersections with o of a line along p+s*d in the interval s in [a,b], returning
841 the length of the s intersection in this interval. (Note: o must not be a compound object.) */
intersect_line_segment_with_object(vector3 p,vector3 d,geometric_object o,double a,double b)842 double intersect_line_segment_with_object(vector3 p, vector3 d, geometric_object o, double a,
843 double b) {
844 if (o.which_subclass == GEOM PRISM) {
845 return intersect_line_segment_with_prism(o.subclass.prism_data, p, d, a, b);
846 }
847 else {
848 double s[2];
849 if (2 == intersect_line_with_object(p, d, o, s)) {
850 double ds = (s[0] < s[1] ? MIN(s[1], b) - MAX(s[0], a) : MIN(s[0], b) - MAX(s[1], a));
851 return (ds > 0 ? ds : 0.0);
852 }
853 else
854 return 0.0;
855 }
856 }
857
858 /**************************************************************************/
859
860 /* Given a basis (matrix columns are the basis unit vectors) and the
861 size of the lattice (in basis vectors), returns a new "square"
862 basis. This corresponds to a region of the same volume, but made
863 rectangular, suitable for outputing to an HDF file.
864
865 Given a vector in the range (0..1, 0..1, 0..1), multiplying by
866 the square basis matrix will yield the coordinates of a point
867 in the rectangular volume, given in the lattice basis. */
868
square_basis(matrix3x3 basis,vector3 size)869 matrix3x3 CTLIO square_basis(matrix3x3 basis, vector3 size) {
870 matrix3x3 square;
871
872 square.c0 = basis.c0;
873
874 square.c1 = vector3_minus(basis.c1, vector3_scale(vector3_dot(basis.c0, basis.c1), basis.c1));
875
876 square.c2 = vector3_minus(basis.c2, vector3_scale(vector3_dot(basis.c0, basis.c2), basis.c2));
877 square.c2 = vector3_minus(
878 square.c2, vector3_scale(vector3_dot(basis.c0, square.c2), unit_vector3(square.c2)));
879
880 square.c0 = vector3_scale(size.x, square.c0);
881 square.c1 = vector3_scale(size.y, square.c1);
882 square.c2 = vector3_scale(size.z, square.c2);
883
884 return matrix3x3_mult(matrix3x3_inverse(basis), square);
885 }
886
887 /**************************************************************************/
888
889 /* compute the 3d volume enclosed by a geometric object o. */
890
geom_object_volume(GEOMETRIC_OBJECT o)891 double geom_object_volume(GEOMETRIC_OBJECT o) {
892 switch (o.which_subclass) {
893 case GEOM SPHERE: {
894 number radius = o.subclass.sphere_data->radius;
895 return (1.333333333333333333 * K_PI) * radius * radius * radius;
896 }
897 case GEOM CYLINDER: {
898 number height = o.subclass.cylinder_data->height;
899 number radius = o.subclass.cylinder_data->radius;
900 number radius2 = o.subclass.cylinder_data->which_subclass == CYL CONE
901 ? o.subclass.cylinder_data->subclass.cone_data->radius2
902 : radius;
903 double vol = height * (K_PI / 3) * (radius * radius + radius * radius2 + radius2 * radius2);
904 if (o.subclass.cylinder_data->which_subclass == CYL WEDGE)
905 return vol * fabs(o.subclass.cylinder_data->subclass.wedge_data->wedge_angle) / (2 * K_PI);
906 else
907 return vol;
908 }
909 case GEOM BLOCK: {
910 vector3 size = o.subclass.block_data->size;
911 double vol = size.x * size.y * size.z *
912 fabs(matrix3x3_determinant(geometry_lattice.basis) /
913 matrix3x3_determinant(o.subclass.block_data->projection_matrix));
914 return o.subclass.block_data->which_subclass == BLK BLOCK_SELF ? vol : vol * (K_PI / 6);
915 }
916 case GEOM PRISM: {
917 return get_prism_volume(o.subclass.prism_data);
918 }
919 default: return 0; /* unsupported object types? */
920 }
921 }
922
923 /**************************************************************************/
924 /**************************************************************************/
925
926 /* Fast geometry routines */
927
928 /* Using the above material_of_point routine is way too slow, especially
929 when there are lots of objects to test. Thus, we develop the following
930 replacement routines.
931
932 The basic idea here is twofold. (1) Compute bounding boxes for
933 each geometric object, for which inclusion tests can be computed
934 quickly. (2) Build a tree that recursively breaks down the unit cell
935 in half, allowing us to perform searches in logarithmic time. */
936
937 /**************************************************************************/
938
939 /* geom_box utilities: */
geom_box_union(geom_box * bu,const geom_box * b1,const geom_box * b2)940 static void geom_box_union(geom_box *bu, const geom_box *b1, const geom_box *b2) {
941 bu->low.x = MIN(b1->low.x, b2->low.x);
942 bu->low.y = MIN(b1->low.y, b2->low.y);
943 bu->low.z = MIN(b1->low.z, b2->low.z);
944 bu->high.x = MAX(b1->high.x, b2->high.x);
945 bu->high.y = MAX(b1->high.y, b2->high.y);
946 bu->high.z = MAX(b1->high.z, b2->high.z);
947 }
948
geom_box_intersection(geom_box * bi,const geom_box * b1,const geom_box * b2)949 static void geom_box_intersection(geom_box *bi, const geom_box *b1, const geom_box *b2) {
950 bi->low.x = MAX(b1->low.x, b2->low.x);
951 bi->low.y = MAX(b1->low.y, b2->low.y);
952 bi->low.z = MAX(b1->low.z, b2->low.z);
953 bi->high.x = MIN(b1->high.x, b2->high.x);
954 bi->high.y = MIN(b1->high.y, b2->high.y);
955 bi->high.z = MIN(b1->high.z, b2->high.z);
956 }
957
geom_box_add_pt(geom_box * b,vector3 p)958 static void geom_box_add_pt(geom_box *b, vector3 p) {
959 b->low.x = MIN(b->low.x, p.x);
960 b->low.y = MIN(b->low.y, p.y);
961 b->low.z = MIN(b->low.z, p.z);
962 b->high.x = MAX(b->high.x, p.x);
963 b->high.y = MAX(b->high.y, p.y);
964 b->high.z = MAX(b->high.z, p.z);
965 }
966
967 #define BETWEEN(x, low, high) ((x) >= (low) && (x) <= (high))
968
geom_box_contains_point(const geom_box * b,vector3 p)969 static int geom_box_contains_point(const geom_box *b, vector3 p) {
970 return (BETWEEN(p.x, b->low.x, b->high.x) && BETWEEN(p.y, b->low.y, b->high.y) &&
971 BETWEEN(p.z, b->low.z, b->high.z));
972 }
973
974 /* return whether or not the given two boxes intersect */
geom_boxes_intersect(const geom_box * b1,const geom_box * b2)975 static int geom_boxes_intersect(const geom_box *b1, const geom_box *b2) {
976 /* true if the x, y, and z ranges all intersect. */
977 return (
978 (BETWEEN(b1->low.x, b2->low.x, b2->high.x) || BETWEEN(b1->high.x, b2->low.x, b2->high.x) ||
979 BETWEEN(b2->low.x, b1->low.x, b1->high.x)) &&
980 (BETWEEN(b1->low.y, b2->low.y, b2->high.y) || BETWEEN(b1->high.y, b2->low.y, b2->high.y) ||
981 BETWEEN(b2->low.y, b1->low.y, b1->high.y)) &&
982 (BETWEEN(b1->low.z, b2->low.z, b2->high.z) || BETWEEN(b1->high.z, b2->low.z, b2->high.z) ||
983 BETWEEN(b2->low.z, b1->low.z, b1->high.z)));
984 }
985
geom_box_shift(geom_box * b,vector3 shiftby)986 static void geom_box_shift(geom_box *b, vector3 shiftby) {
987 b->low = vector3_plus(b->low, shiftby);
988 b->high = vector3_plus(b->high, shiftby);
989 }
990
991 /**************************************************************************/
992
993 /* Computing a bounding box for a geometric object: */
994
995 /* compute | (b x c) / (a * (b x c)) |, for use below */
compute_dot_cross(vector3 a,vector3 b,vector3 c)996 static number compute_dot_cross(vector3 a, vector3 b, vector3 c) {
997 vector3 bxc = vector3_cross(b, c);
998 return fabs(vector3_norm(bxc) / vector3_dot(a, bxc));
999 }
1000
1001 /* Compute a bounding box for the object o, preferably the smallest
1002 bounding box. The box is a parallelepiped with axes given by
1003 the geometry lattice vectors, and its corners are given in the
1004 lattice basis.
1005
1006 Requires that geometry_lattice global has been initialized,
1007 etcetera. */
geom_get_bounding_box(geometric_object o,geom_box * box)1008 void geom_get_bounding_box(geometric_object o, geom_box *box) {
1009 geom_fix_object_ptr(&o);
1010
1011 /* initialize to empty box at the center of the object: */
1012 box->low = box->high = o.center;
1013
1014 switch (o.which_subclass) {
1015 case GEOM GEOMETRIC_OBJECT_SELF: break;
1016 case GEOM SPHERE: {
1017 /* Find the parallelepiped that the sphere inscribes.
1018 The math comes out surpisingly simple--try it! */
1019
1020 number radius = o.subclass.sphere_data->radius;
1021 /* actually, we could achieve the same effect here
1022 by inverting the geometry_lattice.basis matrix... */
1023 number r1 =
1024 compute_dot_cross(geometry_lattice.b1, geometry_lattice.b2, geometry_lattice.b3) * radius;
1025 number r2 =
1026 compute_dot_cross(geometry_lattice.b2, geometry_lattice.b3, geometry_lattice.b1) * radius;
1027 number r3 =
1028 compute_dot_cross(geometry_lattice.b3, geometry_lattice.b1, geometry_lattice.b2) * radius;
1029 box->low.x -= r1;
1030 box->low.y -= r2;
1031 box->low.z -= r3;
1032 box->high.x += r1;
1033 box->high.y += r2;
1034 box->high.z += r3;
1035 break;
1036 }
1037 case GEOM CYLINDER: {
1038 /* Find the bounding boxes of the two (circular) ends of
1039 the cylinder, then take the union. Again, the math
1040 for finding the bounding parallelepiped of a circle
1041 comes out suprisingly simple in the end. Proof left
1042 as an exercise for the reader. */
1043
1044 number radius = o.subclass.cylinder_data->radius;
1045 number h = o.subclass.cylinder_data->height * 0.5;
1046 vector3 axis = /* cylinder axis in cartesian coords */
1047 matrix3x3_vector3_mult(geometry_lattice.basis, o.subclass.cylinder_data->axis);
1048 vector3 e12 = vector3_cross(geometry_lattice.basis1, geometry_lattice.basis2);
1049 vector3 e23 = vector3_cross(geometry_lattice.basis2, geometry_lattice.basis3);
1050 vector3 e31 = vector3_cross(geometry_lattice.basis3, geometry_lattice.basis1);
1051 number elen2, eproj;
1052 number r1, r2, r3;
1053 geom_box tmp_box;
1054
1055 /* Find bounding box dimensions, in lattice coords,
1056 for the circular ends of the cylinder: */
1057
1058 elen2 = vector3_dot(e23, e23);
1059 eproj = vector3_dot(e23, axis);
1060 r1 = fabs(sqrt(fabs(elen2 - eproj * eproj)) / vector3_dot(e23, geometry_lattice.b1));
1061
1062 elen2 = vector3_dot(e31, e31);
1063 eproj = vector3_dot(e31, axis);
1064 r2 = fabs(sqrt(fabs(elen2 - eproj * eproj)) / vector3_dot(e31, geometry_lattice.b2));
1065
1066 elen2 = vector3_dot(e12, e12);
1067 eproj = vector3_dot(e12, axis);
1068 r3 = fabs(sqrt(fabs(elen2 - eproj * eproj)) / vector3_dot(e12, geometry_lattice.b3));
1069
1070 /* Get axis in lattice coords: */
1071 axis = o.subclass.cylinder_data->axis;
1072
1073 tmp_box = *box; /* set tmp_box to center of object */
1074
1075 /* bounding box for -h*axis cylinder end: */
1076 box->low.x -= h * axis.x + r1 * radius;
1077 box->low.y -= h * axis.y + r2 * radius;
1078 box->low.z -= h * axis.z + r3 * radius;
1079 box->high.x -= h * axis.x - r1 * radius;
1080 box->high.y -= h * axis.y - r2 * radius;
1081 box->high.z -= h * axis.z - r3 * radius;
1082
1083 if (o.subclass.cylinder_data->which_subclass == CYL CONE)
1084 radius = fabs(o.subclass.cylinder_data->subclass.cone_data->radius2);
1085
1086 /* bounding box for +h*axis cylinder end: */
1087 tmp_box.low.x += h * axis.x - r1 * radius;
1088 tmp_box.low.y += h * axis.y - r2 * radius;
1089 tmp_box.low.z += h * axis.z - r3 * radius;
1090 tmp_box.high.x += h * axis.x + r1 * radius;
1091 tmp_box.high.y += h * axis.y + r2 * radius;
1092 tmp_box.high.z += h * axis.z + r3 * radius;
1093
1094 geom_box_union(box, box, &tmp_box);
1095 break;
1096 }
1097 case GEOM BLOCK: {
1098 /* blocks are easy: just enlarge the box to be big enough to
1099 contain all 8 corners of the block. */
1100
1101 vector3 s1 = vector3_scale(o.subclass.block_data->size.x, o.subclass.block_data->e1);
1102 vector3 s2 = vector3_scale(o.subclass.block_data->size.y, o.subclass.block_data->e2);
1103 vector3 s3 = vector3_scale(o.subclass.block_data->size.z, o.subclass.block_data->e3);
1104 vector3 corner =
1105 vector3_plus(o.center, vector3_scale(-0.5, vector3_plus(s1, vector3_plus(s2, s3))));
1106
1107 geom_box_add_pt(box, corner);
1108 geom_box_add_pt(box, vector3_plus(corner, s1));
1109 geom_box_add_pt(box, vector3_plus(corner, s2));
1110 geom_box_add_pt(box, vector3_plus(corner, s3));
1111 geom_box_add_pt(box, vector3_plus(corner, vector3_plus(s1, s2)));
1112 geom_box_add_pt(box, vector3_plus(corner, vector3_plus(s1, s3)));
1113 geom_box_add_pt(box, vector3_plus(corner, vector3_plus(s3, s2)));
1114 geom_box_add_pt(box, vector3_plus(corner, vector3_plus(s1, vector3_plus(s2, s3))));
1115 break;
1116 }
1117 case GEOM PRISM: {
1118 get_prism_bounding_box(o.subclass.prism_data, box);
1119 break;
1120 }
1121 case GEOM COMPOUND_GEOMETRIC_OBJECT: {
1122 int i;
1123 int n = o.subclass.compound_geometric_object_data->component_objects.num_items;
1124 geometric_object *os = o.subclass.compound_geometric_object_data->component_objects.items;
1125 for (i = 0; i < n; ++i) {
1126 geom_box boxi;
1127 geom_get_bounding_box(os[i], &boxi);
1128 geom_box_shift(&boxi, o.center);
1129 geom_box_union(box, box, &boxi);
1130 }
1131 break;
1132 }
1133 }
1134 }
1135
1136 /**************************************************************************/
1137 /* Compute the fraction of a box's volume (or area/length in 2d/1d) that
1138 overlaps an object. Instead of a box, we also allow an ellipsoid
1139 inscribed inside the box (or a skewed ellipsoid if the box is not
1140 orthogonal). */
1141
1142 typedef struct {
1143 geometric_object o;
1144 vector3 p, dir;
1145 int pdim[2]; /* the (up to two) integration directions */
1146 double scx[2]; /* scale factor (e.g. sign flip) for x coordinates */
1147 unsigned dim;
1148 double a0, b0; /* box limits along analytic direction */
1149 int is_ellipsoid; /* 0 for box, 1 for ellipsoid */
1150 double winv[2], c[2]; /* ellipsoid width-inverses/centers in int. dirs */
1151 double w0, c0; /* width/center along analytic direction */
1152 } overlap_data;
1153
overlap_integrand(integer ndim,number * x,void * data_)1154 static double overlap_integrand(integer ndim, number *x, void *data_) {
1155 overlap_data *data = (overlap_data *)data_;
1156 double s[2];
1157 const double *scx = data->scx;
1158 vector3 p = data->p;
1159 double a0 = data->a0, b0 = data->b0;
1160 double scale_result = 1.0;
1161
1162 if (ndim > 0) {
1163 switch (data->pdim[0]) {
1164 case 0: p.x = scx[0] * x[0]; break;
1165 case 1: p.y = scx[0] * x[0]; break;
1166 case 2: p.z = scx[0] * x[0]; break;
1167 }
1168 if (ndim > 1) {
1169 switch (data->pdim[1]) {
1170 case 0: p.x = scx[1] * x[1]; break;
1171 case 1: p.y = scx[1] * x[1]; break;
1172 case 2: p.z = scx[1] * x[1]; break;
1173 }
1174 }
1175 }
1176
1177 if (data->is_ellipsoid && ndim > 0) {
1178 /* compute width of ellipsoid at this point, along the
1179 analytic-intersection direction */
1180 double dx = (x[0] - data->c[0]) * data->winv[0];
1181 double w = 1.0 - dx * dx;
1182 if (ndim > 1) { /* rescale 2nd dimension to stay inside ellipsoid */
1183 double x1;
1184 if (w < 0) return 0.0; /* outside the ellipsoid */
1185 scale_result = sqrt(w);
1186 x1 = data->c[1] + (x[1] - data->c[1]) * scale_result;
1187 switch (data->pdim[1]) {
1188 case 0: p.x = scx[1] * x1; break;
1189 case 1: p.y = scx[1] * x1; break;
1190 case 2: p.z = scx[1] * x1; break;
1191 }
1192 dx = (x1 - data->c[1]) * data->winv[1];
1193 w -= dx * dx;
1194 }
1195 if (w < 0) return 0.0; /* outside the ellipsoid */
1196 w = data->w0 * sqrt(w);
1197 a0 = data->c0 - w;
1198 b0 = data->c0 + w;
1199 }
1200
1201 return intersect_line_segment_with_object(p, data->dir, data->o, a0, b0) * scale_result;
1202 }
1203
overlap_with_object(geom_box b,int is_ellipsoid,geometric_object o,number tol,integer maxeval)1204 number overlap_with_object(geom_box b, int is_ellipsoid, geometric_object o, number tol,
1205 integer maxeval) {
1206 overlap_data data;
1207 int empty_x = b.low.x == b.high.x;
1208 int empty_y = b.low.y == b.high.y;
1209 int empty_z = b.low.z == b.high.z;
1210 double V0 = ((empty_x ? 1 : b.high.x - b.low.x) * (empty_y ? 1 : b.high.y - b.low.y) *
1211 (empty_z ? 1 : b.high.z - b.low.z));
1212 vector3 ex = {1, 0, 0}, ey = {0, 1, 0}, ez = {0, 0, 1};
1213 geom_box bb;
1214 double xmin[2] = {0, 0}, xmax[2] = {0, 0}, esterr;
1215 int errflag;
1216 unsigned i;
1217
1218 geom_get_bounding_box(o, &bb);
1219 if (!is_ellipsoid && !empty_x && !empty_y && !empty_z && /* todo: optimize 1d and 2d cases */
1220 bb.low.x >= b.low.x && bb.high.x <= b.high.x && bb.low.y >= b.low.y &&
1221 bb.high.y <= b.high.y && bb.low.z >= b.low.z && bb.high.z <= b.high.z)
1222 return geom_object_volume(o) /
1223 (V0 * fabs(matrix3x3_determinant(
1224 geometry_lattice.basis))); /* o is completely contained within b */
1225 geom_box_intersection(&bb, &b, &bb);
1226 if (bb.low.x > bb.high.x || bb.low.y > bb.high.y || bb.low.z > bb.high.z ||
1227 (!empty_x && bb.low.x == bb.high.x) || (!empty_y && bb.low.y == bb.high.y) ||
1228 (!empty_z && bb.low.z == bb.high.z))
1229 return 0.0;
1230
1231 data.winv[0] = data.winv[1] = data.w0 = 1.0;
1232 data.c[0] = data.c[1] = data.c0 = 0;
1233
1234 data.o = o;
1235 data.p.x = data.p.y = data.p.z = 0;
1236 data.dim = 0;
1237 if (!empty_x) {
1238 data.dir = ex;
1239 data.a0 = bb.low.x;
1240 data.b0 = bb.high.x;
1241 data.w0 = 0.5 * (b.high.x - b.low.x);
1242 data.c0 = 0.5 * (b.high.x + b.low.x);
1243 if (!empty_y) {
1244 xmin[data.dim] = bb.low.y;
1245 xmax[data.dim] = bb.high.y;
1246 data.winv[data.dim] = 2.0 / (b.high.y - b.low.y);
1247 data.c[data.dim] = 0.5 * (b.high.y + b.low.y);
1248 data.pdim[data.dim++] = 1;
1249 }
1250 if (!empty_z) {
1251 xmin[data.dim] = bb.low.z;
1252 xmax[data.dim] = bb.high.z;
1253 data.winv[data.dim] = 2.0 / (b.high.z - b.low.z);
1254 data.c[data.dim] = 0.5 * (b.high.z + b.low.z);
1255 data.pdim[data.dim++] = 2;
1256 }
1257 }
1258 else if (!empty_y) {
1259 data.dir = ey;
1260 data.a0 = bb.low.y;
1261 data.b0 = bb.high.y;
1262 data.w0 = 0.5 * (b.high.y - b.low.y);
1263 data.c0 = 0.5 * (b.high.y + b.low.y);
1264 if (!empty_x) {
1265 xmin[data.dim] = bb.low.x;
1266 xmax[data.dim] = bb.high.x;
1267 data.winv[data.dim] = 2.0 / (b.high.x - b.low.x);
1268 data.c[data.dim] = 0.5 * (b.high.x + b.low.x);
1269 data.pdim[data.dim++] = 0;
1270 }
1271 if (!empty_z) {
1272 xmin[data.dim] = bb.low.z;
1273 xmax[data.dim] = bb.high.z;
1274 data.winv[data.dim] = 2.0 / (b.high.z - b.low.z);
1275 data.c[data.dim] = 0.5 * (b.high.z + b.low.z);
1276 data.pdim[data.dim++] = 2;
1277 }
1278 }
1279 else if (!empty_z) {
1280 data.dir = ez;
1281 data.a0 = bb.low.z;
1282 data.b0 = bb.high.z;
1283 data.w0 = 0.5 * (b.high.z - b.low.z);
1284 data.c0 = 0.5 * (b.high.z + b.low.z);
1285 if (!empty_x) {
1286 xmin[data.dim] = bb.low.x;
1287 xmax[data.dim] = bb.high.x;
1288 data.winv[data.dim] = 2.0 / (b.high.x - b.low.x);
1289 data.c[data.dim] = 0.5 * (b.high.x + b.low.x);
1290 data.pdim[data.dim++] = 0;
1291 }
1292 if (!empty_y) {
1293 xmin[data.dim] = bb.low.y;
1294 xmax[data.dim] = bb.high.y;
1295 data.winv[data.dim] = 2.0 / (b.high.y - b.low.y);
1296 data.c[data.dim] = 0.5 * (b.high.y + b.low.y);
1297 data.pdim[data.dim++] = 1;
1298 }
1299 }
1300 else
1301 return 1.0;
1302
1303 #if 1
1304 /* To maintain mirror symmetries through the x/y/z axes, we flip
1305 the integration range whenever xmax < 0. (This is in case
1306 the integration routine is not fully symmetric, which may
1307 happen(?) due to the upper bound on the #evaluations.)*/
1308 for (i = 0; i < data.dim; ++i) {
1309 if (xmax[i] < 0) {
1310 double xm = xmin[i];
1311 data.scx[i] = -1;
1312 xmin[i] = -xmax[i];
1313 xmax[i] = -xm;
1314 data.c[i] = -data.c[i];
1315 }
1316 else
1317 data.scx[i] = 1;
1318 }
1319 #else
1320 for (i = 0; i < data.dim; ++i)
1321 data.scx[i] = 1;
1322 #endif
1323
1324 if ((data.is_ellipsoid = is_ellipsoid)) { /* data for ellipsoid calc. */
1325 if (data.dim == 1)
1326 V0 *= K_PI / 4;
1327 else if (data.dim == 2)
1328 V0 *= K_PI / 6;
1329 }
1330
1331 return adaptive_integration(overlap_integrand, xmin, xmax, data.dim, &data, 0.0, tol, maxeval,
1332 &esterr, &errflag) /
1333 V0;
1334 }
1335
box_overlap_with_object(geom_box b,geometric_object o,number tol,integer maxeval)1336 number box_overlap_with_object(geom_box b, geometric_object o, number tol, integer maxeval) {
1337 return overlap_with_object(b, 0, o, tol, maxeval);
1338 }
1339
ellipsoid_overlap_with_object(geom_box b,geometric_object o,number tol,integer maxeval)1340 number ellipsoid_overlap_with_object(geom_box b, geometric_object o, number tol, integer maxeval) {
1341 return overlap_with_object(b, 1, o, tol, maxeval);
1342 }
1343
range_overlap_with_object(vector3 low,vector3 high,geometric_object o,number tol,integer maxeval)1344 number CTLIO range_overlap_with_object(vector3 low, vector3 high, geometric_object o, number tol,
1345 integer maxeval) {
1346 geom_box b;
1347 b.low = low;
1348 b.high = high;
1349 return box_overlap_with_object(b, o, tol, maxeval);
1350 }
1351
1352 /**************************************************************************/
1353
1354 /* geom_box_tree: a tree of boxes and the objects contained within
1355 them. The tree recursively partitions the unit cell, allowing us
1356 to perform binary searches for the object containing a given point. */
1357
destroy_geom_box_tree(geom_box_tree t)1358 void destroy_geom_box_tree(geom_box_tree t) {
1359 if (t) {
1360 destroy_geom_box_tree(t->t1);
1361 destroy_geom_box_tree(t->t2);
1362 if (t->nobjects && t->objects) FREE(t->objects);
1363 FREE1(t);
1364 }
1365 }
1366
1367 /* return whether the object o, shifted by the vector shiftby,
1368 possibly intersects b. Upon return, obj_b is the bounding
1369 box for o. */
object_in_box(geometric_object o,vector3 shiftby,geom_box * obj_b,const geom_box * b)1370 static int object_in_box(geometric_object o, vector3 shiftby, geom_box *obj_b, const geom_box *b) {
1371 geom_get_bounding_box(o, obj_b);
1372 geom_box_shift(obj_b, shiftby);
1373 return geom_boxes_intersect(obj_b, b);
1374 }
1375
new_geom_box_tree(void)1376 static geom_box_tree new_geom_box_tree(void) {
1377 geom_box_tree t;
1378
1379 t = MALLOC1(struct geom_box_tree_struct);
1380 CHECK(t, "out of memory");
1381 t->t1 = t->t2 = NULL;
1382 t->nobjects = 0;
1383 t->objects = NULL;
1384 return t;
1385 }
1386
1387 /* Divide b into b1 and b2, cutting b in two along the axis
1388 divide_axis (0 = x, 1 = y, 2 = z) at divide_point. */
divide_geom_box(const geom_box * b,int divide_axis,number divide_point,geom_box * b1,geom_box * b2)1389 static void divide_geom_box(const geom_box *b, int divide_axis, number divide_point, geom_box *b1,
1390 geom_box *b2) {
1391 *b1 = *b2 = *b;
1392 switch (divide_axis) {
1393 case 0: b1->high.x = b2->low.x = divide_point; break;
1394 case 1: b1->high.y = b2->low.y = divide_point; break;
1395 case 2: b1->high.z = b2->low.z = divide_point; break;
1396 }
1397 }
1398
1399 #define VEC_I(v, i) ((i) == 0 ? (v).x : ((i) == 1 ? (v).y : (v).z))
1400 #define SMALL 1.0e-7
1401
1402 /* Find the best place (best_partition) to "cut" along the axis
1403 divide_axis in order to maximally divide the objects between
1404 the partitions. Upon return, n1 and n2 are the number of objects
1405 below and above the partition, respectively. */
find_best_partition(int nobjects,const geom_box_object * objects,int divide_axis,number * best_partition,int * n1,int * n2)1406 static void find_best_partition(int nobjects, const geom_box_object *objects, int divide_axis,
1407 number *best_partition, int *n1, int *n2) {
1408 number cur_partition;
1409 int i, j, cur_n1, cur_n2;
1410
1411 *n1 = *n2 = nobjects + 1;
1412 *best_partition = 0;
1413
1414 /* Search for the best partition, by checking all possible partitions
1415 either just above the high end of an object or just below the
1416 low end of an object. */
1417
1418 for (i = 0; i < nobjects; ++i) {
1419 cur_partition = VEC_I(objects[i].box.high, divide_axis) * (1 + SMALL);
1420 cur_n1 = cur_n2 = 0;
1421 for (j = 0; j < nobjects; ++j) {
1422 double low = VEC_I(objects[j].box.low, divide_axis);
1423 double high = VEC_I(objects[j].box.high, divide_axis);
1424 cur_n1 += low <= cur_partition;
1425 cur_n2 += high >= cur_partition;
1426 }
1427 CHECK(cur_n1 + cur_n2 >= nobjects, "assertion failure 1 in find_best_partition");
1428 if (MAX(cur_n1, cur_n2) < MAX(*n1, *n2)) {
1429 *best_partition = cur_partition;
1430 *n1 = cur_n1;
1431 *n2 = cur_n2;
1432 }
1433 }
1434 for (i = 0; i < nobjects; ++i) {
1435 cur_partition = VEC_I(objects[i].box.low, divide_axis) * (1 - SMALL);
1436 cur_n1 = cur_n2 = 0;
1437 for (j = 0; j < nobjects; ++j) {
1438 double low = VEC_I(objects[j].box.low, divide_axis);
1439 double high = VEC_I(objects[j].box.high, divide_axis);
1440 cur_n1 += low <= cur_partition;
1441 cur_n2 += high >= cur_partition;
1442 }
1443 CHECK(cur_n1 + cur_n2 >= nobjects, "assertion failure 2 in find_best_partition");
1444 if (MAX(cur_n1, cur_n2) < MAX(*n1, *n2)) {
1445 *best_partition = cur_partition;
1446 *n1 = cur_n1;
1447 *n2 = cur_n2;
1448 }
1449 }
1450 }
1451
1452 /* divide_geom_box_tree: recursively divide t in two, each time
1453 dividing along the axis that maximally partitions the boxes,
1454 and only stop partitioning when partitioning doesn't help any
1455 more. Upon return, t points to the partitioned tree. */
divide_geom_box_tree(geom_box_tree t)1456 static void divide_geom_box_tree(geom_box_tree t) {
1457 int division_nobjects[3][2] = {{0, 0}, {0, 0}, {0, 0}};
1458 number division_point[3];
1459 int best = -1;
1460 int i, j, n1, n2;
1461
1462 if (!t) return;
1463 if (t->t1 || t->t2) { /* this node has already been divided */
1464 divide_geom_box_tree(t->t1);
1465 divide_geom_box_tree(t->t2);
1466 return;
1467 }
1468
1469 if (t->nobjects <= 2) return; /* no point in partitioning */
1470
1471 /* Try partitioning along each dimension, counting the
1472 number of objects in the partitioned boxes and finding
1473 the best partition. */
1474 for (i = 0; i < dimensions; ++i) {
1475 if (VEC_I(t->b.high, i) == VEC_I(t->b.low, i)) continue; /* skip empty dimensions */
1476 find_best_partition(t->nobjects, t->objects, i, &division_point[i], &division_nobjects[i][0],
1477 &division_nobjects[i][1]);
1478 if (best < 0 || MAX(division_nobjects[i][0], division_nobjects[i][1]) <
1479 MAX(division_nobjects[best][0], division_nobjects[best][1]))
1480 best = i;
1481 }
1482
1483 /* don't do anything if division makes the worst case worse or if
1484 it fails to improve the best case: */
1485 if (best < 0 || MAX(division_nobjects[best][0], division_nobjects[best][1]) + 1 > t->nobjects ||
1486 MIN(division_nobjects[best][0], division_nobjects[best][1]) + 1 >= t->nobjects)
1487 return; /* division didn't help us */
1488
1489 divide_geom_box(&t->b, best, division_point[best], &t->b1, &t->b2);
1490 t->t1 = new_geom_box_tree();
1491 t->t2 = new_geom_box_tree();
1492 t->t1->b = t->b1;
1493 t->t2->b = t->b2;
1494
1495 t->t1->nobjects = division_nobjects[best][0];
1496 t->t1->objects = MALLOC(geom_box_object, t->t1->nobjects);
1497 CHECK(t->t1->objects, "out of memory");
1498
1499 t->t2->nobjects = division_nobjects[best][1];
1500 t->t2->objects = MALLOC(geom_box_object, t->t2->nobjects);
1501 CHECK(t->t2->objects, "out of memory");
1502
1503 for (j = n1 = n2 = 0; j < t->nobjects; ++j) {
1504 if (geom_boxes_intersect(&t->b1, &t->objects[j].box)) {
1505 CHECK(n1 < t->t1->nobjects, "BUG in divide_geom_box_tree");
1506 t->t1->objects[n1++] = t->objects[j];
1507 }
1508 if (geom_boxes_intersect(&t->b2, &t->objects[j].box)) {
1509 CHECK(n2 < t->t2->nobjects, "BUG in divide_geom_box_tree");
1510 t->t2->objects[n2++] = t->objects[j];
1511 }
1512 }
1513 CHECK(j == t->nobjects && n1 == t->t1->nobjects && n2 == t->t2->nobjects,
1514 "BUG in divide_geom_box_tree: wrong nobjects");
1515
1516 t->nobjects = 0;
1517 FREE(t->objects);
1518 t->objects = NULL;
1519
1520 divide_geom_box_tree(t->t1);
1521 divide_geom_box_tree(t->t2);
1522 }
1523
create_geom_box_tree(void)1524 geom_box_tree create_geom_box_tree(void) {
1525 geom_box b0;
1526 b0.low = vector3_plus(geometry_center, vector3_scale(-0.5, geometry_lattice.size));
1527 b0.high = vector3_plus(geometry_center, vector3_scale(0.5, geometry_lattice.size));
1528 return create_geom_box_tree0(geometry, b0);
1529 }
1530
num_objects_in_box(const geometric_object * o,vector3 shiftby,const geom_box * b)1531 static int num_objects_in_box(const geometric_object *o, vector3 shiftby, const geom_box *b) {
1532 if (o->which_subclass == GEOM COMPOUND_GEOMETRIC_OBJECT) {
1533 int n = o->subclass.compound_geometric_object_data->component_objects.num_items;
1534 geometric_object *os = o->subclass.compound_geometric_object_data->component_objects.items;
1535 int i, sum = 0;
1536 shiftby = vector3_plus(shiftby, o->center);
1537 for (i = 0; i < n; ++i)
1538 sum += num_objects_in_box(os + i, shiftby, b);
1539 return sum;
1540 }
1541 else {
1542 geom_box ob;
1543 return object_in_box(*o, shiftby, &ob, b);
1544 }
1545 }
1546
store_objects_in_box(const geometric_object * o,vector3 shiftby,const geom_box * b,geom_box_object * bo,int precedence)1547 static int store_objects_in_box(const geometric_object *o, vector3 shiftby, const geom_box *b,
1548 geom_box_object *bo, int precedence) {
1549 if (o->which_subclass == GEOM COMPOUND_GEOMETRIC_OBJECT) {
1550 int n = o->subclass.compound_geometric_object_data->component_objects.num_items;
1551 geometric_object *os = o->subclass.compound_geometric_object_data->component_objects.items;
1552 int i, sum = 0;
1553 shiftby = vector3_plus(shiftby, o->center);
1554 for (i = 0; i < n; ++i)
1555 sum += store_objects_in_box(os + i, shiftby, b, bo + sum, precedence - sum);
1556 return sum;
1557 }
1558 else {
1559 geom_box ob;
1560 if (object_in_box(*o, shiftby, &ob, b)) {
1561 bo->box = ob;
1562 bo->o = o;
1563 bo->shiftby = shiftby;
1564 bo->precedence = precedence;
1565 return 1;
1566 }
1567 else
1568 return 0;
1569 }
1570 }
1571
create_geom_box_tree0(geometric_object_list geometry,geom_box b0)1572 geom_box_tree create_geom_box_tree0(geometric_object_list geometry, geom_box b0) {
1573 geom_box_tree t = new_geom_box_tree();
1574 int i, index;
1575
1576 t->b = b0;
1577
1578 for (i = geometry.num_items - 1; i >= 0; --i) {
1579 vector3 shiftby = {0, 0, 0};
1580 if (ensure_periodicity) {
1581 LOOP_PERIODIC(shiftby, t->nobjects += num_objects_in_box(geometry.items + i, shiftby, &t->b));
1582 }
1583 else
1584 t->nobjects += num_objects_in_box(geometry.items + i, shiftby, &t->b);
1585 }
1586
1587 t->objects = MALLOC(geom_box_object, t->nobjects);
1588 CHECK(t->objects || t->nobjects == 0, "out of memory");
1589
1590 for (i = geometry.num_items - 1, index = 0; i >= 0; --i) {
1591 vector3 shiftby = {0, 0, 0};
1592 if (ensure_periodicity) {
1593 int precedence = t->nobjects - index;
1594 LOOP_PERIODIC(shiftby, index += store_objects_in_box(geometry.items + i, shiftby, &t->b,
1595 t->objects + index, precedence));
1596 }
1597 else
1598 index += store_objects_in_box(geometry.items + i, shiftby, &t->b, t->objects + index,
1599 t->nobjects - index);
1600 }
1601 CHECK(index == t->nobjects, "bug in create_geom_box_tree0");
1602
1603 divide_geom_box_tree(t);
1604
1605 return t;
1606 }
1607
1608 /* create a new tree from t, pruning all nodes that don't intersect b */
restrict_geom_box_tree(geom_box_tree t,const geom_box * b)1609 geom_box_tree restrict_geom_box_tree(geom_box_tree t, const geom_box *b) {
1610 geom_box_tree tr;
1611 int i, j;
1612
1613 if (!t || !geom_boxes_intersect(&t->b, b)) return NULL;
1614
1615 tr = new_geom_box_tree();
1616
1617 for (i = 0, j = 0; i < t->nobjects; ++i)
1618 if (geom_boxes_intersect(&t->objects[i].box, b)) ++j;
1619 tr->nobjects = j;
1620 tr->objects = MALLOC(geom_box_object, tr->nobjects);
1621 CHECK(tr->objects || tr->nobjects == 0, "out of memory");
1622
1623 for (i = 0, j = 0; i < t->nobjects; ++i)
1624 if (geom_boxes_intersect(&t->objects[i].box, b)) tr->objects[j++] = t->objects[i];
1625
1626 tr->t1 = restrict_geom_box_tree(t->t1, b);
1627 tr->t2 = restrict_geom_box_tree(t->t2, b);
1628
1629 if (tr->nobjects == 0) {
1630 if (tr->t1 && !tr->t2) {
1631 geom_box_tree tr0 = tr;
1632 tr = tr->t1;
1633 FREE1(tr0);
1634 }
1635 else if (tr->t2 && !tr->t1) {
1636 geom_box_tree tr0 = tr;
1637 tr = tr->t2;
1638 FREE1(tr0);
1639 }
1640 }
1641
1642 return tr;
1643 }
1644
1645 /**************************************************************************/
1646
1647 /* recursively search the tree for the given point, returning the
1648 subtree (if any) that contains it and the index oindex of the
1649 object in that tree. The input value of oindex indicates the
1650 starting object to search in t (0 to search all). */
tree_search(vector3 p,geom_box_tree t,int * oindex)1651 static geom_box_tree tree_search(vector3 p, geom_box_tree t, int *oindex) {
1652 int i;
1653 geom_box_tree gbt;
1654
1655 if (!t || !geom_box_contains_point(&t->b, p)) return NULL;
1656
1657 for (i = *oindex; i < t->nobjects; ++i)
1658 if (geom_box_contains_point(&t->objects[i].box, p) &&
1659 point_in_fixed_objectp(vector3_minus(p, t->objects[i].shiftby), *t->objects[i].o)) {
1660 *oindex = i;
1661 return t;
1662 }
1663
1664 *oindex = 0;
1665 gbt = tree_search(p, t->t1, oindex);
1666 if (!gbt) gbt = tree_search(p, t->t2, oindex);
1667 return gbt;
1668 }
1669
1670 /* shift p to be within the unit cell of the lattice (centered on the
1671 origin) */
shift_to_unit_cell(vector3 p)1672 vector3 shift_to_unit_cell(vector3 p) {
1673 while (p.x >= 0.5 * geometry_lattice.size.x)
1674 p.x -= geometry_lattice.size.x;
1675 while (p.x < -0.5 * geometry_lattice.size.x)
1676 p.x += geometry_lattice.size.x;
1677 while (p.y >= 0.5 * geometry_lattice.size.y)
1678 p.y -= geometry_lattice.size.y;
1679 while (p.y < -0.5 * geometry_lattice.size.y)
1680 p.y += geometry_lattice.size.y;
1681 while (p.z >= 0.5 * geometry_lattice.size.z)
1682 p.z -= geometry_lattice.size.z;
1683 while (p.z < -0.5 * geometry_lattice.size.z)
1684 p.z += geometry_lattice.size.z;
1685 return p;
1686 }
1687
object_of_point_in_tree(vector3 p,geom_box_tree t,vector3 * shiftby,int * precedence)1688 const geometric_object *object_of_point_in_tree(vector3 p, geom_box_tree t, vector3 *shiftby,
1689 int *precedence) {
1690 int oindex = 0;
1691 t = tree_search(p, t, &oindex);
1692 if (t) {
1693 geom_box_object *gbo = t->objects + oindex;
1694 *shiftby = gbo->shiftby;
1695 *precedence = gbo->precedence;
1696 return gbo->o;
1697 }
1698 else {
1699 shiftby->x = shiftby->y = shiftby->z = 0;
1700 *precedence = 0;
1701 return 0;
1702 }
1703 }
1704
material_of_unshifted_point_in_tree_inobject(vector3 p,geom_box_tree t,boolean * inobject)1705 material_type material_of_unshifted_point_in_tree_inobject(vector3 p, geom_box_tree t,
1706 boolean *inobject) {
1707 int oindex = 0;
1708 t = tree_search(p, t, &oindex);
1709 if (t) {
1710 *inobject = 1;
1711 return (t->objects[oindex].o->material);
1712 }
1713 else {
1714 *inobject = 0;
1715 return default_material;
1716 }
1717 }
1718
material_of_point_in_tree_inobject(vector3 p,geom_box_tree t,boolean * inobject)1719 material_type material_of_point_in_tree_inobject(vector3 p, geom_box_tree t, boolean *inobject) {
1720 /* backwards compatibility */
1721 return material_of_unshifted_point_in_tree_inobject(shift_to_unit_cell(p), t, inobject);
1722 }
1723
material_of_point_in_tree(vector3 p,geom_box_tree t)1724 material_type material_of_point_in_tree(vector3 p, geom_box_tree t) {
1725 boolean inobject;
1726 return material_of_point_in_tree_inobject(p, t, &inobject);
1727 }
1728
geom_tree_search_next(vector3 p,geom_box_tree t,int * oindex)1729 geom_box_tree geom_tree_search_next(vector3 p, geom_box_tree t, int *oindex) {
1730 *oindex += 1; /* search starting at next oindex */
1731 return tree_search(p, t, oindex);
1732 }
1733
geom_tree_search(vector3 p,geom_box_tree t,int * oindex)1734 geom_box_tree geom_tree_search(vector3 p, geom_box_tree t, int *oindex) {
1735 *oindex = -1; /* search all indices > -1 */
1736 return geom_tree_search_next(p, t, oindex);
1737 }
1738
1739 /**************************************************************************/
1740
1741 /* convert a vector p in the given object to some coordinate
1742 in [0,1]^3 that is a more "natural" map of the object interior. */
to_geom_box_coords(vector3 p,geom_box_object * gbo)1743 vector3 to_geom_box_coords(vector3 p, geom_box_object *gbo) {
1744 return to_geom_object_coords(vector3_minus(p, gbo->shiftby), *gbo->o);
1745 }
1746
1747 /**************************************************************************/
1748
display_geom_box_tree(int indentby,geom_box_tree t)1749 void display_geom_box_tree(int indentby, geom_box_tree t) {
1750 int i;
1751
1752 if (!t) return;
1753 ctl_printf("%*sbox (%g..%g, %g..%g, %g..%g)\n", indentby, "", t->b.low.x, t->b.high.x, t->b.low.y,
1754 t->b.high.y, t->b.low.z, t->b.high.z);
1755 for (i = 0; i < t->nobjects; ++i) {
1756 ctl_printf("%*sbounding box (%g..%g, %g..%g, %g..%g)\n", indentby + 5, "",
1757 t->objects[i].box.low.x, t->objects[i].box.high.x, t->objects[i].box.low.y,
1758 t->objects[i].box.high.y, t->objects[i].box.low.z, t->objects[i].box.high.z);
1759 ctl_printf("%*sshift object by (%g, %g, %g)\n", indentby + 5, "", t->objects[i].shiftby.x,
1760 t->objects[i].shiftby.y, t->objects[i].shiftby.z);
1761 display_geometric_object_info(indentby + 5, *t->objects[i].o);
1762 }
1763 display_geom_box_tree(indentby + 5, t->t1);
1764 display_geom_box_tree(indentby + 5, t->t2);
1765 }
1766
1767 /**************************************************************************/
1768
1769 /* Computing tree statistics (depth and number of nodes): */
1770
1771 /* helper function for geom_box_tree_stats */
get_tree_stats(geom_box_tree t,int * depth,int * nobjects)1772 static void get_tree_stats(geom_box_tree t, int *depth, int *nobjects) {
1773 if (t) {
1774 int d1, d2;
1775
1776 *nobjects += t->nobjects;
1777 d1 = d2 = *depth + 1;
1778 get_tree_stats(t->t1, &d1, nobjects);
1779 get_tree_stats(t->t2, &d2, nobjects);
1780 *depth = MAX(d1, d2);
1781 }
1782 }
1783
geom_box_tree_stats(geom_box_tree t,int * depth,int * nobjects)1784 void geom_box_tree_stats(geom_box_tree t, int *depth, int *nobjects) {
1785 *depth = *nobjects = 0;
1786 get_tree_stats(t, depth, nobjects);
1787 }
1788
1789 /**************************************************************************/
1790
1791 #ifndef LIBCTLGEOM
1792
get_grid_size(void)1793 vector3 get_grid_size(void) {
1794 return ctl_convert_vector3_to_c(gh_call0(gh_lookup("get-grid-size")));
1795 }
1796
get_resolution(void)1797 vector3 get_resolution(void) {
1798 return ctl_convert_vector3_to_c(gh_call0(gh_lookup("get-resolution")));
1799 }
1800
get_grid_size_n(int * nx,int * ny,int * nz)1801 void get_grid_size_n(int *nx, int *ny, int *nz) {
1802 vector3 grid_size;
1803 grid_size = get_grid_size();
1804 *nx = (int)grid_size.x;
1805 *ny = (int)grid_size.y;
1806 *nz = (int)grid_size.z;
1807 }
1808
1809 #endif
1810
1811 /**************************************************************************/
1812
1813 /* constructors for the geometry types (ugh, wish these
1814 could be automatically generated from geom.scm) */
1815
make_geometric_object(material_type material,vector3 center)1816 geometric_object make_geometric_object(material_type material, vector3 center) {
1817 geometric_object o;
1818 material_type_copy(&material, &o.material);
1819 o.center = center;
1820 o.which_subclass = GEOM GEOMETRIC_OBJECT_SELF;
1821 return o;
1822 }
1823
make_cylinder(material_type material,vector3 center,number radius,number height,vector3 axis)1824 geometric_object make_cylinder(material_type material, vector3 center, number radius, number height,
1825 vector3 axis) {
1826 geometric_object o = make_geometric_object(material, center);
1827 o.which_subclass = GEOM CYLINDER;
1828 o.subclass.cylinder_data = MALLOC1(cylinder);
1829 CHECK(o.subclass.cylinder_data, "out of memory");
1830 o.subclass.cylinder_data->radius = radius;
1831 o.subclass.cylinder_data->height = height;
1832 o.subclass.cylinder_data->axis = axis;
1833 o.subclass.cylinder_data->which_subclass = CYL CYLINDER_SELF;
1834 geom_fix_object_ptr(&o);
1835 return o;
1836 }
1837
make_cone(material_type material,vector3 center,number radius,number height,vector3 axis,number radius2)1838 geometric_object make_cone(material_type material, vector3 center, number radius, number height,
1839 vector3 axis, number radius2) {
1840 geometric_object o = make_cylinder(material, center, radius, height, axis);
1841 o.subclass.cylinder_data->which_subclass = CYL CONE;
1842 o.subclass.cylinder_data->subclass.cone_data = MALLOC1(cone);
1843 CHECK(o.subclass.cylinder_data->subclass.cone_data, "out of memory");
1844 o.subclass.cylinder_data->subclass.cone_data->radius2 = radius2;
1845 return o;
1846 }
1847
make_wedge(material_type material,vector3 center,number radius,number height,vector3 axis,number wedge_angle,vector3 wedge_start)1848 geometric_object make_wedge(material_type material, vector3 center, number radius, number height,
1849 vector3 axis, number wedge_angle, vector3 wedge_start) {
1850 geometric_object o = make_cylinder(material, center, radius, height, axis);
1851 o.subclass.cylinder_data->which_subclass = CYL WEDGE;
1852 o.subclass.cylinder_data->subclass.wedge_data = MALLOC1(wedge);
1853 CHECK(o.subclass.cylinder_data->subclass.wedge_data, "out of memory");
1854 o.subclass.cylinder_data->subclass.wedge_data->wedge_angle = wedge_angle;
1855 o.subclass.cylinder_data->subclass.wedge_data->wedge_start = wedge_start;
1856 geom_fix_object_ptr(&o);
1857 return o;
1858 }
1859
make_sphere(material_type material,vector3 center,number radius)1860 geometric_object make_sphere(material_type material, vector3 center, number radius) {
1861 geometric_object o = make_geometric_object(material, center);
1862 o.which_subclass = GEOM SPHERE;
1863 o.subclass.sphere_data = MALLOC1(sphere);
1864 CHECK(o.subclass.sphere_data, "out of memory");
1865 o.subclass.sphere_data->radius = radius;
1866 return o;
1867 }
1868
make_block(material_type material,vector3 center,vector3 e1,vector3 e2,vector3 e3,vector3 size)1869 geometric_object make_block(material_type material, vector3 center, vector3 e1, vector3 e2,
1870 vector3 e3, vector3 size) {
1871 geometric_object o = make_geometric_object(material, center);
1872 o.which_subclass = GEOM BLOCK;
1873 o.subclass.block_data = MALLOC1(block);
1874 CHECK(o.subclass.block_data, "out of memory");
1875 o.subclass.block_data->e1 = e1;
1876 o.subclass.block_data->e2 = e2;
1877 o.subclass.block_data->e3 = e3;
1878 o.subclass.block_data->size = size;
1879 o.subclass.block_data->which_subclass = BLK BLOCK_SELF;
1880 geom_fix_object_ptr(&o);
1881 return o;
1882 }
1883
make_ellipsoid(material_type material,vector3 center,vector3 e1,vector3 e2,vector3 e3,vector3 size)1884 geometric_object make_ellipsoid(material_type material, vector3 center, vector3 e1, vector3 e2,
1885 vector3 e3, vector3 size) {
1886 geometric_object o = make_block(material, center, e1, e2, e3, size);
1887 o.subclass.block_data->which_subclass = BLK ELLIPSOID;
1888 o.subclass.block_data->subclass.ellipsoid_data = MALLOC1(ellipsoid);
1889 CHECK(o.subclass.block_data->subclass.ellipsoid_data, "out of memory");
1890 o.subclass.block_data->subclass.ellipsoid_data->inverse_semi_axes.x = 2.0 / size.x;
1891 o.subclass.block_data->subclass.ellipsoid_data->inverse_semi_axes.y = 2.0 / size.y;
1892 o.subclass.block_data->subclass.ellipsoid_data->inverse_semi_axes.z = 2.0 / size.z;
1893 return o;
1894 }
1895
1896 /***************************************************************
1897 * The remainder of this file implements geometric primitives for prisms.
1898 * A prism is a planar polygon, consisting of 3 or more user-specified
1899 * vertices (the "bottom_vertices), extruded through a given thickness
1900 * (the "height") in the direction of a given unit vector (the "axis")
1901 * with the walls of the extrusion tapering at a given angle angle
1902 * (the "sidewall_angle).
1903 * Most calculations are done in the "prism coordinate system",
1904 * in which the prism floor lies in the XY plane with centroid
1905 * at the origin and the prism axis is the positive Z-axis.
1906 * Some variable naming conventions:
1907 * -- Suffix 'p' or '_p' on variable names identifies variables
1908 * storing coordinates or vector components in the prism system.
1909 * Suffix 'c' or '_c' (or no suffix) corresponds to coodinates/components
1910 * in ordinary 3d space. ('c' stands for 'cartesian').
1911 * -- We use the term 'vertex' for points in 3-space, stored as vector3
1912 * quantities with variable names beginning with 'p' or 'v'. For 3D
1913 * direction vectors we use variable names beginning with 'd'.
1914 * -- We use the term 'node' for points in 2-space, stored as vector3
1915 * quantities (with the z component unused) with variables beginning with 'q'.
1916 * For 2D direction vectors we use variable names beginning with 'u'.
1917 * homer reid 4/2018
1918 ***************************************************************/
1919
1920 /***************************************************************/
1921 /* given coordinates of a point in the prism coordinate system,*/
1922 /* return cartesian coordinates of that point */
1923 /***************************************************************/
prism_coordinate_p2c(prism * prsm,vector3 pp)1924 vector3 prism_coordinate_p2c(prism *prsm, vector3 pp) {
1925 return vector3_plus(prsm->centroid, matrix3x3_vector3_mult(prsm->m_p2c, pp));
1926 }
1927
prism_vector_p2c(prism * prsm,vector3 vp)1928 vector3 prism_vector_p2c(prism *prsm, vector3 vp) {
1929 return matrix3x3_vector3_mult(prsm->m_p2c, vp);
1930 }
1931
prism_coordinate_c2p(prism * prsm,vector3 pc)1932 vector3 prism_coordinate_c2p(prism *prsm, vector3 pc) {
1933 return matrix3x3_vector3_mult(prsm->m_c2p, vector3_minus(pc, prsm->centroid));
1934 }
1935
prism_vector_c2p(prism * prsm,vector3 vc)1936 vector3 prism_vector_c2p(prism *prsm, vector3 vc) {
1937 return matrix3x3_vector3_mult(prsm->m_c2p, vc);
1938 }
1939
1940 /***************************************************************/
1941 /* given 2D points q0,q1,q2 and a 2D vector u, determine */
1942 /* whether or not the line q0 + s*u intersects the line */
1943 /* segment q1--q2. */
1944 /* algorithm: solve the 2x2 linear system q0+s*u = q1+t*(q2-q1)*/
1945 /* for the scalar quantities s, t; intersection corresponds to */
1946 /* 0 <= t < 1. */
1947 /* return values: */
1948 /* ** case 1: u is not parallel to q1--q2 ** */
1949 /* NON_INTERSECTING: test negative */
1950 /* INTERSECTING: test positive */
1951 /* ** case 2: u is parallel to q1--q2 ** */
1952 /* IN_SEGMENT: q0 lies on line segment q1--q2 */
1953 /* ON_RAY: q0 does not lie on q1--q2, but there is a*/
1954 /* *positive* value of s such that q0+s*u */
1955 /* lies on q1--q2 */
1956 /* NON_INTERSECTING neither of the above */
1957 /***************************************************************/
1958 #define THRESH 1.0e-5
1959 #define NON_INTERSECTING 0
1960 #define INTERSECTING 1
1961 #define IN_SEGMENT 2
1962 #define ON_RAY 3
intersect_line_with_segment(vector3 q0,vector3 q1,vector3 q2,vector3 u,double * s)1963 int intersect_line_with_segment(vector3 q0, vector3 q1, vector3 q2, vector3 u, double *s) {
1964 /* ||ux q1x-q2x|| |s| = | q1x-q0x | */
1965 /* ||uy q1y-q2y|| |t| = | q1y-q0y | */
1966 double M00 = u.x, M01 = q1.x - q2.x;
1967 double M10 = u.y, M11 = q1.y - q2.y;
1968 double RHSx = q1.x - q0.x;
1969 double RHSy = q1.y - q0.y;
1970 double DetM = M00 * M11 - M01 * M10;
1971 double L2 = M01 * M01 + M11 * M11; // squared length of edge, used to set length scale
1972 if (fabs(DetM) < 1.0e-10 * L2) { // d zero or nearly parallel to edge
1973 if (vector3_nearly_equal(q0, q1, 1e-12) || vector3_nearly_equal(q0, q2, 1e-12)) return IN_SEGMENT;
1974 double q01x = q0.x - q1.x, q01y = q0.y - q1.y, q01 = sqrt(q01x * q01x + q01y * q01y);
1975 double q02x = q0.x - q2.x, q02y = q0.y - q2.y, q02 = sqrt(q02x * q02x + q02y * q02y);
1976 double dot = q01x * q02x + q01y * q02y;
1977 if (fabs(dot) < (1.0 - THRESH) * q01 * q02)
1978 return NON_INTERSECTING;
1979 else if (dot < 0.0) {
1980 *s = 0.0;
1981 return IN_SEGMENT;
1982 }
1983 else if ((u.x * q01x + u.y * q01y) < 0.0) {
1984 *s = fmin(q01, q02) / sqrt(u.x * u.x + u.y * u.y);
1985 return ON_RAY;
1986 }
1987 return NON_INTERSECTING;
1988 }
1989
1990 float t = (M00 * RHSy - M10 * RHSx) / DetM;
1991 if (s) *s = (M11 * RHSx - M01 * RHSy) / DetM;
1992
1993 // the plumb line intersects the segment if 0<=t<=1, with t==0,1
1994 // corresponding to the endpoints; for our purposes we count
1995 // the intersection if the plumb line runs through the t==0 vertex, but
1996 // NOT the t==1 vertex, to avoid double-counting for complete polygons.
1997 return (t < -THRESH || t >= (1 - THRESH)) ? NON_INTERSECTING : INTERSECTING;
1998 }
1999
2000 // like the previous routine, but only count intersections if s>=0
intersect_ray_with_segment(vector3 q0,vector3 q1,vector3 q2,vector3 u,double * s)2001 boolean intersect_ray_with_segment(vector3 q0, vector3 q1, vector3 q2, vector3 u, double *s) {
2002 double ss;
2003 int status = intersect_line_with_segment(q0, q1, q2, u, &ss);
2004 if (status == INTERSECTING && ss < 0.0) return NON_INTERSECTING;
2005 if (s) *s = ss;
2006 return status;
2007 }
2008
2009 /***************************************************************/
2010 /* 2D point-in-polygon test: return 1 if q0 lies within the */
2011 /* polygon with the given vertices, 0 otherwise. */
2012 // method: cast a plumb line in the positive x direction from */
2013 /* q0 to infinity and count the number of edges intersected; */
2014 /* point lies in polygon iff this is number is odd. */
2015 /***************************************************************/
2016 /* Implementation of: */
2017 /* */
2018 /* M. Galetzka and P. Glauner, "A Simple and Correct Even-Odd */
2019 /* Algorithm for the Point-in-Polygon Problem for Complex */
2020 /* Polygons", Proceedings of the 12th International Joint */
2021 /* Conference on Computer Vision, Imaging and Computer */
2022 /* Graphics Theory and Applications (VISIGRAPP 2017), Volume */
2023 /* 1: GRAPP, Porto, Portugal, 2017. */
2024 /***************************************************************/
2025
node_in_or_on_polygon(vector3 q0,vector3 * nodes,int num_nodes,boolean include_boundaries)2026 boolean node_in_or_on_polygon(vector3 q0, vector3 *nodes, int num_nodes,
2027 boolean include_boundaries) {
2028 // Create axes
2029 vector3 xAxis = {1.0, 0.0, 0.0};
2030
2031 // Initial start point
2032 vector3 startPoint;
2033 vector3 endPoint;
2034
2035 int startNodePosition = -1;
2036 int nn, edges_crossed = 0;
2037
2038 // Is q0 on a vertex or edge?
2039 // Transform coordinate system of nodes such that q0 is at 0|0
2040 for (nn = 0; nn < num_nodes; nn++) {
2041 int status = intersect_ray_with_segment(q0, nodes[nn], nodes[(nn + 1) % num_nodes], unit_vector3(vector3_minus(nodes[(nn + 1) % num_nodes], nodes[nn])), 0);
2042 if (status == IN_SEGMENT) { return include_boundaries; }
2043
2044 // Find start point which is not on the x axis (from q0)
2045 if (fabs(nodes[nn].y - q0.y) > THRESH) {
2046 startNodePosition = nn;
2047 startPoint = nodes[startNodePosition];
2048 }
2049 }
2050
2051 // No start point found and point is not on an edge or node
2052 // --> point is outside
2053 if (startNodePosition == -1) { return 0; }
2054
2055 int checkedPoints = 0;
2056 nn = startNodePosition;
2057
2058 // Consider all edges
2059 while (checkedPoints < num_nodes) {
2060 int savedIndex = (nn + 1) % num_nodes;
2061 double savedX = nodes[savedIndex].x;
2062
2063 // Move to next point which is not on the x-axis
2064 do {
2065 nn = (nn + 1) % num_nodes;
2066 checkedPoints++;
2067 } while (fabs(nodes[nn].y - q0.y) < THRESH);
2068 // Found end point
2069 endPoint = nodes[nn];
2070
2071 // Only intersect lines that cross the x-axis (don't need to correct for rounding
2072 // error in the if statement because startPoint and endPoint are screened to
2073 // never lie on the x-axis)
2074 if ((startPoint.y - q0.y) * (endPoint.y - q0.y) < 0) {
2075 // No nodes have been skipped and the successor node
2076 // has been chose as the end point
2077 if (savedIndex == nn) {
2078 int status = intersect_ray_with_segment(q0, startPoint, endPoint, xAxis, 0);
2079 if (status == INTERSECTING) { edges_crossed++; }
2080 }
2081 // If at least one node on the right side has been skipped,
2082 // the original edge would have been intersected
2083 // --> intersect with full x-axis
2084 else if (savedX > q0.x + THRESH) {
2085 int status = intersect_line_with_segment(q0, startPoint, endPoint, xAxis, 0);
2086 if (status == INTERSECTING) { edges_crossed++; }
2087 }
2088 }
2089 // End point is the next start point
2090 startPoint = endPoint;
2091 }
2092
2093 // Odd count --> in the polygon (1)
2094 // Even count --> outside (0)
2095 return edges_crossed % 2;
2096 }
2097
node_in_polygon(double q0x,double q0y,vector3 * nodes,int num_nodes)2098 boolean node_in_polygon(double q0x, double q0y, vector3 *nodes, int num_nodes) {
2099 vector3 q0;
2100 q0.x = q0x;
2101 q0.y = q0y;
2102 q0.z = 0.0;
2103 return node_in_or_on_polygon(q0, nodes, num_nodes, 1);
2104 }
2105
2106 /***************************************************************/
2107 /* return 1 or 0 if pc lies inside or outside the prism */
2108 /***************************************************************/
point_in_or_on_prism(prism * prsm,vector3 pc,boolean include_boundaries)2109 boolean point_in_or_on_prism(prism *prsm, vector3 pc, boolean include_boundaries) {
2110 double height = prsm->height;
2111 vector3 pp = prism_coordinate_c2p(prsm, pc);
2112 if (pp.z < 0.0 || pp.z > prsm->height) return 0;
2113 int num_nodes = prsm->vertices_p.num_items;
2114 vector3 nodes[num_nodes];
2115 int nv;
2116 for (nv = 0; nv < num_nodes; nv++) {
2117 nodes[nv] = vector3_plus(prsm->vertices_p.items[nv], vector3_scale(pp.z, prsm->top_polygon_diff_vectors_scaled_p.items[nv]));
2118 }
2119 return node_in_or_on_polygon(pp, nodes, num_nodes, include_boundaries);
2120 }
2121
point_in_prism(prism * prsm,vector3 pc)2122 boolean point_in_prism(prism *prsm, vector3 pc) {
2123 // by default, points on polygon edges are considered to lie inside the
2124 // polygon; this can be reversed by setting the environment variable
2125 // LIBCTL_EXCLUDE_BOUNDARIES=1
2126 static boolean include_boundaries = 1, init = 0;
2127 if (init == 0) {
2128 init = 1;
2129 char *s = getenv("LIBCTL_EXCLUDE_BOUNDARIES");
2130 if (s && s[0] == '1') include_boundaries = 0;
2131 }
2132 return point_in_or_on_prism(prsm, pc, include_boundaries);
2133 }
2134
2135 // comparator for qsort
dcmp(const void * pd1,const void * pd2)2136 static int dcmp(const void *pd1, const void *pd2) {
2137 double d1 = *((double *)pd1), d2 = *((double *)pd2);
2138 return (d1 < d2) ? -1.0 : (d1 > d2) ? 1.0 : 0.0;
2139 }
2140
2141 /******************************************************************/
2142 /* 3D line-prism intersection: compute all values of s at which */
2143 /* the line p+s*d intersects a prism face. */
2144 /* pc, dc = cartesian coordinates of p, cartesian components of d */
2145 /* slist is a caller-allocated buffer with enough room for */
2146 /* at least num_vertices+2 doubles. on return it contains */
2147 /* the intersection s-values sorted in ascending order. */
2148 /* the return value is the number of intersections. */
2149 /******************************************************************/
intersect_line_with_prism(prism * prsm,vector3 pc,vector3 dc,double * slist)2150 int intersect_line_with_prism(prism *prsm, vector3 pc, vector3 dc, double *slist) {
2151 vector3 pp = prism_coordinate_c2p(prsm, pc);
2152 vector3 dp = prism_vector_c2p(prsm, dc);
2153 vector3 *vps_bottom = prsm->vertices_p.items;
2154 vector3 *vps_top = prsm->vertices_top_p.items;
2155 int num_vertices = prsm->vertices_p.num_items;
2156 double height = prsm->height;
2157
2158 // identify intersections with prism side faces
2159 double tus_tolerance = 1e-8;
2160 int num_intersections = 0;
2161 int nv;
2162 for (nv = 0; nv < num_vertices; nv++) {
2163 int nvp1 = nv + 1;
2164 if (nvp1 == num_vertices) nvp1 = 0;
2165
2166 // checks if dp is parallel to the plane of the prism side face under consideration
2167 vector3 v1 = vector3_minus(vps_bottom[nvp1], vps_bottom[nv]);
2168 vector3 v2 = vector3_minus(vps_top[nv], vps_bottom[nv]);
2169 double dot_tolerance = 1e-6;
2170 if (fabs(vector3_dot(dp, vector3_cross(v1, v2))) <= dot_tolerance) continue;
2171
2172 // to find the intersection point pp + s*dp between the line and the
2173 // prism side face, we will solve the vector equation
2174 // pp + s*dp = o + t*v1 + u*v2
2175 // where o is vps_bottom[nv], v1 is vps_bottom[nvp1]-vps_bottom[nv],
2176 // v2 is vps_top[nv]-vps_bottom[nv], and 0 <= t <= 1, 0 <= u <= 1.
2177 matrix3x3 M;
2178 M.c0 = v1;
2179 M.c1 = v2;
2180 M.c2 = vector3_scale(-1, dp);
2181 vector3 RHS = vector3_minus(pp, vps_bottom[nv]);
2182 vector3 tus = matrix3x3_vector3_mult(matrix3x3_inverse(M), RHS);
2183 if (tus.x < -tus_tolerance || tus.x > 1+tus_tolerance || tus.y < -tus_tolerance || tus.y > 1+tus_tolerance) continue;
2184 double s = tus.z;
2185 slist[num_intersections++] = s;
2186 }
2187
2188 // identify intersections with prism ceiling and floor faces
2189 int LowerUpper;
2190 if (fabs(dp.z) > 1.0e-7 * vector3_norm(dp))
2191 for (LowerUpper = 0; LowerUpper < 2; LowerUpper++) {
2192 double z0p = LowerUpper ? height : 0.0;
2193 double s = (z0p - pp.z) / dp.z;
2194 vector3 *vps = LowerUpper ? vps_top : vps_bottom;
2195 if (!node_in_polygon(pp.x + s * dp.x, pp.y + s * dp.y, vps, num_vertices)) continue;
2196 slist[num_intersections++] = s;
2197 }
2198
2199 qsort((void *)slist, num_intersections, sizeof(double), dcmp);
2200 // if num_intersections is zero then just return that
2201 if (num_intersections == 0) return num_intersections;
2202 else {
2203 // remove duplicates from slist
2204 double duplicate_tolerance = 1e-3;
2205 int num_unique_elements = 1;
2206 double slist_unique[num_vertices+2];
2207 slist_unique[0] = slist[0];
2208 for (nv = 1; nv < num_intersections; nv++) {
2209 if (fabs(slist[nv] - slist[nv-1]) > duplicate_tolerance*fabs(slist[nv])) {
2210 slist_unique[num_unique_elements] = slist[nv];
2211 num_unique_elements++;
2212 }
2213 }
2214 slist = slist_unique;
2215 num_intersections = num_unique_elements;
2216 return num_intersections;
2217 }
2218 }
2219
2220 /***************************************************************/
2221 /***************************************************************/
2222 /***************************************************************/
intersect_line_segment_with_prism(prism * prsm,vector3 pc,vector3 dc,double a,double b)2223 double intersect_line_segment_with_prism(prism *prsm, vector3 pc, vector3 dc, double a, double b) {
2224 double *slist = prsm->workspace.items;
2225 int num_intersections = intersect_line_with_prism(prsm, pc, dc, slist);
2226
2227 // na=smallest index such that slist[na] > a
2228 int na = -1;
2229 int ns;
2230 for (ns = 0; na == -1 && ns < num_intersections; ns++)
2231 if (slist[ns] > a) na = ns;
2232
2233 if (na == -1) return 0.0;
2234
2235 int inside = ((na % 2) == 0 ? 0 : 1);
2236 double last_s = a;
2237 double ds = 0.0;
2238 for (ns = na; ns < num_intersections; ns++) {
2239 double this_s = fmin(b, slist[ns]);
2240 if (inside) ds += (this_s - last_s);
2241 if (b < slist[ns]) break;
2242 inside = (1 - inside);
2243 last_s = this_s;
2244 }
2245 return ds > 0.0 ? ds : 0.0;
2246 }
2247
2248 /***************************************************************/
2249 /* compute the minimum distance from a 3D point p to the */
2250 /* line segment with endpoints v1,v2. */
2251 /* algorithm: let pLine = v1 + d*(v2-v1) be the point on the */
2252 /* line closest to p; d is defined by minimizing |p-pLine|^2. */
2253 /* --> |p-v1|^2 + d^2 |v2-v1|^2 - 2*d*dot(p-v1,v2-v1) = min */
2254 /* --> 2d |v2-v1|^2 - 2*dot(p-v1,v2-v1) = 0 */
2255 /* --> d = dot(p-v1,v2-v1) / |v2-v1|^2 */
2256 /***************************************************************/
min_distance_to_line_segment(vector3 p,vector3 v1,vector3 v2)2257 double min_distance_to_line_segment(vector3 p, vector3 v1, vector3 v2) {
2258 vector3 v2mv1 = vector3_minus(v2, v1);
2259 vector3 pmv1 = vector3_minus(p, v1);
2260 double d = vector3_dot(v2mv1, pmv1) / vector3_dot(v2mv1, v2mv1);
2261 if (d < 0.0) d = 0.0; // if pProj lies outside the line segment,
2262 if (d > 1.0) d = 1.0; // displace it to whichever vertex is closer
2263 vector3 pLine = vector3_plus(v1, vector3_scale(d, v2mv1));
2264 return vector3_norm(vector3_minus(p, pLine));
2265 }
2266
2267 /***************************************************************/
2268 /* compute the projection of a 3D point p into the plane */
2269 /* that contains the three points {o, o+v1, o+v2} and has */
2270 /* normal vector v3. */
2271 /* algorithm: solve a 3x3 system to compute the projection of */
2272 /* p into the plane (call it pPlane) */
2273 /* pPlane = p-s*v3 = o + t*v1 + u*v2 */
2274 /* where v3 is the normal to the plane and s,t,u */
2275 /* are unknowns. */
2276 /* the return value is the value of s (where pPlane = p-s*v3), */
2277 /* i.e. the minimum distance from p to the plane. */
2278 /* if in_quadrilateral is non-null it is set to 0 */
2279 /* or 1 according as pPlane does or does not lie in the */
2280 /* quadrilateral with vertices (o, o+v1, o+v2, o+v1+v2). */
2281 /***************************************************************/
normal_distance_to_plane(vector3 p,vector3 o,vector3 v1,vector3 v2,vector3 v3,int * in_quadrilateral)2282 double normal_distance_to_plane(vector3 p, vector3 o, vector3 v1, vector3 v2, vector3 v3,
2283 int *in_quadrilateral) {
2284 CHECK((vector3_norm(v3) > 1.0e-6), "degenerate plane in project_point_into_plane");
2285 matrix3x3 M;
2286 M.c0 = v1;
2287 M.c1 = v2;
2288 M.c2 = v3;
2289 vector3 RHS = vector3_minus(p, o);
2290 vector3 tus = matrix3x3_vector3_mult(matrix3x3_inverse(M), RHS); // "t, u, s"
2291 float t = tus.x, u = tus.y, s = tus.z;
2292 if (in_quadrilateral)
2293 *in_quadrilateral = ((0.0 <= t && t <= 1.0 && 0.0 <= u && u <= 1.0) ? 1 : 0);
2294 return s;
2295 }
2296
2297 // like normal_distance_to_plane, but if pPlane (projection of point into plane)
2298 // lies outside the quadrilateral {o,o+v1,o+v2,o+v1+v2} then take into account
2299 // the in-plane distance from pPlane to the quadrilateral
min_distance_to_quadrilateral(vector3 p,vector3 o,vector3 v1,vector3 v2,vector3 v3)2300 double min_distance_to_quadrilateral(vector3 p, vector3 o, vector3 v1, vector3 v2, vector3 v3) {
2301 int inside;
2302 double s = normal_distance_to_plane(p, o, v1, v2, v3, &inside);
2303 if (inside == 1) return s;
2304 vector3 pPlane = vector3_minus(p, vector3_scale(s, v3));
2305 vector3 p01 = vector3_plus(o, v1);
2306 vector3 p10 = vector3_plus(o, v2);
2307 vector3 p11 = vector3_plus(p01, v2);
2308 double d = min_distance_to_line_segment(pPlane, o, p01);
2309 d = fmin(d, min_distance_to_line_segment(pPlane, o, p10));
2310 d = fmin(d, min_distance_to_line_segment(pPlane, p01, p11));
2311 d = fmin(d, min_distance_to_line_segment(pPlane, p11, p10));
2312 return sqrt(s * s + d * d);
2313 }
2314
2315 // fc==0/1 for floor/ceiling
min_distance_to_prism_roof_or_ceiling(vector3 pp,prism * prsm,int fc)2316 double min_distance_to_prism_roof_or_ceiling(vector3 pp, prism *prsm, int fc) {
2317 int num_vertices = prsm->vertices_p.num_items, i;
2318 vector3 op = {0.0, 0.0, 0.0}; // origin of floor/ceiling
2319 vector3 vps[num_vertices];
2320 if (fc == 1) {
2321 memcpy(vps, prsm->vertices_top_p.items, num_vertices * sizeof(vector3));
2322 for (i = 0; i < num_vertices; i++) {
2323 vps[i].z = 0;
2324 }
2325 op.z = prsm->height;
2326 }
2327 else {
2328 memcpy(vps, prsm->vertices_p.items, num_vertices * sizeof(vector3));
2329 }
2330 vector3 zhatp = {0, 0, 1.0};
2331 double s = normal_distance_to_plane(pp, op, vps[0], vps[1], zhatp, 0);
2332 vector3 ppProj =
2333 vector3_minus(pp, vector3_scale(s, zhatp)); // projection of p into plane of floor/ceiling
2334 if (node_in_polygon(ppProj.x, ppProj.y, vps, num_vertices) == 1) return s;
2335
2336 int nv;
2337 double d = min_distance_to_line_segment(ppProj, vps[0], vps[1]);
2338 for (nv = 1; nv < num_vertices; nv++)
2339 d = fmin(d, min_distance_to_line_segment(ppProj, vps[nv], vps[(nv + 1) % num_vertices]));
2340 return sqrt(s * s + d * d);
2341 }
2342
2343 /***************************************************************/
2344 /* find the face of the prism for which the normal distance */
2345 /* from p to the plane of that face is the shortest, then */
2346 /* return the normal vector to that plane. */
2347 /***************************************************************/
normal_to_prism(prism * prsm,vector3 pc)2348 vector3 normal_to_prism(prism *prsm, vector3 pc) {
2349 if (prsm->height == 0.0) return prsm->axis;
2350
2351 double height = prsm->height;
2352 vector3 *vps_bottom = prsm->vertices_p.items;
2353 vector3 *vps_diff_to_top = prsm->top_polygon_diff_vectors_p.items;
2354 int num_vertices = prsm->vertices_p.num_items;
2355
2356 vector3 zhatp = {0.0, 0.0, 1.0};
2357 vector3 axisp = vector3_scale(height, zhatp);
2358 vector3 pp = prism_coordinate_c2p(prsm, pc);
2359
2360 vector3 retval;
2361 double min_distance = HUGE_VAL;
2362 int nv;
2363 // consider side walls
2364 for (nv = 0; nv < num_vertices; nv++) {
2365 int nvp1 = (nv == (num_vertices - 1) ? 0 : nv + 1);
2366 vector3 v0p = vps_bottom[nv];
2367 vector3 v1p = vector3_minus(vps_bottom[nvp1], vps_bottom[nv]);
2368 vector3 v2p = vps_diff_to_top[nv];
2369 vector3 v3p = unit_vector3(vector3_cross(v1p, v2p));
2370 double s = min_distance_to_quadrilateral(pp, v0p, v1p, v2p, v3p);
2371 if (fabs(s) < min_distance) {
2372 min_distance = fabs(s);
2373 retval = v3p;
2374 }
2375 }
2376
2377 int fc; // 'floor / ceiling'
2378 for (fc = 0; fc < 2; fc++) {
2379 double s = min_distance_to_prism_roof_or_ceiling(pp, prsm, fc);
2380 if (fabs(s) < min_distance) {
2381 min_distance = fabs(s);
2382 retval = zhatp;
2383 }
2384 }
2385 return prism_vector_p2c(prsm, retval);
2386 }
2387
2388 /***************************************************************/
2389 /* Compute the area of a polygon using its vertices. */
2390 /***************************************************************/
get_area_of_polygon_from_nodes(vector3 * nodes,int num_nodes)2391 double get_area_of_polygon_from_nodes(vector3 *nodes, int num_nodes){
2392 double area = 0;
2393 int i;
2394 for (i = 0; i < num_nodes; ++i) {
2395 int i1 = (i + 1) % num_nodes;
2396 area += 0.5 * (nodes[i1].x - nodes[i].x) *
2397 (nodes[i1].y + nodes[i].y);
2398 }
2399 return fabs(area);
2400 }
2401
2402 /***************************************************************/
2403 /* This computes the volume of an irregular triangular prism */
2404 /* following the scheme of http://darrenirvine.blogspot.com/2011/10/volume-of-irregular-triangular-prism.html */
2405 /* The two end triangles have points a, b, and c. Angle abc */
2406 /* is right, and angles bac and acb are acute, of course. The */
2407 /* primary constraint for this method is that the lines be- */
2408 /* tween 'a's between triangles, betweens 'b's between */
2409 /* triangles, and between 'c's between triangles must all be */
2410 /* parallel, though the end triangles need not be parallel. */
2411
2412 /***************************************************************/
get_volume_irregular_triangular_prism(vector3 a0,vector3 b0,vector3 c0,vector3 a1,vector3 b1,vector3 c1)2413 double get_volume_irregular_triangular_prism(vector3 a0, vector3 b0, vector3 c0, vector3 a1, vector3 b1, vector3 c1) {
2414 vector3 side_a = vector3_minus(a1, a0);
2415 vector3 side_b = vector3_minus(b1, b0);
2416 vector3 side_c = vector3_minus(c1, c0);
2417 if (vector3_norm(vector3_cross(side_a, side_b)) != 0 ||
2418 vector3_norm(vector3_cross(side_b, side_c)) != 0 ||
2419 vector3_norm(vector3_cross(side_c, side_a)) != 0) {
2420 //throw an error
2421 }
2422 double length_side_a = vector3_norm(side_a);
2423 double length_side_b = vector3_norm(side_b);
2424 double length_side_c = vector3_norm(side_c);
2425 double average_length = (length_side_a + length_side_b + length_side_c) / 3.0;
2426
2427 vector3 point_on_plane = a0;
2428 vector3 plane_normal_vector = unit_vector3(side_a);
2429 vector3 plane_to_a1 = vector3_minus(a1, point_on_plane);
2430 vector3 plane_to_b1 = vector3_minus(b1, point_on_plane);
2431 vector3 plane_to_c1 = vector3_minus(c1, point_on_plane);
2432 vector3 proj_plane_to_a1_on_normal_vector = vector3_scale(vector3_dot(plane_normal_vector, plane_to_a1), plane_normal_vector);
2433 vector3 proj_plane_to_b1_on_normal_vector = vector3_scale(vector3_dot(plane_normal_vector, plane_to_b1), plane_normal_vector);
2434 vector3 proj_plane_to_c1_on_normal_vector = vector3_scale(vector3_dot(plane_normal_vector, plane_to_c1), plane_normal_vector);
2435 vector3 a1_on_plane = vector3_minus(a1, proj_plane_to_a1_on_normal_vector);
2436 vector3 b1_on_plane = vector3_minus(b1, proj_plane_to_b1_on_normal_vector);
2437 vector3 c1_on_plane = vector3_minus(c1, proj_plane_to_c1_on_normal_vector);
2438 double base = vector3_norm(vector3_minus(c1_on_plane, b1_on_plane));
2439 double height = vector3_norm(vector3_minus(a1_on_plane, b1_on_plane));
2440 double cross_sectional_area = fabs(0.5 * base * height);
2441
2442 return fabs(average_length * cross_sectional_area);
2443 }
2444
2445 /***************************************************************/
2446 /***************************************************************/
2447 /***************************************************************/
get_prism_volume(prism * prsm)2448 double get_prism_volume(prism *prsm) {
2449 if (prsm->sidewall_angle == 0.0) {
2450 return get_area_of_polygon_from_nodes(prsm->vertices_p.items, prsm->vertices_p.num_items) * fabs(prsm->height);
2451 }
2452 else {
2453 int num_vertices = prsm->vertices_p.num_items, nv;
2454 double bottom_polygon_area = get_area_of_polygon_from_nodes(prsm->vertices_p.items, prsm->vertices_p.num_items);
2455 double top_polygon_area = get_area_of_polygon_from_nodes(prsm->vertices_top_p.items, prsm->vertices_top_p.num_items);
2456 double volume;
2457 vector3 *wedges_a;
2458 wedges_a = (vector3 *)malloc(num_vertices * sizeof(vector3));
2459 CHECK(wedges_a, "out of memory");
2460 vector3 *wedges_b;
2461 wedges_b = (vector3 *)malloc(num_vertices * sizeof(vector3));
2462 CHECK(wedges_b, "out of memory");
2463 vector3 *wedges_c;
2464 wedges_c = (vector3 *)malloc(num_vertices * sizeof(vector3));
2465 CHECK(wedges_c, "out of memory");
2466 if (bottom_polygon_area > top_polygon_area) {
2467 volume = fabs(top_polygon_area * prsm->height);
2468 memcpy(wedges_a, prsm->vertices_top_p.items, num_vertices * sizeof(vector3));
2469 memcpy(wedges_b, prsm->vertices_top_p.items, num_vertices * sizeof(vector3));
2470 for (nv = 0; nv < num_vertices; nv++) {
2471 wedges_b[nv].z = 0.0;
2472 }
2473 memcpy(wedges_c, prsm->vertices_p.items, num_vertices * sizeof(vector3));
2474 }
2475 else {
2476 volume = fabs(bottom_polygon_area * prsm->height);
2477 memcpy(wedges_a, prsm->vertices_p.items, num_vertices * sizeof(vector3));
2478 memcpy(wedges_b, prsm->vertices_p.items, num_vertices * sizeof(vector3));
2479 for (nv = 0; nv < num_vertices; nv++) {
2480 wedges_b[nv].z = prsm->height;
2481 }
2482 memcpy(wedges_c, prsm->vertices_top_p.items, num_vertices * sizeof(vector3));
2483 }
2484 for (nv = 0; nv < num_vertices; nv++) {
2485 int nvp1 = (nv + 1 == num_vertices ? 0 : nv + 1);
2486 volume += get_volume_irregular_triangular_prism(wedges_a[nv], wedges_b[nv], wedges_c[nv], wedges_a[nvp1], wedges_b[nvp1], wedges_c[nvp1]);
2487 }
2488 return volume;
2489 }
2490 }
2491
2492 /***************************************************************/
2493 /***************************************************************/
2494 /***************************************************************/
get_prism_bounding_box(prism * prsm,geom_box * box)2495 void get_prism_bounding_box(prism *prsm, geom_box *box) {
2496 vector3 *vertices = prsm->vertices.items;
2497 vector3 *vertices_top = prsm->vertices_top.items;
2498 int num_vertices = prsm->vertices.num_items;
2499 box->low = box->high = vertices[0];
2500 int nv, fc;
2501 for (nv = 0; nv < num_vertices; nv++)
2502 for (fc = 0; fc < 2; fc++) // 'floor,ceiling'
2503 {
2504 vector3 v;
2505 if (fc == 0) v = vertices[nv];
2506 if (fc == 1) v = vertices_top[nv];
2507
2508 box->low.x = fmin(box->low.x, v.x);
2509 box->low.y = fmin(box->low.y, v.y);
2510 box->low.z = fmin(box->low.z, v.z);
2511
2512 box->high.x = fmax(box->high.x, v.x);
2513 box->high.y = fmax(box->high.y, v.y);
2514 box->high.z = fmax(box->high.z, v.z);
2515 }
2516 }
2517
2518 /***************************************************************/
2519 /***************************************************************/
2520 /***************************************************************/
display_prism_info(int indentby,geometric_object * o)2521 void display_prism_info(int indentby, geometric_object *o) {
2522 prism *prsm = o->subclass.prism_data;
2523
2524 vector3 *vs = prsm->vertices.items;
2525 int num_vertices = prsm->vertices.num_items;
2526
2527 ctl_printf("%*s height %g, axis (%g,%g,%g), sidewall angle: %g radians, %i vertices:\n", indentby, "", prsm->height,
2528 prsm->axis.x, prsm->axis.y, prsm->axis.z, prsm->sidewall_angle, num_vertices);
2529 int nv;
2530 for (nv = 0; nv < num_vertices; nv++)
2531 ctl_printf("%*s (%g,%g,%g)\n", indentby, "", vs[nv].x, vs[nv].y, vs[nv].z);
2532 }
2533
2534 /***************************************************************/
2535 // like vector3_equal but tolerant of small floating-point discrepancies
2536 /***************************************************************/
vector3_nearly_equal(vector3 v1,vector3 v2,double tolerance)2537 int vector3_nearly_equal(vector3 v1, vector3 v2, double tolerance) {
2538 return (vector3_norm(vector3_minus(v1, v2)) <= tolerance * vector3_norm(v1));
2539 }
2540 matrix3x3 sidewall_scaling_matrix;
2541
2542 /***************************************************************/
2543 /* return the unit normal vector to the triangle with the given*/
2544 /* vertices */
2545 /***************************************************************/
triangle_normal(vector3 v1,vector3 v2,vector3 v3)2546 vector3 triangle_normal(vector3 v1, vector3 v2, vector3 v3) {
2547 vector3 nv = vector3_cross(vector3_minus(v2, v1), vector3_minus(v3, v1));
2548 double nvnorm = vector3_norm(nv);
2549 // if (area) *area += 0.5*nvnorm;
2550 return unit_vector3(nv);
2551 }
2552
2553 /***************************************************************/
2554 /* On entry, the only fields in o->prism that are assumed to */
2555 /* be initialized are: vertices, height, (optionally) */
2556 /* axis, and sidewall_angle. If axis has not been initialized */
2557 /* (i.e. it is set to its default value, which is the zero */
2558 /* vector) then the prism axis is automatically computed as */
2559 /* the normal to the vertex plane. If o->center is equal to */
2560 /* auto_center on entry, then it is set to the prism center, */
2561 /* as computed from the vertices, axis, and height. Otherwise, */
2562 /* the prism is rigidly translated to center it at the */
2563 /* specified value of o->center. */
2564 /***************************************************************/
2565 // special vector3 that signifies 'no value specified'
2566 vector3 auto_center = {NAN, NAN, NAN};
init_prism(geometric_object * o)2567 void init_prism(geometric_object *o) {
2568 prism *prsm = o->subclass.prism_data;
2569 vector3 *vertices = prsm->vertices.items;
2570 int num_vertices = prsm->vertices.num_items;
2571 CHECK(num_vertices >= 3, "fewer than 3 vertices in init_prism");
2572
2573 // compute centroid of vertices
2574 vector3 centroid = {0.0, 0.0, 0.0};
2575 int nv;
2576 for (nv = 0; nv < num_vertices; nv++)
2577 centroid = vector3_plus(centroid, vertices[nv]);
2578 prsm->centroid = centroid = vector3_scale(1.0 / ((double)num_vertices), centroid);
2579
2580 // make sure all vertices lie in a plane, i.e. that the normal
2581 // vectors to all triangles (v_n, v_{n+1}, centroid) agree.
2582 int plane_normal_set = 0;
2583 vector3 plane_normal;
2584 double tol = 1.0e-6;
2585 for (nv = 0; nv < num_vertices; nv++) {
2586 int nvp1 = (nv + 1) % num_vertices;
2587 vector3 tri_normal = triangle_normal(centroid, vertices[nv], vertices[nvp1]);
2588 if (vector3_norm(tri_normal) == 0.0) // vertices collinear with centroid
2589 continue;
2590 if (!plane_normal_set) {
2591 plane_normal = tri_normal;
2592 plane_normal_set = 1;
2593 }
2594 else {
2595 boolean normals_agree =
2596 (vector3_nearly_equal(plane_normal, tri_normal, tol) ||
2597 vector3_nearly_equal(plane_normal, vector3_scale(-1.0, tri_normal), tol));
2598 CHECK(normals_agree, "non-coplanar vertices in init_prism");
2599 }
2600 }
2601
2602 // if no prism axis was specified, set the prism axis equal to the
2603 // normal to the vertex plane.
2604 // if a prism axis was specified, check that it agrees up to sign
2605 // with the normal to the vertex plane.
2606 if (vector3_norm(prsm->axis) == 0.0)
2607 prsm->axis = plane_normal;
2608 else {
2609 prsm->axis = unit_vector3(prsm->axis);
2610 boolean axis_normal_to_plane =
2611 (vector3_nearly_equal(prsm->axis, plane_normal, tol) ||
2612 vector3_nearly_equal(prsm->axis, vector3_scale(-1.0, plane_normal), tol));
2613 CHECK(axis_normal_to_plane, "axis not normal to vertex plane in init_prism");
2614 }
2615
2616 // set current_center=prism center as determined by vertices and height.
2617 // if the center of the geometric object was left unspecified,
2618 // set it to current_center; otherwise displace the entire prism
2619 // so that it is centered at the specified center.
2620 vector3 current_center = vector3_plus(centroid, vector3_scale(0.5 * prsm->height, prsm->axis));
2621 if (isnan(o->center.x) && isnan(o->center.y) && isnan(o->center.z)) // center == auto-center
2622 o->center = current_center;
2623 else {
2624 vector3 shift = vector3_minus(o->center, current_center);
2625 for (nv = 0; nv < num_vertices; nv++)
2626 vertices[nv] = vector3_plus(vertices[nv], shift);
2627 centroid = vector3_plus(centroid, shift);
2628 }
2629
2630 // compute rotation matrix that operates on a vector of cartesian coordinates
2631 // to yield the coordinates of the same point in the prism coordinate system.
2632 // the prism coordinate system is a right-handed coordinate system
2633 // in which the prism lies in the xy plane (extrusion axis is the positive z-axis)
2634 // with centroid at the origin.
2635 // note: the prism *centroid* is the center of mass of the planar vertex polygon,
2636 // i.e. it is a point lying on the bottom surface of the prism.
2637 // This is the origin of coordinates in the prism system.
2638 // The *center* of the geometric object is the center of mass of the
2639 // 3D prism. So center = centroid + 0.5*height*zHat.
2640 vector3 x0hat = {1.0, 0.0, 0.0}, y0hat = {0.0, 1.0, 0.0}, z0hat = {0.0, 0.0, 1.0};
2641 vector3 xhat, yhat, zhat = prsm->axis;
2642 if (vector3_nearly_equal(zhat, x0hat, tol)) {
2643 xhat = y0hat;
2644 yhat = z0hat;
2645 }
2646 else if (vector3_nearly_equal(zhat, y0hat, tol)) {
2647 xhat = z0hat;
2648 yhat = x0hat;
2649 }
2650 else if (vector3_nearly_equal(zhat, z0hat, tol)) {
2651 xhat = x0hat;
2652 yhat = y0hat;
2653 }
2654 else {
2655 xhat = unit_vector3(vector3_minus(vertices[1], vertices[0]));
2656 yhat = unit_vector3(vector3_cross(zhat, xhat));
2657 }
2658 matrix3x3 m_p2c = {xhat, yhat, zhat};
2659 prsm->m_p2c = m_p2c;
2660 prsm->m_c2p = matrix3x3_inverse(m_p2c);
2661
2662 // compute vertices in prism coordinate system
2663 prsm->vertices_p.num_items = num_vertices;
2664 prsm->vertices_p.items = (vector3 *)malloc(num_vertices * sizeof(vector3));
2665 for (nv = 0; nv < num_vertices; nv++)
2666 prsm->vertices_p.items[nv] = prism_coordinate_c2p(prsm, vertices[nv]);
2667
2668 // Calculate difference vertices of the top polygon and vectors between bottom
2669 // polygon and the top polygon, where:
2670 // * the bottom polygon is the one passed in to the the make_prism() function,
2671 // stored in vertices and vertices_p,
2672 // * the top polygon is the top surface (parallel to the bottom polygon) resulting
2673 // from the extrusion of the bottom polygon. Whether or not the extrusion tapers
2674 // is dependent on the value of sidewall_angle.
2675 //
2676 // The top polygon is calculated by first copying the values of vertices_p into
2677 // vertices_top_p, except z=prsm->height for all top vertices. If prsm->sidewall_angle
2678 // is equal to zero, then no further calculations are performed on the top vertices.
2679 // If not, we know that all EDGES of the the top polygon will be offset so that in the
2680 // xy plane they are parallel to the edges of the bottom polygon. The offset amount is
2681 // determined by the sidewall angle and the height of the prism. To perform the
2682 // calculation, each of the edges of the top polygon (without an offset) are stored in
2683 // an array of edges (edge is a struct defined if prsm->sidewall_angle!=0 containing
2684 // the endpoints a1 a2, with a third vector v defined a2-a1). Then the vector normal to
2685 // v is calculated, and the offset vector. A test is performed to determine in which
2686 // direction (the direction of +offset or -offset) from the edge we can find points
2687 // inside the polygon by performing a node_in_or_on_polygon test at a finite distance
2688 // away from the midpoint of the edge:
2689 // edge.a1 + 0.5*edge.v + 1e-3*offset.
2690 // This information is used to determine in which direction the offset of the edge is
2691 // applied, in conjunction with whether prsm->sidewall_angle is positive or negative
2692 // (if positive, the offset will be applied in towards the points where
2693 // node_in_or_on_polygon is true, else the offset will be applied out away from those
2694 // points). After the offsets are applied to the edges, the intersections between the
2695 // new edges are calculated, which are the new values of vertices_top_p.
2696 //
2697 // Some side notes on the difference vectors:
2698 // * The value of each of the top polygon vertices can be found
2699 // vertices_p + top_polygon_diff_vectors_p
2700 // vertices + top_polygon_diff_vectors
2701 // * A linearly interpolated value of the polygon vertices between the bottom
2702 // polygon and the top can be found
2703 // vertices_p + top_polygon_diff_vectors_scaled_p * z
2704 number theta = (K_PI/2) - fabs(prsm->sidewall_angle);
2705 prsm->vertices_top_p.num_items = num_vertices;
2706 prsm->vertices_top_p.items = (vector3 *)malloc(num_vertices * sizeof(vector3));
2707 CHECK(prsm->vertices_top_p.items, "out of memory");
2708 memcpy(prsm->vertices_top_p.items, prsm->vertices_p.items, num_vertices * sizeof(vector3));
2709 for (nv = 0; nv < num_vertices; nv++) {
2710 prsm->vertices_top_p.items[nv].z = prsm->height;
2711 }
2712
2713 if (prsm->sidewall_angle != 0.0) {
2714 typedef struct {
2715 vector3 a1, a2, v; // v will be defined as a2 - a1
2716 } edge;
2717
2718 // find the point at the bottom left corner of the polygon
2719 double smallest_x = HUGE_VAL;
2720 double smallest_y = HUGE_VAL;
2721 int index_for_point_a = -1;
2722 int index_for_point_b = -1;
2723 int index_for_point_c = -1;
2724 for (nv = 0; nv < num_vertices; nv++) {
2725 double current_x = prsm->vertices_p.items[nv].x;
2726 double current_y = prsm->vertices_p.items[nv].y;
2727 if (current_x < smallest_x) {
2728 smallest_x = current_x;
2729 smallest_y = current_y;
2730 index_for_point_b = nv;
2731 }
2732 else if (current_x == smallest_x && current_y < smallest_y) {
2733 smallest_y = current_y;
2734 index_for_point_b = nv;
2735 }
2736 }
2737 if (index_for_point_b == -1) {
2738 exit(EXIT_FAILURE);
2739 }
2740 else {
2741 index_for_point_a = (index_for_point_b + 1 == num_vertices ? 0 : index_for_point_b + 1);
2742 index_for_point_c = (index_for_point_b - 1 == -1 ? num_vertices - 1 : index_for_point_b - 1);
2743 }
2744 // find orientation of the polygon
2745 vector3 A = prsm->vertices_p.items[index_for_point_a];
2746 vector3 B = prsm->vertices_p.items[index_for_point_b];
2747 vector3 C = prsm->vertices_p.items[index_for_point_c];
2748 double orientation_number = (B.x - A.x)*(C.y - A.y)-(C.x - A.x)*(B.y - A.y);
2749 int orientation_positive_or_negative = (orientation_number < 0 ? 0 : 1);
2750
2751 edge *top_polygon_edges;
2752 top_polygon_edges = (edge *)malloc(num_vertices * sizeof(edge));
2753 number w = prsm->height / tan(theta);
2754
2755 for (nv = 0; nv < num_vertices; nv++) {
2756 top_polygon_edges[nv].a1 = prsm->vertices_top_p.items[(nv - 1 == -1 ? num_vertices - 1 : nv - 1)];
2757 top_polygon_edges[nv].a2 = prsm->vertices_top_p.items[nv];
2758 top_polygon_edges[nv].v = vector3_minus(top_polygon_edges[nv].a2, top_polygon_edges[nv].a1);
2759
2760 vector3 normal_vector = (orientation_positive_or_negative ? unit_vector3(vector3_cross(top_polygon_edges[nv].v, zhat)) : unit_vector3(vector3_cross(top_polygon_edges[nv].v, vector3_scale(-1, zhat))));
2761
2762 // positive sidewall angles means the prism tapers in towards the rest of the prism body
2763 // negative sidewall angles means the prism tapers out away from the rest of the prism body
2764 vector3 offset = vector3_scale(prsm->sidewall_angle > 0 ? w : -w, normal_vector);
2765 top_polygon_edges[nv].a1 = vector3_plus(top_polygon_edges[nv].a1, offset);
2766 top_polygon_edges[nv].a2 = vector3_plus(top_polygon_edges[nv].a2, offset);
2767 }
2768
2769 for (nv = 0; nv < num_vertices; nv++) {
2770 number x1 = top_polygon_edges[nv].a1.x;
2771 number y1 = top_polygon_edges[nv].a1.y;
2772 number x2 = top_polygon_edges[nv].a2.x;
2773 number y2 = top_polygon_edges[nv].a2.y;
2774 number x3 = top_polygon_edges[(nv + 1 == num_vertices ? 0 : nv + 1)].a1.x;
2775 number y3 = top_polygon_edges[(nv + 1 == num_vertices ? 0 : nv + 1)].a1.y;
2776 number x4 = top_polygon_edges[(nv + 1 == num_vertices ? 0 : nv + 1)].a2.x;
2777 number y4 = top_polygon_edges[(nv + 1 == num_vertices ? 0 : nv + 1)].a2.y;
2778
2779 // Intersection point calculated with https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection#Given_two_points_on_each_line
2780 number px = ((x1*y2-y1*x2)*(x3-x4)-(x1-x2)*(x3*y4-y3*x4)) / ((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4));
2781 number py = ((x1*y2-y1*x2)*(y3-y4)-(y1-y2)*(x3*y4-y3*x4)) / ((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4));
2782 prsm->vertices_top_p.items[nv].x = px;
2783 prsm->vertices_top_p.items[nv].y = py;
2784 }
2785 }
2786
2787 prsm->top_polygon_diff_vectors_p.num_items = num_vertices;
2788 prsm->top_polygon_diff_vectors_p.items = (vector3 *)malloc(num_vertices * sizeof(vector3));
2789 CHECK(prsm->top_polygon_diff_vectors_p.items, "out of memory");
2790 for (nv = 0; nv < num_vertices; nv++) {
2791 prsm->top_polygon_diff_vectors_p.items[nv] = vector3_minus(prsm->vertices_top_p.items[nv], prsm->vertices_p.items[nv]);
2792 }
2793
2794 prsm->top_polygon_diff_vectors_scaled_p.num_items = num_vertices;
2795 prsm->top_polygon_diff_vectors_scaled_p.items = (vector3 *)malloc(num_vertices * sizeof(vector3));
2796 CHECK(prsm->top_polygon_diff_vectors_scaled_p.items, "out of memory");
2797 for (nv = 0; nv < num_vertices; nv++) {
2798 prsm->top_polygon_diff_vectors_scaled_p.items[nv] = vector3_scale(1/prsm->height, prsm->top_polygon_diff_vectors_p.items[nv]);
2799 }
2800
2801 prsm->vertices_top.num_items = num_vertices;
2802 prsm->vertices_top.items = (vector3 *)malloc(num_vertices * sizeof(vector3));
2803 CHECK(prsm->vertices_top.items, "out of memory");
2804 for (nv = 0; nv < num_vertices; nv++) {
2805 prsm->vertices_top.items[nv] = prism_coordinate_p2c(prsm, prsm->vertices_top_p.items[nv]);
2806 }
2807
2808 // workspace is an internally-stored double-valued array of length num_vertices+2
2809 // that is used by some geometry routines
2810 prsm->workspace.num_items = num_vertices + 2;
2811 prsm->workspace.items = (double *)malloc((num_vertices + 2) * sizeof(double));
2812 }
2813
2814 /***************************************************************/
2815 /* routines called from C++ or python codes to create prisms */
2816 /***************************************************************/
2817 // prism with center determined automatically from vertices, height, and axis
make_prism(material_type material,const vector3 * vertices,int num_vertices,double height,vector3 axis)2818 geometric_object make_prism(material_type material, const vector3 *vertices, int num_vertices,
2819 double height, vector3 axis) {
2820 return make_prism_with_center(material, auto_center, vertices, num_vertices, height, axis);
2821 }
2822
2823 // prism in which all vertices are translated to ensure that the prism is centered at center.
make_prism_with_center(material_type material,vector3 center,const vector3 * vertices,int num_vertices,double height,vector3 axis)2824 geometric_object make_prism_with_center(material_type material, vector3 center, const vector3 *vertices,
2825 int num_vertices, double height, vector3 axis) {
2826 return make_slanted_prism_with_center(material, center, vertices, num_vertices, height, axis, 0);
2827 }
2828
2829 // slanted prism with center determined automatically from vertices, height, axis, and sidewall_angle
make_slanted_prism(material_type material,const vector3 * vertices,int num_vertices,double height,vector3 axis,double sidewall_angle)2830 geometric_object make_slanted_prism(material_type material, const vector3 *vertices,
2831 int num_vertices, double height, vector3 axis, double sidewall_angle) {
2832 return make_slanted_prism_with_center(material, auto_center, vertices, num_vertices, height, axis, sidewall_angle);
2833 }
2834
2835 // Have both make_prism_with_center and make_slanted_prism_with_center keep the same parameters to maintain ABI
2836 // compatibility, though make_prism_with_center just calls make_slanted_prism_with_center with the sidewall angle equal
2837 // to zero. To make a slanted prism, the user will have to call make_slanted_prism for now.
make_slanted_prism_with_center(material_type material,vector3 center,const vector3 * vertices,int num_vertices,double height,vector3 axis,double sidewall_angle)2838 geometric_object make_slanted_prism_with_center(material_type material, vector3 center, const vector3 *vertices,
2839 int num_vertices, double height, vector3 axis, double sidewall_angle) {
2840 geometric_object o = make_geometric_object(material, center);
2841 o.which_subclass = GEOM PRISM;
2842 prism *prsm = o.subclass.prism_data = MALLOC1(prism);
2843 CHECK(prsm, "out of memory");
2844 prsm->vertices.num_items = num_vertices;
2845 prsm->vertices.items = (vector3 *)malloc(num_vertices * sizeof(vector3));
2846 CHECK(prsm->vertices.items, "out of memory");
2847 memcpy(prsm->vertices.items, vertices, num_vertices * sizeof(vector3));
2848 prsm->height = height;
2849 prsm->axis = axis;
2850 prsm->sidewall_angle = sidewall_angle;
2851 init_prism(&o);
2852 return o;
2853 }
2854