Hello community, here is the log from the commit of package libctl for openSUSE:Factory checked in at 2019-03-18 10:41:23 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Comparing /work/SRC/openSUSE:Factory/libctl (Old) and /work/SRC/openSUSE:Factory/.libctl.new.28833 (New) ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Package is "libctl" Mon Mar 18 10:41:23 2019 rev:3 rq:685482 version:4.2.0 Changes: -------- --- /work/SRC/openSUSE:Factory/libctl/libctl.changes 2018-11-26 10:28:35.525148618 +0100 +++ /work/SRC/openSUSE:Factory/.libctl.new.28833/libctl.changes 2019-03-18 10:41:23.835270180 +0100 @@ -1,0 +2,6 @@ +Tue Feb 19 09:33:56 UTC 2019 - [email protected] + +- Update to version 4.2.0 + * Better handling of center parameter of prisms + +------------------------------------------------------------------- Old: ---- libctl-4.1.4.tar.gz New: ---- libctl-4.2.0.tar.gz ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Other differences: ------------------ ++++++ libctl.spec ++++++ --- /var/tmp/diff_new_pack.Ca69nA/_old 2019-03-18 10:41:24.347269555 +0100 +++ /var/tmp/diff_new_pack.Ca69nA/_new 2019-03-18 10:41:24.347269555 +0100 @@ -1,7 +1,7 @@ # # spec file for package libctl # -# Copyright (c) 2018 SUSE LINUX GmbH, Nuernberg, Germany. +# Copyright (c) 2019 SUSE LINUX GmbH, Nuernberg, Germany. # # All modifications and additions to the file contributed by third parties # remain the property of their copyright owners, unless otherwise agreed @@ -17,7 +17,7 @@ Name: libctl -Version: 4.1.4 +Version: 4.2.0 Release: 0 %define somajor 5 Summary: A guile Library for Scientific Simulations ++++++ libctl-4.1.4.tar.gz -> libctl-4.2.0.tar.gz ++++++ diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/libctl-4.1.4/NEWS.md new/libctl-4.2.0/NEWS.md --- old/libctl-4.1.4/NEWS.md 2018-11-16 22:03:44.000000000 +0100 +++ new/libctl-4.2.0/NEWS.md 2019-01-07 20:44:56.000000000 +0100 @@ -1,5 +1,13 @@ # Libctl Release Notes +## libctl 4.2.0 + +* Better handling of `center` parameter of prisms, allowing this + to be optionally specified (#35). Deprecates old `geom_fix_object` + and `geom_fix_objects0` functions in favor of `geom_fix_object_ptr` + and `geom_fix_object_list`. (In particular, the old `geom_fix_object` routine will not work properly for prisms in + which the center was not computed.) + ## libctl 4.1.4 11/16/18 diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/libctl-4.1.4/configure new/libctl-4.2.0/configure --- old/libctl-4.1.4/configure 2018-11-16 22:04:33.000000000 +0100 +++ new/libctl-4.2.0/configure 2019-01-07 20:48:52.000000000 +0100 @@ -1,6 +1,6 @@ #! /bin/sh # Guess values for system-dependent variables and create Makefiles. -# Generated by GNU Autoconf 2.69 for libctl 4.1.4. +# Generated by GNU Autoconf 2.69 for libctl 4.2.0. # # Report bugs to <[email protected]>. # @@ -590,8 +590,8 @@ # Identity of this package. PACKAGE_NAME='libctl' PACKAGE_TARNAME='libctl' -PACKAGE_VERSION='4.1.4' -PACKAGE_STRING='libctl 4.1.4' +PACKAGE_VERSION='4.2.0' +PACKAGE_STRING='libctl 4.2.0' PACKAGE_BUGREPORT='[email protected]' PACKAGE_URL='' @@ -1335,7 +1335,7 @@ # Omit some internal or obsolete options to make the list less imposing. # This message is too long to be a string in the A/UX 3.1 sh. cat <<_ACEOF -\`configure' configures libctl 4.1.4 to adapt to many kinds of systems. +\`configure' configures libctl 4.2.0 to adapt to many kinds of systems. Usage: $0 [OPTION]... [VAR=VALUE]... @@ -1405,7 +1405,7 @@ if test -n "$ac_init_help"; then case $ac_init_help in - short | recursive ) echo "Configuration of libctl 4.1.4:";; + short | recursive ) echo "Configuration of libctl 4.2.0:";; esac cat <<\_ACEOF @@ -1522,7 +1522,7 @@ test -n "$ac_init_help" && exit $ac_status if $ac_init_version; then cat <<\_ACEOF -libctl configure 4.1.4 +libctl configure 4.2.0 generated by GNU Autoconf 2.69 Copyright (C) 2012 Free Software Foundation, Inc. @@ -1975,7 +1975,7 @@ This file contains any messages produced by compilers while running configure, to aid debugging if configure makes a mistake. -It was created by libctl $as_me 4.1.4, which was +It was created by libctl $as_me 4.2.0, which was generated by GNU Autoconf 2.69. Invocation command line was $ $0 $@ @@ -2354,7 +2354,7 @@ # Shared-library version number; indicates api compatibility, and is # not the same as the "public" version number. (Don't worry about this # except for public releases.) -SHARED_VERSION_INFO="7:2:0" # CURRENT:REVISION:AGE +SHARED_VERSION_INFO="7:3:0" # CURRENT:REVISION:AGE am__api_version='1.16' @@ -2871,7 +2871,7 @@ # Define the identity of the package. PACKAGE='libctl' - VERSION='4.1.4' + VERSION='4.2.0' cat >>confdefs.h <<_ACEOF @@ -18374,7 +18374,7 @@ # report actual input values of CONFIG_FILES etc. instead of their # values after options handling. ac_log=" -This file was extended by libctl $as_me 4.1.4, which was +This file was extended by libctl $as_me 4.2.0, which was generated by GNU Autoconf 2.69. Invocation command line was CONFIG_FILES = $CONFIG_FILES @@ -18440,7 +18440,7 @@ cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1 ac_cs_config="`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`" ac_cs_version="\\ -libctl config.status 4.1.4 +libctl config.status 4.2.0 configured by $0, generated by GNU Autoconf 2.69, with options \\"\$ac_cs_config\\" diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/libctl-4.1.4/configure.ac new/libctl-4.2.0/configure.ac --- old/libctl-4.1.4/configure.ac 2018-11-16 22:03:32.000000000 +0100 +++ new/libctl-4.2.0/configure.ac 2019-01-07 20:45:12.000000000 +0100 @@ -1,5 +1,5 @@ # Process this file with autoconf to produce a configure script. -AC_INIT(libctl, 4.1.4, [email protected]) +AC_INIT(libctl, 4.2.0, [email protected]) AC_CONFIG_SRCDIR([src/ctl.c]) AC_CONFIG_HEADERS([config.h src/ctl.h]) AC_CONFIG_MACRO_DIR([m4]) @@ -8,7 +8,7 @@ # Shared-library version number; indicates api compatibility, and is # not the same as the "public" version number. (Don't worry about this # except for public releases.) -SHARED_VERSION_INFO="7:2:0" # CURRENT:REVISION:AGE +SHARED_VERSION_INFO="7:3:0" # CURRENT:REVISION:AGE AM_INIT_AUTOMAKE([foreign]) AC_SUBST(SHARED_VERSION_INFO) diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/libctl-4.1.4/utils/ctlgeom-types.h new/libctl-4.2.0/utils/ctlgeom-types.h --- old/libctl-4.1.4/utils/ctlgeom-types.h 2018-11-16 22:04:50.000000000 +0100 +++ new/libctl-4.2.0/utils/ctlgeom-types.h 2019-01-07 20:49:24.000000000 +0100 @@ -90,8 +90,10 @@ typedef struct prism_struct { vector3_list vertices; - vector3 centroid; number height; + vector3 axis; + vector3_list vertices_p; + vector3 centroid; number_list workspace; matrix3x3 m_c2p; matrix3x3 m_p2c; diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/libctl-4.1.4/utils/ctlgeom.h new/libctl-4.2.0/utils/ctlgeom.h --- old/libctl-4.1.4/utils/ctlgeom.h 2018-06-14 17:46:09.000000000 +0200 +++ new/libctl-4.2.0/utils/ctlgeom.h 2019-01-07 20:41:38.000000000 +0100 @@ -39,6 +39,14 @@ # define MATERIAL_TYPE void* #endif +/* Where possible (e.g. for gcc >= 3.1), enable a compiler warning + for code that uses a deprecated function */ +#if defined(__GNUC__) && (__GNUC__ > 3 || (__GNUC__==3 && __GNUC_MINOR__ > 0)) +# define CTLGEOM_DEPRECATED __attribute__((deprecated)) +#else +# define CTLGEOM_DEPRECATED +#endif + #ifdef __cplusplus extern "C" { #endif /* __cplusplus */ @@ -51,9 +59,11 @@ #endif extern void geom_initialize(void); -extern void geom_fix_object(GEOMETRIC_OBJECT o); +extern void geom_fix_object_ptr(GEOMETRIC_OBJECT *o); +extern void geom_fix_object(GEOMETRIC_OBJECT o) CTLGEOM_DEPRECATED; extern void geom_fix_objects(void); -extern void geom_fix_objects0(GEOMETRIC_OBJECT_LIST geometry); +extern void geom_fix_objects0(GEOMETRIC_OBJECT_LIST geometry) CTLGEOM_DEPRECATED; +extern void geom_fix_object_list(GEOMETRIC_OBJECT_LIST geometry); extern void geom_fix_lattice(void); extern void geom_fix_lattice0(LATTICE *L); extern void geom_cartesian_lattice(void); @@ -147,10 +157,18 @@ GEOMETRIC_OBJECT make_ellipsoid(MATERIAL_TYPE material, vector3 center, vector3 e1, vector3 e2, vector3 e3, vector3 size); + +// prism with `center` field computed automatically from vertices, height, axis GEOMETRIC_OBJECT make_prism(MATERIAL_TYPE material, const vector3 *vertices, int num_vertices, double height, vector3 axis); +// as make_prism, but with a rigid translation so that the prism is centered at center +GEOMETRIC_OBJECT make_prism_with_center(MATERIAL_TYPE material, vector3 center, + const vector3 *vertices, int num_vertices, + double height, vector3 axis); + + int vector3_nearly_equal(vector3 v1, vector3 v2, double tolerance); /**************************************************************************/ diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/libctl-4.1.4/utils/geom-ctl-io.c new/libctl-4.2.0/utils/geom-ctl-io.c --- old/libctl-4.1.4/utils/geom-ctl-io.c 2018-11-16 22:04:50.000000000 +0100 +++ new/libctl-4.2.0/utils/geom-ctl-io.c 2019-01-07 20:49:24.000000000 +0100 @@ -58,8 +58,17 @@ o->vertices.items[i_t] = o0->vertices.items[i_t]; } } - o->centroid = o0->centroid; o->height = o0->height; + o->axis = o0->axis; + { + int i_t; + o->vertices_p.num_items = o0->vertices_p.num_items; + o->vertices_p.items = ((vector3 *) malloc(sizeof(vector3) * (o->vertices_p.num_items))); + for (i_t = 0; i_t < o->vertices_p.num_items; i_t++) { + o->vertices_p.items[i_t] = o0->vertices_p.items[i_t]; + } + } + o->centroid = o0->centroid; { int i_t; o->workspace.num_items = o0->workspace.num_items; @@ -219,10 +228,21 @@ return 0; } } - if (!vector3_equal(o->centroid, o0->centroid)) - return 0; if (o->height != o0->height) return 0; + if (!vector3_equal(o->axis, o0->axis)) + return 0; + { + int i_t; + if (o->vertices_p.num_items != o0->vertices_p.num_items) + return 0; + for (i_t = 0; i_t < o->vertices_p.num_items; i_t++) { + if (!vector3_equal(o->vertices_p.items[i_t], o0->vertices_p.items[i_t])) + return 0; + } + } + if (!vector3_equal(o->centroid, o0->centroid)) + return 0; { int i_t; if (o->workspace.num_items != o0->workspace.num_items) @@ -383,6 +403,12 @@ free(o.vertices.items); { int index_t; + for (index_t = 0; index_t < o.vertices_p.num_items; index_t++) { + } + } + free(o.vertices_p.items); + { + int index_t; for (index_t = 0; index_t < o.workspace.num_items; index_t++) { } } diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/libctl-4.1.4/utils/geom.c new/libctl-4.2.0/utils/geom.c --- old/libctl-4.1.4/utils/geom.c 2018-11-08 21:37:49.000000000 +0100 +++ new/libctl-4.2.0/utils/geom.c 2019-01-07 20:39:42.000000000 +0100 @@ -66,12 +66,13 @@ #define MIN(a,b) ((a) < (b) ? (a) : (b)) // forward declarations of prism-related routines, at the bottom of this file +static boolean node_in_polygon(double qx, double qy, vector3 *nodes, int num_nodes); +static boolean point_in_prism(prism *prsm, vector3 pc); +static vector3 normal_to_prism(prism *prsm, vector3 pc); +static double intersect_line_segment_with_prism(prism *prsm, vector3 pc, vector3 dc, double a, double b); static void get_prism_bounding_box(prism *prsm, geom_box *box); -static double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b); -static boolean point_in_polygon(double px, double py, vector3 *vertices, int num_vertices); -static boolean point_in_prism(prism *prsm, vector3 p); -static vector3 normal_to_prism(prism *prsm, vector3 p); -static void display_prism_info(int indentby, prism *prsm); +static void display_prism_info(int indentby, geometric_object *o); +static void init_prism(geometric_object *o); /**************************************************************************/ /* If v is a vector in the lattice basis, normalize v so that @@ -97,7 +98,12 @@ v); } -/* "Fix" the parameters of the given object to account for the +/* geom_fix_object_ptr is called after an object's externally-configurable parameters + have been initialized, but before any actual geometry calculations are done; + it is an opportunity to (re)compute internal data fields (such as cached + rotation matrices) that depend on externally-configurable parameters. + + One example: "Fix" the parameters of the given object to account for the geometry_lattice basis, which may be non-orthogonal. In particular, this means that the normalization of several unit vectors, such as the cylinder or block axes, needs to be changed. @@ -105,59 +111,52 @@ Unfortunately, we can't do this stuff at object-creation time in Guile, because the geometry_lattice variable may not have been assigned to its final value at that point. */ -void geom_fix_object(geometric_object o) +void geom_fix_object_ptr(geometric_object *o) { - switch(o.which_subclass) { + switch(o->which_subclass) { case GEOM CYLINDER: - lattice_normalize(&o.subclass.cylinder_data->axis); - if (o.subclass.cylinder_data->which_subclass == CYL WEDGE) { - vector3 a = o.subclass.cylinder_data->axis; - vector3 s = o.subclass.cylinder_data->subclass.wedge_data->wedge_start; + lattice_normalize(&o->subclass.cylinder_data->axis); + if (o->subclass.cylinder_data->which_subclass == CYL WEDGE) { + vector3 a = o->subclass.cylinder_data->axis; + vector3 s = o->subclass.cylinder_data->subclass.wedge_data->wedge_start; double p = vector3_dot(s, matrix3x3_vector3_mult(geometry_lattice.metric, a)); - o.subclass.cylinder_data->subclass.wedge_data->e1 = + o->subclass.cylinder_data->subclass.wedge_data->e1 = vector3_minus(s, vector3_scale(p, a)); - lattice_normalize(&o.subclass.cylinder_data->subclass.wedge_data->e1); - o.subclass.cylinder_data->subclass.wedge_data->e2 = + lattice_normalize(&o->subclass.cylinder_data->subclass.wedge_data->e1); + o->subclass.cylinder_data->subclass.wedge_data->e2 = cartesian_to_lattice( - vector3_cross(lattice_to_cartesian(o.subclass.cylinder_data->axis), - lattice_to_cartesian(o.subclass.cylinder_data->subclass.wedge_data->e1))); + vector3_cross(lattice_to_cartesian(o->subclass.cylinder_data->axis), + lattice_to_cartesian(o->subclass.cylinder_data->subclass.wedge_data->e1))); } break; case GEOM BLOCK: { matrix3x3 m; - lattice_normalize(&o.subclass.block_data->e1); - lattice_normalize(&o.subclass.block_data->e2); - lattice_normalize(&o.subclass.block_data->e3); - m.c0 = o.subclass.block_data->e1; - m.c1 = o.subclass.block_data->e2; - m.c2 = o.subclass.block_data->e3; - o.subclass.block_data->projection_matrix = matrix3x3_inverse(m); + lattice_normalize(&o->subclass.block_data->e1); + lattice_normalize(&o->subclass.block_data->e2); + lattice_normalize(&o->subclass.block_data->e3); + m.c0 = o->subclass.block_data->e1; + m.c1 = o->subclass.block_data->e2; + m.c2 = o->subclass.block_data->e3; + o->subclass.block_data->projection_matrix = matrix3x3_inverse(m); break; } case GEOM PRISM: { - // recompute centroid of prism floor polygon from prism center, - // which may have shifted since prism was initialized - prism *prsm = o.subclass.prism_data; - vector3 zhat = prsm->m_p2c.c2; - double height = prsm->height; - prsm->centroid = vector3_plus(o.center, vector3_scale(-0.5*height,zhat)); - break; + init_prism(o); + break; } case GEOM COMPOUND_GEOMETRIC_OBJECT: { int i; - int n = o.subclass.compound_geometric_object_data - ->component_objects.num_items; - geometric_object *os = o.subclass.compound_geometric_object_data - ->component_objects.items; + int n = o->subclass.compound_geometric_object_data->component_objects.num_items; + geometric_object *os = o->subclass.compound_geometric_object_data->component_objects.items; for (i = 0; i < n; ++i) { #if MATERIAL_TYPE_ABSTRACT if (os[i].material.which_subclass == MAT MATERIAL_TYPE_SELF) - material_type_copy(&o.material, &os[i].material); + material_type_copy(&o->material, &os[i].material); #endif - geom_fix_object(os[i]); + geom_fix_object_ptr(os + i); } break; } @@ -166,14 +165,25 @@ } } +// deprecated API — doesn't work for prisms +void geom_fix_object(geometric_object o) +{ + geom_fix_object_ptr(&o); +} + /* fix all objects in the geometry list as described in geom_fix_object, above */ -void geom_fix_objects0(geometric_object_list geometry) +void geom_fix_object_list(geometric_object_list geometry) { int index; for (index = 0; index < geometry.num_items; ++index) - geom_fix_object(geometry.items[index]); + geom_fix_object_ptr(geometry.items + index); +} + +void geom_fix_objects0(geometric_object_list geometry) +{ + geom_fix_object_list(geometry); } void geom_fix_objects(void) @@ -240,7 +250,7 @@ boolean CTLIO point_in_objectp(vector3 p, geometric_object o) { - geom_fix_object(o); + geom_fix_object_ptr(&o); return point_in_fixed_objectp(p, o); } @@ -318,6 +328,7 @@ return(a*a + b*b + c*c <= 1.0); } } + break; // never get here but silence compiler warning } case GEOM PRISM: { @@ -420,7 +431,7 @@ vector3 CTLIO normal_to_object(vector3 p, geometric_object o) { - geom_fix_object(o); + geom_fix_object_ptr(&o); return normal_to_fixed_object(p, o); } @@ -429,6 +440,7 @@ vector3 r = vector3_minus(p,o.center); switch (o.which_subclass) { + case GEOM CYLINDER: { vector3 rm = matrix3x3_vector3_mult(geometry_lattice.metric, r); @@ -449,7 +461,8 @@ return vector3_minus(r, vector3_scale(proj + prad * (o.subclass.cylinder_data->subclass.cone_data->radius2 - radius) / height, o.subclass.cylinder_data->axis)); else return vector3_minus(r, vector3_scale(proj, o.subclass.cylinder_data->axis)); - } + } // case GEOM CYLINDER + case GEOM BLOCK: { vector3 proj = @@ -467,8 +480,10 @@ return matrix3x3_row2(o.subclass.block_data->projection_matrix); else return matrix3x3_row3(o.subclass.block_data->projection_matrix); - } + } // case BLK BLOCK_SELF + case BLK ELLIPSOID: + default: { vector3 isa = o.subclass.block_data->subclass.ellipsoid_data->inverse_semi_axes; @@ -477,16 +492,20 @@ proj.z *= isa.z * isa.z; return matrix3x3_transpose_vector3_mult( o.subclass.block_data->projection_matrix, proj); - } - } // switch + } // case BLK ELLIPSOID + + } // switch (o.subclass.block_data->which_subclass) + } // case GEOM BLOCK + case GEOM PRISM: - { - return normal_to_prism(o.subclass.prism_data, p); - } + return normal_to_prism(o.subclass.prism_data, p); + default: - return r; - } + return r; + } // switch (o.which_subclass) + + return r; // never get here } /**************************************************************************/ @@ -549,7 +568,7 @@ boolean CTLIO point_in_periodic_objectp(vector3 p, geometric_object o) { - geom_fix_object(o); + geom_fix_object_ptr(&o); return point_in_periodic_fixed_objectp(p, o); } @@ -646,7 +665,7 @@ void CTLIO display_geometric_object_info(int indentby, geometric_object o) { - geom_fix_object(o); + geom_fix_object_ptr(&o); printf("%*s", indentby, ""); switch (o.which_subclass) { case GEOM CYLINDER: @@ -724,7 +743,7 @@ o.subclass.block_data->e3.z); break; case GEOM PRISM: - display_prism_info(indentby, o.subclass.prism_data); + display_prism_info(indentby, &o); break; case GEOM COMPOUND_GEOMETRIC_OBJECT: { @@ -882,8 +901,10 @@ ++ret; } return ret; - } - case BLK ELLIPSOID: + } // case BLK BLOCK_SELF: + + case BLK ELLIPSOID: + default: { vector3 isa = o.subclass.block_data->subclass.ellipsoid_data->inverse_semi_axes; double a, b2, c, discrim; @@ -905,10 +926,11 @@ s[1] = (b2 - discrim) / a; return 2; } - } - } - } // case GEOM BLOCK + } // case BLK BLOCK_SELF, default + } // switch (o.subclass.block_data->which_subclass) + + } // case GEOM BLOCK default: return 0; } @@ -1070,7 +1092,7 @@ etcetera. */ void geom_get_bounding_box(geometric_object o, geom_box *box) { - geom_fix_object(o); + geom_fix_object_ptr(&o); /* initialize to empty box at the center of the object: */ box->low = box->high = o.center; @@ -1996,7 +2018,7 @@ o.subclass.cylinder_data->height = height; o.subclass.cylinder_data->axis = axis; o.subclass.cylinder_data->which_subclass = CYL CYLINDER_SELF; - geom_fix_object(o); + geom_fix_object_ptr(&o); return o; } @@ -2022,7 +2044,7 @@ CHECK(o.subclass.cylinder_data->subclass.wedge_data, "out of memory"); o.subclass.cylinder_data->subclass.wedge_data->wedge_angle = wedge_angle; o.subclass.cylinder_data->subclass.wedge_data->wedge_start = wedge_start; - geom_fix_object(o); + geom_fix_object_ptr(&o); return o; } @@ -2050,7 +2072,7 @@ o.subclass.block_data->e3 = e3; o.subclass.block_data->size = size; o.subclass.block_data->which_subclass = BLK BLOCK_SELF; - geom_fix_object(o); + geom_fix_object_ptr(&o); return o; } @@ -2072,79 +2094,59 @@ } /*************************************************************** - * the remainder of this file is the content of `prism.c` - * implementing geometric primitives for prisms - * - * prism.c -- geometry routines for prisms. - * a prism is a planar polygon, consisting of 3 or more user-specified + * The remainder of this file implements geometric primitives for prisms. + * A prism is a planar polygon, consisting of 3 or more user-specified * vertices, extruded through a given thickness (the "height") in the - * direction of a given axis. - * most calculations are done in the "prism coordinate system", + * direction of a given unit vector (the "axis.") + * Most calculations are done in the "prism coordinate system", * in which the prism floor lies in the XY plane with centroid * at the origin and the prism axis is the positive Z-axis. + * Some variable naming conventions: + * -- Suffix 'p' or '_p' on variable names identifies variables + * storing coordinates or vector components in the prism system. + * Suffix 'c' or '_c' (or no suffix) corresponds to coodinates/components + * in ordinary 3d space. ('c' stands for 'cartesian'). + * -- We use the term 'vertex' for points in 3-space, stored as vector3 + * quantities with variable names beginning with 'p' or 'v'. For 3D + * direction vectors we use variable names beginning with 'd'. + * -- We use the term 'node' for points in 2-space, stored as vector3 + * quantities (with the z component unused) with variables beginning with 'q'. + * For 2D direction vectors we use variable names beginning with 'u'. * homer reid 4/2018 - * ***************************************************************/ /***************************************************************/ /* given coordinates of a point in the prism coordinate system,*/ /* return cartesian coordinates of that point */ /***************************************************************/ -vector3 prism_coordinate_p2c(prism *prsm, vector3 vp) -{ return vector3_plus(prsm->centroid, matrix3x3_vector3_mult(prsm->m_p2c,vp)); } +vector3 prism_coordinate_p2c(prism *prsm, vector3 pp) +{ return vector3_plus(prsm->centroid, matrix3x3_vector3_mult(prsm->m_p2c,pp)); } vector3 prism_vector_p2c(prism *prsm, vector3 vp) { return matrix3x3_vector3_mult(prsm->m_p2c, vp); } -vector3 prism_coordinate_c2p(prism *prsm, vector3 vc) -{ return matrix3x3_vector3_mult(prsm->m_c2p, vector3_minus(vc,prsm->centroid)); } +vector3 prism_coordinate_c2p(prism *prsm, vector3 pc) +{ return matrix3x3_vector3_mult(prsm->m_c2p, vector3_minus(pc,prsm->centroid)); } vector3 prism_vector_c2p(prism *prsm, vector3 vc) { return matrix3x3_vector3_mult(prsm->m_c2p, vc); } /***************************************************************/ -/***************************************************************/ -/***************************************************************/ -void get_prism_bounding_box(prism *prsm, geom_box *box) -{ - vector3 *vertices = prsm->vertices.items; - int num_vertices = prsm->vertices.num_items; - double height = prsm->height; - - box->low = box->high = prism_coordinate_p2c(prsm, vertices[0]); - int nv, fc; - for(nv=0; nv<num_vertices; nv++) - for(fc=0; fc<2; fc++) // 'floor,ceiling' - { vector3 vp = vertices[nv]; - if (fc==1) vp.z=height; - vector3 vc = prism_coordinate_p2c(prsm, vp); - - box->low.x = fmin(box->low.x, vc.x); - box->low.y = fmin(box->low.y, vc.y); - box->low.z = fmin(box->low.z, vc.z); - - box->high.x = fmax(box->high.x, vc.x); - box->high.y = fmax(box->high.y, vc.y); - box->high.z = fmax(box->high.z, vc.z); - } -} - -/***************************************************************/ -/* given 2D points p,v1,v2 and a 2D direction vector d, */ -/* determine whether or not the line p + s*d intersects the */ -/* line segment v1--v2. */ -/* algorithm: solve the 2x2 linear system p+s*d = a+t*b */ -/* where s,t are scalars and p,d,a,b are 2-vectors with */ -/* a=v1, b=v2-v1. intersection then corresponds to 0 <= t < 1. */ +/* given 2D points q0,q1,q2 and a 2D vector u, determine */ +/* whether or not the line q0 + s*u intersects the line */ +/* segment q1--q2. */ +/* algorithm: solve the 2x2 linear system q0+s*u = q1+t*(q2-q1)*/ +/* for the scalar quantities s, t; intersection corresponds to */ +/* 0 <= t < 1. */ /* return values: */ -/* ** case 1: d is not parallel to v1--v2 ** */ -/* NON_INTERSECTING test negative */ -/* INTERSECTING test positive */ -/* ** case 2: d is parallel to v1--v2 ** */ -/* IN_SEGMENT p lies on line segment v1--v2 */ -/* ON_RAY p does not lie on v1--v2, but there is a */ -/* *positive* value of s such that p+s*d */ -/* lies on v1--v2 */ +/* ** case 1: u is not parallel to q1--q2 ** */ +/* NON_INTERSECTING: test negative */ +/* INTERSECTING: test positive */ +/* ** case 2: u is parallel to q1--q2 ** */ +/* IN_SEGMENT: q0 lies on line segment q1--q2 */ +/* ON_RAY: q0 does not lie on q1--q2, but there is a*/ +/* *positive* value of s such that q0+s*u */ +/* lies on q1--q2 */ /* NON_INTERSECTING neither of the above */ /***************************************************************/ #define THRESH 1.0e-5 @@ -2152,29 +2154,29 @@ #define INTERSECTING 1 #define IN_SEGMENT 2 #define ON_RAY 3 -int intersect_line_with_segment(double px, double py, double dx, double dy, - vector3 v1, vector3 v2, double *s) +int intersect_line_with_segment(vector3 q0, vector3 q1, vector3 q2, vector3 u, double *s) { - double ax = v1.x, ay = v1.y; - double bx = v2.x-v1.x, by = v2.y-v1.y; - double M00 = dx, M10 = dy; - double M01 = -1.0*bx, M11 = -1.0*by; - double RHSx = ax - px, RHSy = ay - py; + /* ||ux q1x-q2x|| |s| = | q1x-q0x | */ + /* ||uy q1y-q2y|| |t| = | q1y-q0y | */ + double M00 = u.x, M01=q1.x-q2.x; + double M10 = u.y, M11=q1.y-q2.y; + double RHSx = q1.x-q0.x; + double RHSy = q1.y-q0.y; double DetM = M00*M11 - M01*M10; - double L2 = bx*bx + by*by; // squared length of edge + double L2 = M01*M01 + M11*M11; // squared length of edge, used to set length scale if ( fabs(DetM) < 1.0e-10*L2 ) { // d zero or nearly parallel to edge - double pmv1x = px-v1.x, pmv1y = py-v1.y, npmv1 = sqrt(pmv1x*pmv1x + pmv1y*pmv1y); - double pmv2x = px-v2.x, pmv2y = py-v2.y, npmv2 = sqrt(pmv2x*pmv2x + pmv2y*pmv2y); - double dot = pmv1x*pmv2x + pmv1y*pmv2y; - if ( fabs(dot) < (1.0-THRESH)*npmv1*npmv2 ) + double q01x = q0.x-q1.x, q01y = q0.y-q1.y, q01 = sqrt(q01x*q01x+q01y*q01y); + double q02x = q0.x-q2.x, q02y = q0.y-q2.y, q02 = sqrt(q02x*q02x+q02y*q02y); + double dot = q01x*q02x + q01y*q02y; + if ( fabs(dot) < (1.0-THRESH)*q01*q02 ) return NON_INTERSECTING; else if (dot<0.0) { *s=0.0; return IN_SEGMENT; } - else if ( (dx*pmv1x + dy*pmv1y) < 0.0 ) - { *s = fmin(npmv1, npmv2) / sqrt(dx*dx + dy*dy); + else if ( (u.x*q01x + u.y*q01y) < 0.0 ) + { *s = fmin(q01, q02) / sqrt(u.x*u.x + u.y*u.y); return ON_RAY; } return NON_INTERSECTING; @@ -2191,11 +2193,11 @@ } // like the previous routine, but only count intersections if s>=0 -boolean intersect_ray_with_segment(double px, double py, double dx, double dy, - vector3 v1, vector3 v2, double *s) +boolean intersect_ray_with_segment(vector3 q0, vector3 q1, vector3 q2, vector3 u, + double *s) { double ss; - int status=intersect_line_with_segment(px,py,dx,dy,v1,v2,&ss); + int status=intersect_line_with_segment(q0,q1,q2,u,&ss); if (status==INTERSECTING && ss<0.0) return NON_INTERSECTING; if (s) *s=ss; @@ -2203,58 +2205,59 @@ } /***************************************************************/ -/* 2D point-in-polygon test: return 1 if the point lies within */ -/* the polygon with the given vertices, 0 otherwise. */ +/* 2D point-in-polygon test: return 1 if q0 lies within the */ +/* polygon with the given vertices, 0 otherwise. */ // method: cast a plumb line in the negative y direction from */ -/* p to infinity and count the number of edges intersected; */ +/* q0 to infinity and count the number of edges intersected; */ /* point lies in polygon iff this is number is odd. */ /***************************************************************/ -boolean point_in_or_on_polygon(double px, double py, - vector3 *vertices, int num_vertices, - boolean include_boundaries) +boolean node_in_or_on_polygon(vector3 q0, vector3 *nodes, int num_nodes, + boolean include_boundaries) { - double dx=0.0, dy=-1.0; - int nv, edges_crossed=0; - for(nv=0; nv<num_vertices; nv++) - { int status = intersect_ray_with_segment(px, py, dx, dy, - vertices[nv], vertices[(nv+1)%num_vertices], 0); + vector3 u = {0.0, -1.0, 0.0}; + int nn, edges_crossed=0; + for(nn=0; nn<num_nodes; nn++) + { int status = intersect_ray_with_segment(q0, nodes[nn], nodes[(nn+1)%num_nodes], u, 0); if (status==IN_SEGMENT) return include_boundaries; - else if (status==INTERSECTING) + else if (status==INTERSECTING) edges_crossed++; else if (status==ON_RAY) - { vector3 vm1 = vertices[ (nv==0 ? num_vertices-1 : nv-1) ]; - vector3 v = vertices[ nv ]; - vector3 vp1 = vertices[ (nv+1) % num_vertices ]; - vector3 vp2 = vertices[ (nv+2) % num_vertices ]; - int last_status = intersect_ray_with_segment(px, py, dx, dy, vm1, v, 0); + { vector3 nm1 = nodes[ (nn==0 ? num_nodes-1 : nn-1) ]; + vector3 n0 = nodes[ nn ]; + vector3 np1 = nodes[ (nn+1) % num_nodes ]; + vector3 np2 = nodes[ (nn+2) % num_nodes ]; + int last_status = intersect_ray_with_segment(q0, nm1, n0, u, 0); if (last_status==INTERSECTING) edges_crossed--; - int next_status = intersect_ray_with_segment(px, py, dx, dy, vp1, vp2, 0); + int next_status = intersect_ray_with_segment(q0, np1, np2, u, 0); if (next_status==INTERSECTING) edges_crossed--; } } return (edges_crossed % 2); } -boolean point_in_polygon(double px, double py, - vector3 *vertices, int num_vertices) -{ return point_in_or_on_polygon(px, py, vertices, num_vertices, 1); } + +boolean node_in_polygon(double q0x, double q0y, vector3 *nodes, int num_nodes) +{ vector3 q0; + q0.x=q0x; q0.y=q0y; q0.z=0.0; + return node_in_or_on_polygon(q0, nodes, num_nodes, 1); +} /***************************************************************/ -/* return 1 or 0 if xc lies inside or outside the prism */ +/* return 1 or 0 if pc lies inside or outside the prism */ /***************************************************************/ -boolean point_in_or_on_prism(prism *prsm, vector3 xc, boolean include_boundaries) +boolean point_in_or_on_prism(prism *prsm, vector3 pc, boolean include_boundaries) { - vector3 *vertices = prsm->vertices.items; - int num_vertices = prsm->vertices.num_items; - double height = prsm->height; - vector3 xp = prism_coordinate_c2p(prsm, xc); - if ( xp.z<0.0 || xp.z>prsm->height) + double height = prsm->height; + vector3 pp = prism_coordinate_c2p(prsm, pc); + if ( pp.z<0.0 || pp.z>prsm->height ) return 0; - return point_in_or_on_polygon(xp.x, xp.y, vertices, num_vertices, include_boundaries); + vector3 *nodes = prsm->vertices_p.items; + int num_nodes = prsm->vertices_p.num_items; + return node_in_or_on_polygon(pp, nodes, num_nodes, include_boundaries); } -boolean point_in_prism(prism *prsm, vector3 xc) +boolean point_in_prism(prism *prsm, vector3 pc) { // by default, points on polygon edges are considered to lie inside the // polygon; this can be reversed by setting the environment variable @@ -2265,39 +2268,36 @@ char *s=getenv("LIBCTL_EXCLUDE_BOUNDARIES"); if (s && s[0]=='1') include_boundaries=0; } - - return point_in_or_on_prism(prsm, xc, include_boundaries); + return point_in_or_on_prism(prsm, pc, include_boundaries); } // comparator for qsort static int dcmp(const void *pd1, const void *pd2) -{ double d1=*((double *)pd1), d2=*((double *)pd2); +{ double d1=*((double *)pd1), d2=*((double *)pd2); return (d1<d2) ? -1.0 : (d1>d2) ? 1.0 : 0.0; } -/***************************************************************/ -/* compute all values of s at which the line p+s*d intersects */ -/* a prism face. */ -/* slist is a caller-allocated buffer with enough room for */ -/* at least num_vertices+2 doubles. on return it contains */ -/* the intersection s-values sorted in ascending order. */ -/* the return value is the number of intersections. */ -/***************************************************************/ -int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double *slist) -{ - int num_vertices = prsm->vertices.num_items; - vector3 *vertices = prsm->vertices.items; - double height = prsm->height; - - // get coords of p (line origin) and components of d (line slope) - // in prism coordinate system - vector3 p_prsm = prism_coordinate_c2p(prsm,p); - vector3 d_prsm = prism_vector_c2p(prsm,d); +/******************************************************************/ +/* 3D line-prism intersection: compute all values of s at which */ +/* the line p+s*d intersects a prism face. */ +/* pc, dc = cartesian coordinates of p, cartesian components of d */ +/* slist is a caller-allocated buffer with enough room for */ +/* at least num_vertices+2 doubles. on return it contains */ +/* the intersection s-values sorted in ascending order. */ +/* the return value is the number of intersections. */ +/******************************************************************/ +int intersect_line_with_prism(prism *prsm, vector3 pc, vector3 dc, double *slist) +{ + vector3 pp = prism_coordinate_c2p(prsm,pc); + vector3 dp = prism_vector_c2p(prsm,dc); + vector3 *vps = prsm->vertices_p.items; + int num_vertices = prsm->vertices_p.num_items; + double height = prsm->height; - // use length of first edge as a length scale for judging + // use length of first polygon edge as a length scale for judging // lengths to be small or large - double length_scale = vector3_norm(vector3_minus(vertices[1], vertices[0])); + double length_scale = vector3_norm(vector3_minus(vps[1], vps[0])); // identify intersections with prism side faces int num_intersections=0; @@ -2310,32 +2310,30 @@ // intersection of the XY-plane projection of the line with the // polygon edge between vertices (nv,nv+1). double s; - int status=intersect_line_with_segment(p_prsm.x, p_prsm.y, d_prsm.x, d_prsm.y, - vertices[nv], vertices[nvp1], &s); + int status=intersect_line_with_segment(pp, vps[nv], vps[nvp1], dp, &s); if (status==NON_INTERSECTING || status==ON_RAY) continue; // OK, we know the XY-plane projection of the line intersects the polygon edge; // now go back to 3D, ask for the z-coordinate at which the line intersects // the z-axis extrusion of the edge, and determine whether this point lies // inside or outside the height of the prism. - double z_intersect = p_prsm.z + s*d_prsm.z; - double AbsTol = 1.0e-7*length_scale; - double z_min = 0.0 - AbsTol; - double z_max = height + AbsTol; - if ( (z_intersect<z_min) || (z_intersect>z_max) ) + double zp_intersect = pp.z + s*dp.z; + double AbsTol = 1.0e-7*length_scale; + double zp_min = 0.0 - AbsTol; + double zp_max = height + AbsTol; + if ( (zp_intersect<zp_min) || (zp_intersect>zp_max) ) continue; slist[num_intersections++]=s; } // identify intersections with prism ceiling and floor faces - double dz = d_prsm.z; int LowerUpper; - if ( fabs(dz) > 1.0e-7*vector3_norm(d_prsm)) + if ( fabs(dp.z) > 1.0e-7*vector3_norm(dp) ) for(LowerUpper=0; LowerUpper<2; LowerUpper++) - { double z0 = LowerUpper ? height : 0.0; - double s = (z0 - p_prsm.z)/dz; - if (!point_in_polygon(p_prsm.x + s*d_prsm.x, p_prsm.y+s*d_prsm.y, vertices, num_vertices)) + { double z0p = LowerUpper ? height : 0.0; + double s = (z0p - pp.z)/dp.z; + if (!node_in_polygon(pp.x+s*dp.x, pp.y+s*dp.y, vps, num_vertices)) continue; slist[num_intersections++]=s; } @@ -2347,10 +2345,10 @@ /***************************************************************/ /***************************************************************/ /***************************************************************/ -double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b) +double intersect_line_segment_with_prism(prism *prsm, vector3 pc, vector3 dc, double a, double b) { double *slist=prsm->workspace.items; - int num_intersections=intersect_line_with_prism(prsm, p, d, slist); + int num_intersections=intersect_line_with_prism(prsm, pc, dc, slist); // na=smallest index such that slist[na] > a int na=-1; @@ -2431,7 +2429,7 @@ // the in-plane distance from pPlane to the quadrilateral double min_distance_to_quadrilateral(vector3 p, vector3 o, vector3 v1, vector3 v2, vector3 v3) -{ +{ int inside; double s=normal_distance_to_plane(p, o, v1, v2, v3, &inside); if(inside==1) @@ -2447,39 +2445,42 @@ return sqrt(s*s+d*d); } -double min_distance_to_prism_roof_or_ceiling(vector3 p, - vector3 *vertices, int num_vertices, - double height) -{ - vector3 o = {0.0,0.0,0.0}; o.z=height; - vector3 v3 = {0,0,1.0}; - double s = normal_distance_to_plane(p,o,vertices[0],vertices[1],v3,0); - vector3 pPlane = vector3_minus(p, vector3_scale(s,v3) ); - if (point_in_polygon(pPlane.x,pPlane.y,vertices,num_vertices)==1) +// fc==0/1 for floor/ceiling +double min_distance_to_prism_roof_or_ceiling(vector3 pp, prism *prsm, int fc) +{ + vector3 *vps = prsm->vertices_p.items; + int num_vertices = prsm->vertices_p.num_items; + vector3 op = {0.0,0.0,0.0}; if (fc==1) op.z = prsm->height; // origin of floor/ceiling + vector3 zhatp = {0,0,1.0}; + double s = normal_distance_to_plane(pp,op,vps[0],vps[1],zhatp,0); + vector3 ppProj = vector3_minus(pp, vector3_scale(s,zhatp) ); // projection of p into plane of floor/ceiling + if (node_in_polygon(ppProj.x,ppProj.y,vps,num_vertices)==1) return s; int nv; - double d=min_distance_to_line_segment(pPlane,vertices[0],vertices[1] ); + double d=min_distance_to_line_segment(ppProj,vps[0],vps[1] ); for(nv=1; nv<num_vertices; nv++) - d=fmin(d,min_distance_to_line_segment(pPlane,vertices[nv],vertices[(nv+1)%num_vertices])); + d=fmin(d,min_distance_to_line_segment(ppProj,vps[nv],vps[(nv+1)%num_vertices])); return sqrt(s*s+d*d); } /***************************************************************/ /* find the face of the prism for which the normal distance */ -/* from x to the plane of that face is the shortest, then */ +/* from p to the plane of that face is the shortest, then */ /* return the normal vector to that plane. */ /***************************************************************/ -vector3 normal_to_prism(prism *prsm, vector3 xc) +vector3 normal_to_prism(prism *prsm, vector3 pc) { - double height = prsm->height; - vector3 *vertices = prsm->vertices.items; - int num_vertices = prsm->vertices.num_items; + if (prsm->height==0.0) + return prsm->axis; - vector3 xp = prism_coordinate_c2p(prsm, xc); - vector3 zhat = {0,0,1.0}, axis=vector3_scale(height, zhat); - if (height==0.0) - return prism_vector_p2c(prsm, zhat); + double height = prsm->height; + vector3 *vps = prsm->vertices_p.items; + int num_vertices = prsm->vertices_p.num_items; + + vector3 zhatp = {0.0, 0.0, 1.0}; + vector3 axisp = vector3_scale(height, zhatp); + vector3 pp = prism_coordinate_c2p(prsm, pc); vector3 retval; double min_distance=HUGE_VAL; @@ -2488,58 +2489,76 @@ for(nv=0; nv<num_vertices; nv++) { int nvp1 = ( nv==(num_vertices-1) ? 0 : nv+1 ); - vector3 v0 = vertices[nv]; - vector3 v1 = vector3_minus(vertices[nvp1],vertices[nv]); - vector3 v2 = axis; - vector3 v3 = unit_vector3(vector3_cross(v1, v2)); - double s = min_distance_to_quadrilateral(xp, v0, v1, v2, v3); + vector3 v0p = vps[nv]; + vector3 v1p = vector3_minus(vps[nvp1],vps[nv]); + vector3 v2p = axisp; + vector3 v3p = unit_vector3(vector3_cross(v1p, v2p)); + double s = min_distance_to_quadrilateral(pp, v0p, v1p, v2p, v3p); if (fabs(s) < min_distance) { min_distance = fabs(s); - retval = v3; + retval = v3p; } } int fc; // 'floor / ceiling' for(fc=0; fc<2; fc++) - { double s=min_distance_to_prism_roof_or_ceiling(xp, vertices, num_vertices, - fc==0 ? 0.0 : height); + { double s=min_distance_to_prism_roof_or_ceiling(pp, prsm, fc); if (fabs(s) < min_distance) { min_distance = fabs(s); - retval = zhat; + retval = zhatp; } } return prism_vector_p2c(prsm, retval); } - /***************************************************************/ /***************************************************************/ /***************************************************************/ -void display_prism_info(int indentby, prism *prsm) +void get_prism_bounding_box(prism *prsm, geom_box *box) { vector3 *vertices = prsm->vertices.items; int num_vertices = prsm->vertices.num_items; - double height = prsm->height; - vector3 z0 = {0.0, 0.0, 1.0}; - vector3 axis = prism_vector_p2c(prsm,z0); + box->low = box->high = vertices[0]; + int nv, fc; + for(nv=0; nv<num_vertices; nv++) + for(fc=0; fc<2; fc++) // 'floor,ceiling' + { + vector3 v = vertices[nv]; + if (fc==1) + v = vector3_plus(v, vector3_scale(prsm->height, prsm->axis) ); + + box->low.x = fmin(box->low.x, v.x); + box->low.y = fmin(box->low.y, v.y); + box->low.z = fmin(box->low.z, v.z); + + box->high.x = fmax(box->high.x, v.x); + box->high.y = fmax(box->high.y, v.y); + box->high.z = fmax(box->high.z, v.z); + } +} + +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ +void display_prism_info(int indentby, geometric_object *o) +{ + prism *prsm = o->subclass.prism_data; + + vector3 *vs = prsm->vertices.items; + int num_vertices = prsm->vertices.num_items; printf("%*s height %g, axis (%g,%g,%g), %i vertices:\n", - indentby, "", height,axis.x,axis.y,axis.z,num_vertices); - matrix3x3 m_p2c = matrix3x3_inverse(prsm->m_c2p); + indentby, "",prsm->height,prsm->axis.x,prsm->axis.y,prsm->axis.z,num_vertices); int nv; for(nv=0; nv<num_vertices; nv++) - { vector3 v = matrix3x3_vector3_mult(m_p2c, vertices[nv]); - printf("%*s (%g,%g,%g)\n",indentby,"",v.x,v.y,v.z); - } + printf("%*s (%g,%g,%g)\n",indentby,"",vs[nv].x,vs[nv].y,vs[nv].z); } /***************************************************************/ // like vector3_equal but tolerant of small floating-point discrepancies /***************************************************************/ int vector3_nearly_equal(vector3 v1, vector3 v2, double tolerance) -{ - return (vector3_norm( vector3_minus(v1,v2) ) <= tolerance*vector3_norm(v1)); -} +{ return (vector3_norm( vector3_minus(v1,v2) ) <= tolerance*vector3_norm(v1)); } /***************************************************************/ /* return the unit normal vector to the triangle with the given*/ @@ -2554,32 +2573,82 @@ } /***************************************************************/ -/***************************************************************/ -/***************************************************************/ -geometric_object make_prism(material_type material, - const vector3 *vertices, int num_vertices, - double height, vector3 axis) +/* On entry, the only fields in o->prism that are assumed to */ +/* be initialized are: vertices, height, and (optionally) axis.*/ +/* If axis has not been initialized (i.e. it is set to its */ +/* default value, which is the zero vector) then the prism axis*/ +/* is automatically computed as the normal to the vertex plane.*/ +/* If o->center is equal to auto_center on entry, then it is */ +/* set to the prism center, as computed from the vertices, */ +/* axis, and height. Otherwise, the prism is rigidly translated*/ +/* to center it at the specified value of o->center. */ +/***************************************************************/ +// special vector3 that signifies 'no value specified' +vector3 auto_center = { NAN, NAN, NAN }; +void init_prism(geometric_object *o) { - CHECK(num_vertices>=3, "fewer than 3 vertices in make_prism"); + prism *prsm = o->subclass.prism_data; + vector3 *vertices = prsm->vertices.items; + int num_vertices = prsm->vertices.num_items; + CHECK(num_vertices>=3, "fewer than 3 vertices in init_prism"); // compute centroid of vertices vector3 centroid = {0.0, 0.0, 0.0}; int nv; for(nv=0; nv<num_vertices; nv++) centroid = vector3_plus(centroid, vertices[nv]); - centroid = vector3_scale(1.0/((double)num_vertices), centroid); + prsm->centroid = centroid = vector3_scale(1.0/((double)num_vertices), centroid); - // make sure all vertices lie in a plane normal to the given axis - vector3 zhat = unit_vector3(axis); + // make sure all vertices lie in a plane, i.e. that the normal + // vectors to all triangles (v_n, v_{n+1}, centroid) agree. + int plane_normal_set=0; + vector3 plane_normal; double tol=1.0e-6; for(nv=0; nv<num_vertices; nv++) { int nvp1 = (nv+1) % num_vertices; - vector3 zhatp = triangle_normal(centroid,vertices[nv],vertices[nvp1]); - boolean axis_normal - = ( vector3_nearly_equal(zhat, zhatp, tol) - || vector3_nearly_equal(zhat, vector3_scale(-1.0,zhatp), tol) + vector3 tri_normal = triangle_normal(centroid,vertices[nv],vertices[nvp1]); + if (vector3_norm(tri_normal)==0.0) // vertices collinear with centroid + continue; + if (!plane_normal_set) + { plane_normal=tri_normal; + plane_normal_set=1; + } + else + { boolean normals_agree + = ( vector3_nearly_equal(plane_normal, tri_normal, tol) + || vector3_nearly_equal(plane_normal, vector3_scale(-1.0,tri_normal), tol) + ); + CHECK(normals_agree, "non-coplanar vertices in init_prism"); + } + } + + // if no prism axis was specified, set the prism axis equal to the + // normal to the vertex plane. + // if a prism axis was specified, check that it agrees up to sign + // with the normal to the vertex plane. + if ( vector3_norm(prsm->axis) == 0.0 ) + prsm->axis = plane_normal; + else + { prsm->axis = unit_vector3(prsm->axis); + boolean axis_normal_to_plane + = ( vector3_nearly_equal(prsm->axis, plane_normal, tol) + || vector3_nearly_equal(prsm->axis, vector3_scale(-1.0,plane_normal), tol) ); - CHECK(axis_normal, "axis not normal to vertex plane in make_prism"); + CHECK(axis_normal_to_plane, "axis not normal to vertex plane in init_prism"); + } + + // set current_center=prism center as determined by vertices and height. + // if the center of the geometric object was left unspecified, + // set it to current_center; otherwise displace the entire prism + // so that it is centered at the specified center. + vector3 current_center = vector3_plus(centroid, vector3_scale(0.5*prsm->height,prsm->axis) ); + if (isnan(o->center.x) && isnan(o->center.y) && isnan(o->center.z)) // center == auto-center + o->center = current_center; + else + { vector3 shift = vector3_minus(o->center, current_center); + for(nv=0; nv<num_vertices; nv++) + vertices[nv]=vector3_plus(vertices[nv],shift); + centroid=vector3_plus(centroid,shift); } // compute rotation matrix that operates on a vector of cartesian coordinates @@ -2591,9 +2660,9 @@ // i.e. it is a point lying on the bottom surface of the prism. // This is the origin of coordinates in the prism system. // The *center* of the geometric object is the center of mass of the - // 3D prism. So center = centroid + 0.5*zHat. + // 3D prism. So center = centroid + 0.5*height*zHat. vector3 x0hat={1.0,0.0,0.0}, y0hat={0.0,1.0,0.0}, z0hat={0.0,0.0,1.0}; - vector3 xhat, yhat; + vector3 xhat, yhat, zhat=prsm->axis; if (vector3_nearly_equal(zhat, x0hat, tol)) { xhat=y0hat; yhat=z0hat; } else if (vector3_nearly_equal(zhat, y0hat, tol)) { xhat=z0hat; yhat=x0hat; } else if (vector3_nearly_equal(zhat, z0hat, tol)) { xhat=x0hat; yhat=y0hat; } @@ -2601,34 +2670,46 @@ { xhat = unit_vector3(vector3_minus(vertices[1],vertices[0])); yhat = unit_vector3(vector3_cross(zhat,xhat)); } - matrix3x3 m_p2c = {xhat, yhat, zhat }; - matrix3x3 m_c2p = matrix3x3_inverse(m_p2c); - - prism *prsm = MALLOC1(prism); - CHECK(prsm, "out of memory"); - prsm->centroid = centroid; - prsm->height = height; + matrix3x3 m_p2c = {xhat, yhat, zhat}; prsm->m_p2c = m_p2c; - prsm->m_c2p = m_c2p; + prsm->m_c2p = matrix3x3_inverse(m_p2c); - // note that the vertices stored in the prism_data structure - // are the vertices in the *prism* coordinate system, which means - // their z-coordinates are all zero and in principle need not be stored - prsm->vertices.num_items = num_vertices; - prsm->vertices.items = (vector3 *)malloc(num_vertices*sizeof(vector3)); + // compute vertices in prism coordinate system + prsm->vertices_p.num_items = num_vertices; + prsm->vertices_p.items = (vector3 *)malloc(num_vertices*sizeof(vector3)); for(nv=0; nv<num_vertices; nv++) - prsm->vertices.items[nv] = prism_coordinate_c2p(prsm,vertices[nv]); + prsm->vertices_p.items[nv] = prism_coordinate_c2p(prsm,vertices[nv]); // workspace is an internally-stored double-valued array of length num_vertices+2 // that is used by some geometry routines prsm->workspace.num_items = num_vertices+2; prsm->workspace.items = (double *)malloc( (num_vertices+2)*sizeof(double) ); +} + +/***************************************************************/ +/* routines called from C++ or python codes to create prisms */ +/***************************************************************/ +// prism with center determined automatically from vertices, height, and axis +geometric_object make_prism(material_type material, + const vector3 *vertices, int num_vertices, + double height, vector3 axis) +{ return make_prism_with_center(material, auto_center, vertices, num_vertices, height, axis); } - // note the center and centroid are different! - vector3 center = vector3_plus(centroid, vector3_scale(0.5*height,zhat) ); +// prism in which all vertices are translated to ensure that the prism is centered at center +geometric_object make_prism_with_center(material_type material, vector3 center, + const vector3 *vertices, int num_vertices, + double height, vector3 axis) +{ geometric_object o=make_geometric_object(material, center); o.which_subclass=GEOM PRISM; - o.subclass.prism_data = prsm; - + prism *prsm = o.subclass.prism_data = MALLOC1(prism); + CHECK(prsm, "out of memory"); + prsm->vertices.num_items = num_vertices; + prsm->vertices.items = (vector3 *)malloc(num_vertices*sizeof(vector3)); + CHECK(prsm->vertices.items, "out of memory"); + memcpy(prsm->vertices.items, vertices, num_vertices*sizeof(vector3)); + prsm->height = height; + prsm->axis = axis; + init_prism(&o); return o; } diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/libctl-4.1.4/utils/geom.scm new/libctl-4.2.0/utils/geom.scm --- old/libctl-4.1.4/utils/geom.scm 2018-07-27 20:46:20.000000000 +0200 +++ new/libctl-4.2.0/utils/geom.scm 2019-01-07 20:32:57.000000000 +0100 @@ -41,6 +41,7 @@ (define-property material nothing MATERIAL-TYPE) (define-property center no-default 'vector3)) + (define-class compound-geometric-object geometric-object (define-property component-objects '() (make-list-type 'geometric-object))) @@ -85,14 +86,45 @@ (object-property-value object 'e2) (object-property-value object 'e3)))))) -(define identity_matrix (matrix3x3 (vector3 1 0 0) - (vector3 0 1 0) +(define identity_matrix (matrix3x3 (vector3 1 0 0) + (vector3 0 1 0) (vector3 0 0 1))) +; some notes regarding prisms: +; (a) When instantiating a prism, typically only the +; fields `vertices`, `height,` and (optionally) `axis` +; will be initialized by the user; all remaining fields are +; derived properties that are computed internally. (So, morally +; they should be thought of as having been declared using +; `define-derived-property` or `define-post-processed-property,` +; except here the code that does the derivation or +; post-processing is implemented in C, not scheme.) +; (b) The suffix _p (for "prism") is used to identify variables +; that store coordinates of points or components of vectors +; in the prism coordinate system. (The prism coordinate system +; is defined by the condition that the prism axis is the z-axis +; and the prism floor lies in the xy plane at z==0.) Variables +; with no suffix refer to quantities in ordinary 3D space. +; (c) "centroid" refers to the centroid of the prism floor polygon; this is +; the origin of the prism coordinate system [i.e. by definition +; we have centroid_p=(0 0 0)]. +; (d) If 'axis' is left unspecified, it is inferred to be the +; normal to the plane of the prism floor, with sign defined +; by a right-hand-rule with respect to the first two vertices, i.e. +; axis = normal_vector( (v1-centroid) x (v2-centroid) ) +; (e) The specification of the prism vertices and height suffices to +; determine the center of the geometric object +; (center = centroid + 0.5*height*axis), so---in contrast to all other +; types of geometric-object---there is no need to specify the `center` +; field when instantiating a prism. (define-class prism geometric-object +; fields to be filled in by users (define-property vertices '() (make-list-type 'vector3)) - (define-property centroid (vector3 0 0 0) 'vector3) (define-property height 0 'number) + (define-property axis (vector3 0 0 0) 'vector3) +; derived fields computed internally + (define-property vertices_p '() (make-list-type 'vector3)) + (define-property centroid (vector3 0 0 0) 'vector3) (define-property workspace '() (make-list-type 'number)) (define-property m_c2p identity_matrix 'matrix3x3) (define-property m_p2c identity_matrix 'matrix3x3)) @@ -104,7 +136,6 @@ (object-property-value object 'size))))) ; **************************************************************** - (define-class lattice no-parent (define-post-processed-property basis1 (vector3 1 0 0) 'vector3 unit-vector3) (define-post-processed-property basis2 (vector3 0 1 0) 'vector3 unit-vector3) @@ -200,6 +231,9 @@ (define-input-var geometry '() (make-list-type 'geometric-object)) (define-input-var ensure-periodicity true 'boolean) +; special vector3 that signifies 'no value specified' +(define auto-center (vector3 (nan) (nan) (nan))) + (define-external-function point-in-object? true false 'boolean 'vector3 'geometric-object) diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/libctl-4.1.4/utils/geomtst.c new/libctl-4.2.0/utils/geomtst.c --- old/libctl-4.1.4/utils/geomtst.c 2018-01-03 19:27:47.000000000 +0100 +++ new/libctl-4.2.0/utils/geomtst.c 2019-01-07 20:42:16.000000000 +0100 @@ -183,7 +183,7 @@ geometry_lattice.basis2 = random_unit_vector3(); geometry_lattice.basis3 = random_unit_vector3(); geom_fix_lattice(); - geom_fix_object(o); + geom_fix_object_ptr(&o); #endif b.low = make_vector3(myurand(-1,0), myurand(-1,0), myurand(-1,0)); diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/libctl-4.1.4/utils/test-prism.c new/libctl-4.2.0/utils/test-prism.c --- old/libctl-4.1.4/utils/test-prism.c 2018-09-07 20:51:38.000000000 +0200 +++ new/libctl-4.2.0/utils/test-prism.c 2019-01-07 20:32:57.000000000 +0100 @@ -81,8 +81,8 @@ /***************************************************************/ void my_get_prism_bounding_box(prism *prsm, geom_box *box) { - vector3 *vertices = prsm->vertices.items; - int num_vertices = prsm->vertices.num_items; + vector3 *vertices = prsm->vertices_p.items; + int num_vertices = prsm->vertices_p.num_items; double height = prsm->height; box->low = box->high = prism_coordinate_p2c(prsm, vertices[0]); @@ -155,8 +155,8 @@ vector3 random_point_on_prism(geometric_object o) { prism *prsm = o.subclass.prism_data; - vector3 *vertices = prsm->vertices.items; - int num_vertices = prsm->vertices.num_items; + vector3 *vertices = prsm->vertices_p.items; + int num_vertices = prsm->vertices_p.num_items; double height = prsm->height; // choose a face @@ -195,8 +195,8 @@ /***************************************************************/ void prism2gnuplot(prism *prsm, char *filename) { - vector3 *vertices = prsm->vertices.items; - int num_vertices = prsm->vertices.num_items; + vector3 *vertices = prsm->vertices_p.items; + int num_vertices = prsm->vertices_p.num_items; double height = prsm->height; FILE *f=fopen(filename,"w"); @@ -227,8 +227,8 @@ /***************************************************************/ void prism2gmsh(prism *prsm, char *filename) { - vector3 *vertices = prsm->vertices.items; - int num_vertices = prsm->vertices.num_items; + vector3 *vertices = prsm->vertices_p.items; + int num_vertices = prsm->vertices_p.num_items; double height = prsm->height; vector3 zhat = prsm->m_p2c.c2; vector3 axis = vector3_scale(height, zhat); @@ -289,10 +289,8 @@ num_adjusted++; } - if (in_block!=in_prism) - { num_failed++; - if (f) fprintf(f,"%i %i %e %e %e \n",in_block,in_prism,p.x,p.y,p.z); - } + if (in_block!=in_prism) num_failed++; + if (f) fprintf(f,"%i %i %e %e %e \n",in_block,in_prism,p.x,p.y,p.z); } if (f) fclose(f);
