From 37afb5bf2712f581b046309945f49a40b016ef0d Mon Sep 17 00:00:00 2001 From: Homer Reid Date: Thu, 1 Mar 2018 20:00:14 -0500 Subject: [PATCH 01/10] add support for prisms (new type of geometric_object) to libctlgeom --- examples/Makefile.am | 5 +- utils/Makefile.am | 7 +- utils/ctlgeom.h | 11 ++ utils/geom.c | 60 ++++++- utils/geom.scm | 9 ++ utils/prism.c | 373 +++++++++++++++++++++++++++++++++++++++++++ utils/test-prism.c | 142 ++++++++++++++++ 7 files changed, 599 insertions(+), 8 deletions(-) create mode 100644 utils/prism.c create mode 100644 utils/test-prism.c diff --git a/examples/Makefile.am b/examples/Makefile.am index 1f6bb8f..8363f92 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -57,7 +57,7 @@ CTL_DEFS = -DCTL_SCM='"'$(LIBCTL_DIR)/base/ctl.scm'"' \ -DVERSION_STRING='"'$(VERSION_STRING)'"' example_SOURCES = $(MY_SOURCES) -nodist_example_SOURCES = main.c geom.c ctl-io.h ctl-io.c +nodist_example_SOURCES = main.c geom.c prism.c ctl-io.h ctl-io.c BUILT_SOURCES = $(nodist_example_SOURCES) example_LDADD = $(MY_LIBS) $(LIBCTL) example_LDFLAGS = $(MY_LDFLAGS) @@ -69,6 +69,9 @@ main.c: $(LIBCTL_DIR)/base/main.c geom.c: $(LIBCTL_DIR)/utils/geom.c cp -f $(LIBCTL_DIR)/utils/geom.c $@ +prism.c: $(LIBCTL_DIR)/utils/prism.c + cp -f $(LIBCTL_DIR)/utils/prism.c $@ + ctl-io.c: $(SPECIFICATION_FILE) $(LIBCTL_DIR)/utils/geom.scm $(GEN_CTL_IO) $(GEN_CTL_IO) --code -o $@ $(SPECIFICATION_FILE) $(LIBCTL_DIR) diff --git a/utils/Makefile.am b/utils/Makefile.am index 8d05579..8fa2acf 100644 --- a/utils/Makefile.am +++ b/utils/Makefile.am @@ -6,7 +6,7 @@ noinst_PROGRAMS = geomtst EXTRA_DIST = gen-ctl-io.in README geom.scm geom-ctl-io-defaults.c nlopt.c -libctlgeom_la_SOURCES = geom.c $(top_srcdir)/src/ctl-math.c $(top_srcdir)/src/integrator.c geom-ctl-io.c ctlgeom-types.h +libctlgeom_la_SOURCES = geom.c $(top_srcdir)/src/ctl-math.c $(top_srcdir)/src/integrator.c geom-ctl-io.c prism.c ctlgeom-types.h libctlgeom_la_LDFLAGS = -version-info @SHARED_VERSION_INFO@ libctlgeom_la_CPPFLAGS = -DLIBCTLGEOM @@ -39,5 +39,10 @@ geom-ctl-io.c: ctl-io.c endif +check_PROGRAMS = test-prism + +test_prism_SOURCES = test-prism.c +test_prism_LDADD = libctlgeom.la + clean-local: rm -f ctl-io.[ch] nlopt-constants.scm diff --git a/utils/ctlgeom.h b/utils/ctlgeom.h index f1f8f85..b6b8c05 100644 --- a/utils/ctlgeom.h +++ b/utils/ctlgeom.h @@ -145,6 +145,17 @@ GEOMETRIC_OBJECT make_ellipsoid(MATERIAL_TYPE material, vector3 center, vector3 e1, vector3 e2, vector3 e3, vector3 size); +// routines in prism.c +GEOMETRIC_OBJECT make_prism(MATERIAL_TYPE material, + vector3 *vertices, int num_vertices, + double height, vector3 axis); + +void get_prism_bounding_box(prism *prsm, geom_box *box); +int intersect_line_with_prism(vector3 p, vector3 d, prism *prsm, double **sarray); +boolean point_in_prism(vector3 p, prism *prsm); +vector3 normal_to_prism(vector3 p, prism *prsm); +void display_prism_info(int indentby, prism *prsm); + /**************************************************************************/ diff --git a/utils/geom.c b/utils/geom.c index 036403f..31471e1 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -21,6 +21,7 @@ #include #include +#include #include #ifndef LIBCTLGEOM @@ -59,6 +60,7 @@ using namespace ctlio; #endif #define K_PI 3.14159265358979323846 +#define CHECK(cond, s) if (!(cond)){fprintf(stderr,s "\n");exit(EXIT_FAILURE);} /**************************************************************************/ @@ -123,6 +125,10 @@ void geom_fix_object(geometric_object o) o.subclass.block_data->projection_matrix = matrix3x3_inverse(m); break; } + case GEOM PRISM: + { + // FIXME not really sure what to do here + } case GEOM COMPOUND_GEOMETRIC_OBJECT: { int i; @@ -205,6 +211,19 @@ void geom_initialize(void) geometry.items = 0; } +/**************************************************************************/ +/* like intersect_line_with_object, but allows for arbitrarily many */ +/* intersections (not just 2) to accommodate non-star-shaped objects. */ +/* if svalues is nonzero, it points on return to a newly allocated buffer */ +/* that must be free()d by the caller. */ +/**************************************************************************/ +#if HAVE_PRISM +int intersect_line_with_prism(vector3 p, vector3 d, prism *prsm, + double **svalues) +{ +} +#endif + /**************************************************************************/ /* Return whether or not the point p (in the lattice basis) is inside @@ -297,6 +316,10 @@ boolean point_in_fixed_pobjectp(vector3 p, geometric_object *o) } } } + case GEOM PRISM: + { + return point_in_prism(r, o->subclass.prism_data); + } case GEOM COMPOUND_GEOMETRIC_OBJECT: { int i; @@ -348,6 +371,8 @@ vector3 to_geom_object_coords(vector3 p, geometric_object o) if (size.z != 0.0) proj.z /= size.z; return vector3_plus(half, proj); } +/* case GEOM PRISM: + NOT YET IMPLEMENTED */ } } @@ -378,6 +403,8 @@ vector3 from_geom_object_coords(vector3 p, geometric_object o) vector3_scale(size.z * p.z, o.subclass.block_data->e3)) )); } +/* case GEOM PRISM: + NOT YET IMPLEMENTED */ } } @@ -448,7 +475,11 @@ vector3 normal_to_fixed_object(vector3 p, geometric_object o) return matrix3x3_transpose_vector3_mult( o.subclass.block_data->projection_matrix, proj); } - } + } // switch + } // case GEOM BLOCK + case GEOM PRISM: + { + return normal_to_prism(r, o.subclass.prism_data); } default: return r; @@ -641,6 +672,9 @@ void CTLIO display_geometric_object_info(int indentby, geometric_object o) break; } break; + case GEOM PRISM: + printf("prism"); + break; case GEOM COMPOUND_GEOMETRIC_OBJECT: printf("compound object"); break; @@ -686,6 +720,9 @@ void CTLIO display_geometric_object_info(int indentby, geometric_object o) o.subclass.block_data->e3.y, o.subclass.block_data->e3.z); break; + case GEOM PRISM: + display_prism_info(indentby, o.subclass.prism_data); + break; case GEOM COMPOUND_GEOMETRIC_OBJECT: { int i; @@ -734,7 +771,7 @@ int intersect_line_with_object(vector3 p, vector3 d, geometric_object o, s[1] = (b2 - discrim) / a; return 2; } - } + } // case GEOM SPHERE case GEOM CYLINDER: { vector3 dm = matrix3x3_vector3_mult(geometry_lattice.metric, d); vector3 pm = matrix3x3_vector3_mult(geometry_lattice.metric, p); @@ -796,7 +833,7 @@ int intersect_line_with_object(vector3 p, vector3 d, geometric_object o, } if (ret == 2 && s[0] == s[1]) ret = 1; return ret; - } + } // case GEOM CYLINDER case GEOM BLOCK: { vector3 dproj = matrix3x3_vector3_mult(o.subclass.block_data->projection_matrix, d); @@ -867,7 +904,15 @@ int intersect_line_with_object(vector3 p, vector3 d, geometric_object o, } } } - } + } // case GEOM BLOCK + case GEOM PRISM: + { double *sprism=0; + int n=intersect_line_with_prism(p, d, o.subclass.prism_data, &sprism); + CHECK( n<=2, "more than 2 intersections in intersect_line_with_object"); + memcpy(s, sprism, n*sizeof(double)); + if (n>0) free(sprism); + return n; + } // case GEOM PRISM default: return 0; } @@ -1140,6 +1185,11 @@ void geom_get_bounding_box(geometric_object o, geom_box *box) vector3_plus(corner, vector3_plus(s1, vector3_plus(s2, s3)))); break; } + case GEOM PRISM: + get_prism_bounding_box(o.subclass.prism_data, box); + box->low=vector3_plus(box->low,o.center); + box->high=vector3_plus(box->high,o.center); + break; case GEOM COMPOUND_GEOMETRIC_OBJECT: { int i; @@ -1411,8 +1461,6 @@ static int object_in_box(geometric_object o, vector3 shiftby, return geom_boxes_intersect(obj_b, b); } -#define CHECK(cond, s) if (!(cond)){fprintf(stderr,s "\n");exit(EXIT_FAILURE);} - static geom_box_tree new_geom_box_tree(void) { geom_box_tree t; diff --git a/utils/geom.scm b/utils/geom.scm index ee0325a..10d22ae 100644 --- a/utils/geom.scm +++ b/utils/geom.scm @@ -85,6 +85,15 @@ (object-property-value object 'e2) (object-property-value object 'e3)))))) +(define identity_matrix (matrix3x3 (vector3 1 0 0) + (vector3 0 1 0) + (vector3 0 0 1))) + +(define-class prism geometric-object + (define-property vertices '() (make-list-type 'vector3)) + (define-property height 0 'number) + (define-property m_c2p identity_matrix 'matrix3x3)) + (define-class ellipsoid block (define-derived-property inverse-semi-axes 'vector3 (lambda (object) diff --git a/utils/prism.c b/utils/prism.c new file mode 100644 index 0000000..083e2ab --- /dev/null +++ b/utils/prism.c @@ -0,0 +1,373 @@ +/* libctl: flexible Guile-based control files for scientific software + * Copyright (C) 1998-2014 Massachusetts Institute of Technology and Steven G. Johnson + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + * + * Steven G. Johnson can be contacted at stevenj@alum.mit.edu. + */ + +/* + * prism.c -- geometry routines for prisms + * a prism is a planar polygon, consisting of any + * number of user-specified vertices, extruded + * through a given thickness (the "height") in the + * direction of a given axis. + */ + +#include +#include +#include + +#ifndef LIBCTLGEOM +# include "ctl-io.h" +#else +# define material_type void* + static void material_type_copy(void **src, void **dest) { *dest = *src; } +#endif +#include + +#ifdef CXX_CTL_IO +using namespace ctlio; +# define CTLIO ctlio:: +# define GEOM geometric_object:: +# define BLK block:: +# define CYL cylinder:: +# define MAT material_type:: +#else +# define CTLIO +# define GEOM +# define BLK +# define CYL +# define MAT +#endif + +#ifdef __cplusplus +# define MALLOC(type,num) (new type[num]) +# define MALLOC1(type) (new type) +# define FREE(p) delete[] (p) +# define FREE1(p) delete (p) +#else +# define MALLOC(type,num) ((type *) malloc(sizeof(type) * (num))) +# define MALLOC1(type) MALLOC(type,1) +# define FREE(p) free(p) +# define FREE1(p) free(p) +#endif + +#define K_PI 3.14159265358979323846 +#define CHECK(cond, s) if (!(cond)){fprintf(stderr,s "\n");exit(EXIT_FAILURE);} + +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ +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 = matrix3x3_vector3_mult(prsm->m_c2p,z0); + + // set x,y coordinates of low, high to bounding box + // of prism base polygon (in prism coordinate system) + vector3 low = vertices[0]; + vector3 high = vertices[0]; + for(int nv=1; nvm_c2p); + box->low = matrix3x3_vector3_mult(m_p2c, low); + box->high = matrix3x3_vector3_mult(m_p2c, high); +} + +/***************************************************************/ +/* find the value of s at which the line p+s*d intersects the */ +/* line segment connecting v1 to v2 (in 2 dimensions) */ +/* 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 */ +/***************************************************************/ +boolean intersect_line_with_segment(double px, double py, double dx, double dy, + vector3 v1, vector3 v2, 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; + double DetM = M00*M11 - M01*M10; + double L2 = bx*bx + by*by; // squared length of edge + if ( fabs(DetM) < 1.0e-10*L2 ) // d zero or nearly parallel to edge-->no intersection + return 0; + + double t = (M00*RHSy-M10*RHSx)/DetM; + if (t<0.0 || t>1.0) // intersection of lines does not lie between vertices + return 0; + + if (s) *s = (M11*RHSx-M01*RHSy)/DetM; + return 1; +} + +// 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) +{ + double ss; + int status=intersect_line_with_segment(px,py,dx,dy,v1,v2,&ss); + if (status==0 && ss<0.0) + return 0; + if (s) *s=ss; + return 1; +} + +/***************************************************************/ +/* 2D point-in-polygon test: return 1 if the point 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; */ +/* point lies in polygon iff this is number is odd. */ +/***************************************************************/ +boolean point_in_polygon(double px, double py, vector3 *vertices, int num_vertices) +{ + double dx=0.0, dy=-1.0; + int num_side_intersections=0; + for(int nv=0; nvvertices.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 + matrix3x3 m_c2p = prsm->m_c2p; + vector3 p_prsm = matrix3x3_vector3_mult(m_c2p, p); + vector3 d_prsm = matrix3x3_vector3_mult(m_c2p, d); + + int num_intersections=0; + double *svalues=0; + + // identify intersections of line with all prism side surfaces + for(int nv=0; nv0.0 && (fabs(z_intersect) > 0.5*height ) ) + continue; + + num_intersections++; + if (sarray) + { svalues = (double *)realloc(svalues, num_intersections*sizeof(double)); + svalues[num_intersections-1] = s; + }; + }; + + // identify intersections of line with prism top and bottom surfaces + double dz=d_prsm.z; + if (dz > 1.0e-10*vector3_norm(d_prsm)) + for(double sign=-1.0; sign<=1.01; sign+=2.0) + { double s = (p_prsm.z - sign*0.5*height)/dz; + if (!point_in_polygon(p_prsm.x + s*d_prsm.x, p_prsm.y+s*d_prsm.y, vertices, num_vertices)) + continue; + num_intersections++; + if (sarray) + { svalues = (double *)realloc(svalues, num_intersections*sizeof(double)); + svalues[num_intersections-1] = s; + }; + }; + + if (sarray) *sarray=svalues; + return num_intersections; + +} + +vector3 normal_to_prism(vector3 p, prism *prsm) +{ +} + +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ +void display_prism_info(int indentby, prism *prsm) +{ + 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 = matrix3x3_vector3_mult(prsm->m_c2p,z0); + + 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); + for(int nv=0; nvvertices.items; + int num_vertices = prsm->vertices.num_items; + double height = prsm->height; + vector3 p_prsm = matrix3x3_vector3_mult(prsm->m_c2p, p); + if ( fabs(p_prsm.z) > 0.5*prsm->height) + return 0; + return point_in_polygon(p_prsm.x, p_prsm.y, vertices, num_vertices); +} + +/***************************************************************/ +// like vector3_equal but tolerant of small floating-point discrepancies +/***************************************************************/ +int vector3_nearly_equal(vector3 v1, vector3 v2) +{ + return (vector3_norm( vector3_minus(v1,v2) ) < 1.0e-7*vector3_norm(v1)); +} + +/***************************************************************/ +/* return the unit normal vector to the triangle with the given*/ +/* vertices */ +/***************************************************************/ +vector3 triangle_normal(vector3 v1, vector3 v2, vector3 v3) +{ + vector3 nv=vector3_cross( vector3_minus(v2,v1), vector3_minus(v3,v1) ); + double nvnorm=vector3_norm(nv); + // if (area) *area += 0.5*nvnorm; + return unit_vector3(nv); +} + +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ +GEOMETRIC_OBJECT make_prism(MATERIAL_TYPE material, + vector3 *vertices, int num_vertices, + double height, vector3 axis) +{ + CHECK(num_vertices>=3, "fewer than 3 vertices in make_prism"); + + // compute centroid of vertices + vector3 centroid = {0.0, 0.0, 0.0}; + for(int nv=0; nvheight = height; + prsm->m_c2p = m_c2p; + + // 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 identical and in principle need not + // be stored + prsm->vertices.num_items = num_vertices; + prsm->vertices.items = (vector3 *)malloc(num_vertices*sizeof(vector3)); + for(int nv=0; nvvertices.items[nv] + =matrix3x3_vector3_mult(m_c2p, vector3_minus(vertices[nv],centroid)); + + vector3 center = vector3_plus(centroid, vector3_scale(0.5*height,zhat) ); + geometric_object o=make_geometric_object(material, center); + o.which_subclass=GEOM PRISM; + o.subclass.prism_data = prsm; + + return o; +} diff --git a/utils/test-prism.c b/utils/test-prism.c new file mode 100644 index 0000000..57de399 --- /dev/null +++ b/utils/test-prism.c @@ -0,0 +1,142 @@ +/************************************************************************/ +/* test-prism.c: unit test for prisms in libctlgeom */ +/************************************************************************/ +#include +#include +#include +#include + +#include "ctlgeom.h" + +#define K_PI 3.141592653589793238462643383279502884197 + +static vector3 make_vector3(double x, double y, double z) +{ + vector3 v; + v.x = x; v.y = y; v.z = z; + return v; +} + +/************************************************************************/ +/* return a uniform random number in [a,b] */ +/************************************************************************/ +static double myurand(double a, double b) +{ + return a + (b-a)*(rand()/((double)RAND_MAX)); +} + +static double mydrand() +{ return myurand(0.0,1.0); } + +/************************************************************************/ +/* return a random unit vector, uniformly distributed over a parallelepiped */ +/************************************************************************/ +vector3 random_point(vector3 min_corner, vector3 max_corner) +{ return make_vector3( myurand(min_corner.x, max_corner.x), + myurand(min_corner.y, max_corner.y), + myurand(min_corner.z, max_corner.z) + ); +} + +/************************************************************************/ +/* return a random unit vector, uniformly distributed over a sphere */ +/************************************************************************/ +vector3 random_unit_vector3(void) +{ + double z, t, r; + vector3 v; + + z = 2*mydrand() - 1; + t = 2*K_PI*mydrand(); + r = sqrt(1 - z*z); + v.x = r * cos(t); + v.y = r * sin(t); + v.z = z; + return v; +} + +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ +void draw_prism(geometric_object o, char *filename) +{ + vector3 center = o.center; + prism *prsm = o.subclass.prism_data; + vector3 *vertices = prsm->vertices.items; + int num_vertices = prsm->vertices.num_items; + double height = prsm->height; + matrix3x3 m_c2p = prsm->m_c2p; + matrix3x3 m_p2c = matrix3x3_inverse(m_c2p); + + FILE *f=fopen(filename,"w"); + for(int nv=0; nv Date: Sun, 4 Mar 2018 05:29:36 -0500 Subject: [PATCH 02/10] updates --- utils/prism.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/prism.c b/utils/prism.c index 083e2ab..9a8fdb9 100644 --- a/utils/prism.c +++ b/utils/prism.c @@ -230,7 +230,7 @@ int intersect_line_with_prism(vector3 p, vector3 d, prism *prsm, double **sarray } vector3 normal_to_prism(vector3 p, prism *prsm) -{ +{ // not implemented yet } /***************************************************************/ From 35dcb8799288eb0aa9a7f918087e2bf8f5e1bbe2 Mon Sep 17 00:00:00 2001 From: Homer Reid Date: Fri, 16 Mar 2018 22:39:02 +0900 Subject: [PATCH 03/10] updates to point_in_prism and other functions --- utils/ctlgeom.h | 10 +- utils/geom.c | 24 +--- utils/geom.scm | 4 +- utils/prism.c | 138 +++++++++--------- utils/test-prism.c | 341 +++++++++++++++++++++++++++++++++++++++++---- 5 files changed, 404 insertions(+), 113 deletions(-) diff --git a/utils/ctlgeom.h b/utils/ctlgeom.h index b6b8c05..8acc197 100644 --- a/utils/ctlgeom.h +++ b/utils/ctlgeom.h @@ -150,10 +150,14 @@ GEOMETRIC_OBJECT make_prism(MATERIAL_TYPE material, vector3 *vertices, int num_vertices, double height, vector3 axis); +vector3 prism_coordinate_p2c(prism *prsm, vector3 vp); +vector3 prism_coordinate_c2p(prism *prsm, vector3 vc); +vector3 prism_vector_p2c(prism *prsm, vector3 vp); +vector3 prism_vector_c2p(prism *prsm, vector3 vc); void get_prism_bounding_box(prism *prsm, geom_box *box); -int intersect_line_with_prism(vector3 p, vector3 d, prism *prsm, double **sarray); -boolean point_in_prism(vector3 p, prism *prsm); -vector3 normal_to_prism(vector3 p, prism *prsm); +int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double **sarray); +boolean point_in_prism(prism *prsm, vector3 p); +vector3 normal_to_prism(prism *prsm, vector3 p); void display_prism_info(int indentby, prism *prsm); diff --git a/utils/geom.c b/utils/geom.c index 31471e1..f9fcd29 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -129,6 +129,7 @@ void geom_fix_object(geometric_object o) { // FIXME not really sure what to do here } + break; case GEOM COMPOUND_GEOMETRIC_OBJECT: { int i; @@ -211,19 +212,6 @@ void geom_initialize(void) geometry.items = 0; } -/**************************************************************************/ -/* like intersect_line_with_object, but allows for arbitrarily many */ -/* intersections (not just 2) to accommodate non-star-shaped objects. */ -/* if svalues is nonzero, it points on return to a newly allocated buffer */ -/* that must be free()d by the caller. */ -/**************************************************************************/ -#if HAVE_PRISM -int intersect_line_with_prism(vector3 p, vector3 d, prism *prsm, - double **svalues) -{ -} -#endif - /**************************************************************************/ /* Return whether or not the point p (in the lattice basis) is inside @@ -318,7 +306,7 @@ boolean point_in_fixed_pobjectp(vector3 p, geometric_object *o) } case GEOM PRISM: { - return point_in_prism(r, o->subclass.prism_data); + return point_in_prism(o->subclass.prism_data, p); } case GEOM COMPOUND_GEOMETRIC_OBJECT: { @@ -479,7 +467,7 @@ vector3 normal_to_fixed_object(vector3 p, geometric_object o) } // case GEOM BLOCK case GEOM PRISM: { - return normal_to_prism(r, o.subclass.prism_data); + return normal_to_prism(o.subclass.prism_data, r); } default: return r; @@ -907,7 +895,7 @@ int intersect_line_with_object(vector3 p, vector3 d, geometric_object o, } // case GEOM BLOCK case GEOM PRISM: { double *sprism=0; - int n=intersect_line_with_prism(p, d, o.subclass.prism_data, &sprism); + int n=intersect_line_with_prism(o.subclass.prism_data, p, d, &sprism); CHECK( n<=2, "more than 2 intersections in intersect_line_with_object"); memcpy(s, sprism, n*sizeof(double)); if (n>0) free(sprism); @@ -1186,10 +1174,10 @@ void geom_get_bounding_box(geometric_object o, geom_box *box) break; } case GEOM PRISM: + { get_prism_bounding_box(o.subclass.prism_data, box); - box->low=vector3_plus(box->low,o.center); - box->high=vector3_plus(box->high,o.center); break; + } case GEOM COMPOUND_GEOMETRIC_OBJECT: { int i; diff --git a/utils/geom.scm b/utils/geom.scm index 10d22ae..42aa40e 100644 --- a/utils/geom.scm +++ b/utils/geom.scm @@ -91,8 +91,10 @@ (define-class prism geometric-object (define-property vertices '() (make-list-type 'vector3)) + (define-property centroid (vector3 0 0 0) 'vector3) (define-property height 0 'number) - (define-property m_c2p identity_matrix 'matrix3x3)) + (define-property m_c2p identity_matrix 'matrix3x3) + (define-property m_p2c identity_matrix 'matrix3x3)) (define-class ellipsoid block (define-derived-property inverse-semi-axes 'vector3 diff --git a/utils/prism.c b/utils/prism.c index 9a8fdb9..4855a43 100644 --- a/utils/prism.c +++ b/utils/prism.c @@ -69,6 +69,22 @@ using namespace ctlio; #define K_PI 3.14159265358979323846 #define CHECK(cond, s) if (!(cond)){fprintf(stderr,s "\n");exit(EXIT_FAILURE);} +/***************************************************************/ +/* 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_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_vector_c2p(prism *prsm, vector3 vc) +{ return matrix3x3_vector3_mult(prsm->m_c2p, vc); } + /***************************************************************/ /***************************************************************/ /***************************************************************/ @@ -77,13 +93,11 @@ 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 = matrix3x3_vector3_mult(prsm->m_c2p,z0); // set x,y coordinates of low, high to bounding box // of prism base polygon (in prism coordinate system) - vector3 low = vertices[0]; - vector3 high = vertices[0]; + vector3 low = vertices[0]; + vector3 high = vertices[0]; for(int nv=1; nvm_c2p); - box->low = matrix3x3_vector3_mult(m_p2c, low); - box->high = matrix3x3_vector3_mult(m_p2c, high); + box->low = prism_coordinate_p2c(prsm, low); + box->high = prism_coordinate_p2c(prsm, high); } /***************************************************************/ @@ -112,8 +125,8 @@ void get_prism_bounding_box(prism *prsm, geom_box *box) boolean intersect_line_with_segment(double px, double py, double dx, double dy, vector3 v1, vector3 v2, double *s) { - double ax = v1.x, ay = v1.y; - double bx = v2.x-v1.x, by = v2.y-v1.y; + 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; @@ -123,6 +136,22 @@ boolean intersect_line_with_segment(double px, double py, double dx, double dy, return 0; double t = (M00*RHSy-M10*RHSx)/DetM; +/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/ +double ss = (M11*RHSx-M01*RHSy)/DetM; +char *str=getenv("PRISM_RAYS"); +if (str&&str[0]=='1') +{ +FILE *f=fopen("/tmp/prism_rays","a"); +fprintf(f,"#(px,py)=(%e,%e) (dx,dy)=(%e,%e) t=%e s=%e (%s)\n",px,py,dx,dy,t,ss,(t<0.0||t>1.0) ? "no" : "yes"); +fprintf(f,"%e %e %e %e \n",px,py,dx,dy); +fprintf(f,"%e %e %e %e \n",px+ss*dx,py+ss*dy,dx,dy); +fprintf(f,"\n"); +fprintf(f,"%e %e %e %e \n",px,py,dx,dy); +fprintf(f,"%e %e %e %e \n",ax+t*bx,ay+t*by,dx,dy); +fprintf(f,"\n"); +fclose(f); +} +/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/ if (t<0.0 || t>1.0) // intersection of lines does not lie between vertices return 0; @@ -136,7 +165,7 @@ boolean intersect_ray_with_segment(double px, double py, double dx, double dy, { double ss; int status=intersect_line_with_segment(px,py,dx,dy,v1,v2,&ss); - if (status==0 && ss<0.0) + if (status==0 || ss<0.0) return 0; if (s) *s=ss; return 1; @@ -155,8 +184,8 @@ boolean point_in_polygon(double px, double py, vector3 *vertices, int num_vertic int num_side_intersections=0; for(int nv=0; nvvertices.num_items; vector3 *vertices = prsm->vertices.items; @@ -175,14 +204,16 @@ int intersect_line_with_prism(vector3 p, vector3 d, prism *prsm, double **sarray // get coords of p (line origin) and components of d (line slope) // in prism coordinate system - matrix3x3 m_c2p = prsm->m_c2p; - vector3 p_prsm = matrix3x3_vector3_mult(m_c2p, p); - vector3 d_prsm = matrix3x3_vector3_mult(m_c2p, d); + vector3 p_prsm = prism_coordinate_c2p(prsm,p); + vector3 d_prsm = prism_vector_c2p(prsm,d); - int num_intersections=0; - double *svalues=0; + // use length of first edge as a length scale for judging + // lengths to be small or large + double length_scale = vector3_norm(vector3_minus(vertices[1], vertices[0])); // identify intersections of line with all prism side surfaces + int num_intersections=0; + double *svalues=0; for(int nv=0; nv0.0 && (fabs(z_intersect) > 0.5*height ) ) + if ( (z_intersect < -1.0e-7*length_scale) || (z_intersect>(1 + 1.0e-7)*length_scale) ) continue; num_intersections++; @@ -229,7 +260,7 @@ int intersect_line_with_prism(vector3 p, vector3 d, prism *prsm, double **sarray } -vector3 normal_to_prism(vector3 p, prism *prsm) +vector3 normal_to_prism(prism *prsm, vector3 p) { // not implemented yet } @@ -257,15 +288,15 @@ void display_prism_info(int indentby, prism *prsm) /***************************************************************/ /***************************************************************/ /***************************************************************/ -boolean point_in_prism(vector3 p, prism *prsm) +boolean point_in_prism(prism *prsm, vector3 xc) { vector3 *vertices = prsm->vertices.items; int num_vertices = prsm->vertices.num_items; double height = prsm->height; - vector3 p_prsm = matrix3x3_vector3_mult(prsm->m_c2p, p); - if ( fabs(p_prsm.z) > 0.5*prsm->height) + vector3 xp = prism_coordinate_c2p(prsm, xc); + if ( fabs(xp.z) > 0.5*prsm->height) return 0; - return point_in_polygon(p_prsm.x, p_prsm.y, vertices, num_vertices); + return point_in_polygon(xp.x, xp.y, vertices, num_vertices); } /***************************************************************/ @@ -303,7 +334,7 @@ GEOMETRIC_OBJECT make_prism(MATERIAL_TYPE material, centroid = vector3_plus(centroid, vertices[nv]); centroid = vector3_scale(1.0/((double)num_vertices), centroid); - // make sure all vertices lie in a plane and that plane is normal to axis + // make sure all vertices lie in a plane normal to the given axis vector3 zhat = unit_vector3(axis); for(int nv=0; nvheight = height; - prsm->m_c2p = m_c2p; + prsm->centroid = centroid; + prsm->height = height; + prsm->m_p2c = m_p2c; + prsm->m_c2p = m_c2p; // 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 identical and in principle need not + // 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)); for(int nv=0; nvvertices.items[nv] - =matrix3x3_vector3_mult(m_c2p, vector3_minus(vertices[nv],centroid)); - + prsm->vertices.items[nv] = prism_coordinate_c2p(prsm,vertices[nv]); + + // note the center and centroid are different! vector3 center = vector3_plus(centroid, vector3_scale(0.5*height,zhat) ); geometric_object o=make_geometric_object(material, center); o.which_subclass=GEOM PRISM; diff --git a/utils/test-prism.c b/utils/test-prism.c index 57de399..1597f5d 100644 --- a/utils/test-prism.c +++ b/utils/test-prism.c @@ -5,11 +5,15 @@ #include #include #include +#include #include "ctlgeom.h" #define K_PI 3.141592653589793238462643383279502884197 +#define NUMPTS 10000 +#define NUMLINES 1000 + static vector3 make_vector3(double x, double y, double z) { vector3 v; @@ -17,6 +21,8 @@ static vector3 make_vector3(double x, double y, double z) return v; } +boolean point_in_polygon(double px, double py, vector3 *vertices, int num_vertices); + /************************************************************************/ /* return a uniform random number in [a,b] */ /************************************************************************/ @@ -56,45 +62,73 @@ vector3 random_unit_vector3(void) } /***************************************************************/ +/* after running this routine to produce a file named MyFile, */ +/* the prism may be plotted in gnuplot like this: */ +/* gnuplot> splot 'MyFile' u 1:2:3 w lp pt 7 ps 1 */ /***************************************************************/ -/***************************************************************/ -void draw_prism(geometric_object o, char *filename) +void prism2gnuplot(prism *prsm, char *filename) { - vector3 center = o.center; - prism *prsm = o.subclass.prism_data; vector3 *vertices = prsm->vertices.items; int num_vertices = prsm->vertices.num_items; double height = prsm->height; - matrix3x3 m_c2p = prsm->m_c2p; - matrix3x3 m_p2c = matrix3x3_inverse(m_c2p); FILE *f=fopen(filename,"w"); for(int nv=0; nvvertices.items; + int num_vertices = prsm->vertices.num_items; + double height = prsm->height; + vector3 zhat = prsm->m_p2c.c2; + vector3 axis = vector3_scale(height, zhat); + + FILE *f=fopen(filename,"w"); + for(int nv=0; nv 1.0e-7 ) + match=0; + + if (match==0) + num_failed++; + } + fclose(f); + + printf("%i/%i lines failed\n",num_failed,NUMLINES); return num_failed; } +void print_usage(char *msg, int print_usage) +{ + if (!msg) + fprintf(stderr,"%s\n",msg); + if (print_usage) + { printf("usage: \n"); + printf(" --vertexfile MyVertices\n"); + printf(" --height height\n"); + printf(" --axis x y z\n"); + printf("\n"); + printf(" --point x y z\n"); + printf(" --pointfile MyPointFile\n"); + }; + exit(1); +} + +void quit(char *msg) +{ print_usage(msg, 0); } + +void usage(char *msg) +{ print_usage(msg, 1); } + +/************************************************************************/ +/************************************************************************/ /************************************************************************/ -int main(void) +int main(int argc, char *argv[]) { srand(time(NULL)); geom_initialize(); - int num_failed = test1(); - return num_failed; + + if (argc<=1) // if no arguments, run unit tests + { + int num_failed_1 = unit_test1(); + int num_failed_2 = unit_test2(); + return num_failed_1 + num_failed_2; + } + + /***************************************************************/ + /* process arguments *******************************************/ + /***************************************************************/ + vector3 *vertices=0; + int num_vertices=0; + char *vertexfile=0, *pointfile=0; + vector3 zhat={0,0,1}; + double height=1.0; + vector3 test_point={0,0,0}; int have_test_point=0; + for(int narg=1; narg=argc) usage("too few arguments to --axis"); + sscanf(argv[narg+1],"%le",&(zhat.x)); + sscanf(argv[narg+2],"%le",&(zhat.y)); + sscanf(argv[narg+3],"%le",&(zhat.z)); + narg+=3; + } + else if (!strcmp(argv[narg],"--point")) + { if (narg+3>=argc) usage("too few arguments to --point"); + sscanf(argv[narg+1],"%le",&(test_point.x)); + sscanf(argv[narg+2],"%le",&(test_point.y)); + sscanf(argv[narg+3],"%le",&(test_point.z)); + have_test_point=1; + narg+=3; + } + else if (!strcmp(argv[narg],"--height")) + sscanf(argv[++narg],"%le",&height); + else if (!strcmp(argv[narg],"--pointfile")) + pointfile=argv[++narg]; + else + usage("unknown argument"); + } + + /***************************************************************/ + /* read vertices from file and create prism ********************/ + /***************************************************************/ + FILE *f=fopen(vertexfile,"r"); + if (!f) usage("could not open vertexfile"); + char Line[100]; + int LineNum=0; + while( fgets(Line,100,f) ) + { if (Line[0]=='\n' || Line[0]=='#') continue; + num_vertices++; + vector3 v; + if ( 3!=sscanf(Line,"%le %le %le\n",&(v.x),&(v.y),&(v.z)) ) + { fprintf(stderr,"bad vertex on line %i of %s",num_vertices,vertexfile); + exit(1); + } + vertices = (vector3 *)realloc(vertices, num_vertices*sizeof(vector3)); + vertices[num_vertices-1]=v; + }; + fclose(f); + + geometric_object the_prism=make_prism(NULL, vertices, num_vertices, height, zhat); + prism *prsm=the_prism.subclass.prism_data; + prism2gmsh(prsm, "test-prism.geo"); + prism2gnuplot(prsm, "test-prism.gp"); + + + /***************************************************************/ + /* read points from file or generate random points and check */ + /* inclusion */ + /***************************************************************/ + geom_box box; + get_prism_bounding_box(prsm, &box); + f=fopen("prism-box.gp","w"); + fprintf(f,"%e %e %e \n",box.low.x,box.low.y,box.low.z); + fprintf(f,"%e %e %e \n",box.high.x,box.low.y,box.low.z); + fprintf(f,"%e %e %e \n",box.high.x,box.high.y,box.low.z); + fprintf(f,"%e %e %e \n",box.low.x,box.high.y,box.low.z); + fprintf(f,"\n\n"); + + fprintf(f,"%e %e %e \n",box.low.x,box.low.y,box.high.z); + fprintf(f,"%e %e %e \n",box.high.x,box.low.y,box.high.z); + fprintf(f,"%e %e %e \n",box.high.x,box.high.y,box.high.z); + fprintf(f,"%e %e %e \n",box.low.x,box.high.y,box.high.z); + fprintf(f,"\n\n"); + + fprintf(f,"%e %e %e \n",box.low.x,box.low.y,box.low.z); + fprintf(f,"%e %e %e \n",box.high.x,box.low.y,box.low.z); + fprintf(f,"%e %e %e \n",box.high.x,box.low.y,box.high.z); + fprintf(f,"%e %e %e \n",box.low.x,box.low.y,box.high.z); + fprintf(f,"\n\n"); + + fprintf(f,"%e %e %e \n",box.low.x,box.high.y,box.low.z); + fprintf(f,"%e %e %e \n",box.high.x,box.high.y,box.low.z); + fprintf(f,"%e %e %e \n",box.high.x,box.high.y,box.high.z); + fprintf(f,"%e %e %e \n",box.low.x,box.high.y,box.high.z); + fprintf(f,"\n\n"); + + fprintf(f,"%e %e %e \n",box.low.x,box.low.y,box.low.z); + fprintf(f,"%e %e %e \n",box.low.x,box.high.y,box.low.z); + fprintf(f,"%e %e %e \n",box.low.x,box.high.y,box.high.z); + fprintf(f,"%e %e %e \n",box.low.x,box.low.y,box.high.z); + fprintf(f,"\n\n"); + + fprintf(f,"%e %e %e \n",box.high.x,box.low.y,box.low.z); + fprintf(f,"%e %e %e \n",box.high.x,box.high.y,box.low.z); + fprintf(f,"%e %e %e \n",box.high.x,box.high.y,box.high.z); + fprintf(f,"%e %e %e \n",box.high.x,box.low.y,box.high.z); + fprintf(f,"\n\n"); + fclose(f); + + vector3 delta = vector3_minus(box.high, box.low); + box.high = vector3_plus(box.high,delta); + box.low = vector3_minus(box.low,delta); + + f = (pointfile ? fopen(pointfile,"r") : 0); + + FILE *gpfile1 = fopen("test-prism.in.gp","w"); + FILE *gpfile2 = fopen("test-prism.out.gp","w"); + FILE *ppfile1 = fopen("test-prism.in.pp","w"); + FILE *ppfile2 = fopen("test-prism.out.pp","w"); + fprintf(ppfile1,"View \"Points inside prism\" {\n"); + fprintf(ppfile2,"View \"Points outside prism\" {\n"); + int Done=0; + int NumPoints=0; + while(!Done) + { + vector3 v; + if (f) + { char Line[100]; + NumPoints++; + if (fgets(Line,100,f)) + sscanf(Line,"%le %le %le\n",&(v.x),&(v.y),&(v.z)); + else + Done=1; + } + else if (have_test_point) + { v=test_point; + Done=1; + } + else + { v.x = myurand(box.low.x, box.high.x); + v.y = myurand(box.low.y, box.high.y); + v.z = myurand(box.low.z, box.high.z); + if (++NumPoints == NUMPTS ) Done=1; + } + + // vector3 vp=prism_coordinate_c2p(prsm,v); + // vp.z=0; + // if (point_in_polygon(vp.x,vp.y,prsm->vertices.items,num_vertices)) + if (point_in_objectp(v,the_prism)) + { fprintf(gpfile1,"%e %e %e\n",v.x,v.y,v.z); + fprintf(ppfile1,"SP(%e,%e,%e) {0};\n",v.x,v.y,v.z); + if (have_test_point) printf("Point in object.\n"); + } + else + { fprintf(gpfile2,"%e %e %e\n",v.x,v.y,v.z); + fprintf(ppfile2,"SP(%e,%e,%e) {0};\n",v.x,v.y,v.z); + if (have_test_point) printf("Point not in object.\n"); + } + } + fprintf(ppfile1,"};\n"); + fprintf(ppfile2,"};\n"); + fclose(gpfile1); + fclose(gpfile2); + fclose(ppfile1); + fclose(ppfile2); } From 21c71ebd1676dc1ff528d34b7da034d43a8eccf9 Mon Sep 17 00:00:00 2001 From: Homer Reid Date: Thu, 22 Mar 2018 09:05:19 +0900 Subject: [PATCH 04/10] implemented intersect_line_with_object for prisms; added working unit test to test-prism.c --- utils/geom.c | 7 ++- utils/prism.c | 143 +++++++++++++++++++++++++++++++++------------ utils/test-prism.c | 89 ++++++++++++++++++++-------- 3 files changed, 176 insertions(+), 63 deletions(-) diff --git a/utils/geom.c b/utils/geom.c index f9fcd29..cce75aa 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -467,7 +467,7 @@ vector3 normal_to_fixed_object(vector3 p, geometric_object o) } // case GEOM BLOCK case GEOM PRISM: { - return normal_to_prism(o.subclass.prism_data, r); + return normal_to_prism(o.subclass.prism_data, p); } default: return r; @@ -736,6 +736,7 @@ void CTLIO display_geometric_object_info(int indentby, geometric_object o) int intersect_line_with_object(vector3 p, vector3 d, geometric_object o, double s[2]) { + vector3 p0=p; p = vector3_minus(p, o.center); s[0] = s[1] = 0; switch (o.which_subclass) { @@ -895,7 +896,9 @@ int intersect_line_with_object(vector3 p, vector3 d, geometric_object o, } // case GEOM BLOCK case GEOM PRISM: { double *sprism=0; - int n=intersect_line_with_prism(o.subclass.prism_data, p, d, &sprism); + int n=intersect_line_with_prism(o.subclass.prism_data, p0, d, &sprism); +if (n>2) + printf("--point %e %e %e --dir %e %e %e\n",p0.x,p0.y,p0.z,d.x,d.y,d.z); CHECK( n<=2, "more than 2 intersections in intersect_line_with_object"); memcpy(s, sprism, n*sizeof(double)); if (n>0) free(sprism); diff --git a/utils/prism.c b/utils/prism.c index 4855a43..6d06324 100644 --- a/utils/prism.c +++ b/utils/prism.c @@ -136,26 +136,11 @@ boolean intersect_line_with_segment(double px, double py, double dx, double dy, return 0; double t = (M00*RHSy-M10*RHSx)/DetM; -/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/ -double ss = (M11*RHSx-M01*RHSy)/DetM; -char *str=getenv("PRISM_RAYS"); -if (str&&str[0]=='1') -{ -FILE *f=fopen("/tmp/prism_rays","a"); -fprintf(f,"#(px,py)=(%e,%e) (dx,dy)=(%e,%e) t=%e s=%e (%s)\n",px,py,dx,dy,t,ss,(t<0.0||t>1.0) ? "no" : "yes"); -fprintf(f,"%e %e %e %e \n",px,py,dx,dy); -fprintf(f,"%e %e %e %e \n",px+ss*dx,py+ss*dy,dx,dy); -fprintf(f,"\n"); -fprintf(f,"%e %e %e %e \n",px,py,dx,dy); -fprintf(f,"%e %e %e %e \n",ax+t*bx,ay+t*by,dx,dy); -fprintf(f,"\n"); -fclose(f); -} -/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/ + if (s) *s = (M11*RHSx-M01*RHSy)/DetM; + if (t<0.0 || t>1.0) // intersection of lines does not lie between vertices return 0; - if (s) *s = (M11*RHSx-M01*RHSy)/DetM; return 1; } @@ -231,7 +216,10 @@ int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double **sarray // 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; - if ( (z_intersect < -1.0e-7*length_scale) || (z_intersect>(1 + 1.0e-7)*length_scale) ) + double AbsTol = 1.0e-7*length_scale; + double z_min = 0.0 - AbsTol; + double z_max = height + AbsTol; + if ( (z_intersectz_max) ) continue; num_intersections++; @@ -242,10 +230,11 @@ int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double **sarray }; // identify intersections of line with prism top and bottom surfaces - double dz=d_prsm.z; - if (dz > 1.0e-10*vector3_norm(d_prsm)) - for(double sign=-1.0; sign<=1.01; sign+=2.0) - { double s = (p_prsm.z - sign*0.5*height)/dz; + double dz = d_prsm.z; + if ( fabs(dz) > 1.0e-7*vector3_norm(d_prsm)) + for(int 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)) continue; num_intersections++; @@ -260,10 +249,102 @@ int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double **sarray } -vector3 normal_to_prism(prism *prsm, vector3 p) -{ // not implemented yet +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ +boolean point_in_prism(prism *prsm, vector3 xc) +{ + 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) + return 0; + return point_in_polygon(xp.x, xp.y, vertices, num_vertices); } +/***************************************************************/ +/* compute the minimum-length vector from p to the plane */ +/* containing o (origin) and spanned by basis vectors v1,v2 */ +/* algorithm: solve the 3x3 system p-s*v3 = o + t*v1 + u*v2 */ +/* where v3 = v1 x v2 and s,t,u are unknowns */ +/***************************************************************/ +vector3 normal_to_plane(vector3 o, vector3 v1, vector3 v2, vector3 p) +{ + vector3 RHS = vector3_minus(p,o); + + // handle the degenerate-to-2D case + if ( (fabs(v1.z) + fabs(v2.z)) < 2.0e-7 && fabs(RHS.z)<1.0e-7 ) + { vector3 zhat={0,0,1}; + vector3 v3 = vector3_cross(zhat, v1); + double M00 = v1.x; double M10 = v3.x; + double M01 = v1.y; double M11 = v3.y; + double DetM = M00*M11 - M01*M10; + if ( fabs(DetM) < 1.0e-10 ) + return vector3_scale(0.0, v3); + // double t = (M00*RHSy-M10*RHSx)/DetM; + double s= (M11*RHS.x-M01*RHS.y)/DetM; + return vector3_scale(-1.0*s,v3); + } + + vector3 v3 = vector3_cross(v1, v2); + matrix3x3 M; + M.c0 = v1; + M.c1 = v2; + M.c2 = vector3_scale(-1.0, v3); + vector3 tus = matrix3x3_vector3_mult(matrix3x3_inverse(M),RHS); + return vector3_scale(-1.0*tus.z, v3); +} + +/***************************************************************/ +/* find the face of the prism for which the normal distance */ +/* from x 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 centroid = prsm->centroid; + double height = prsm->height; + vector3 *vertices = prsm->vertices.items; + int num_vertices = prsm->vertices.num_items; + + vector3 xp = prism_coordinate_c2p(prsm, xc); + vector3 axis = {0,0,0}; axis.z=height; + + vector3 retval; + double min_distance; + // side walls + for(int nv=0; nvvertices.items; - int num_vertices = prsm->vertices.num_items; - double height = prsm->height; - vector3 xp = prism_coordinate_c2p(prsm, xc); - if ( fabs(xp.z) > 0.5*prsm->height) - return 0; - return point_in_polygon(xp.x, xp.y, vertices, num_vertices); -} - /***************************************************************/ // like vector3_equal but tolerant of small floating-point discrepancies /***************************************************************/ diff --git a/utils/test-prism.c b/utils/test-prism.c index 1597f5d..2ec7093 100644 --- a/utils/test-prism.c +++ b/utils/test-prism.c @@ -14,6 +14,10 @@ #define NUMPTS 10000 #define NUMLINES 1000 +#define LX 0.5 +#define LY 1.0 +#define LZ 1.5 + static vector3 make_vector3(double x, double y, double z) { vector3 v; @@ -75,10 +79,10 @@ void prism2gnuplot(prism *prsm, char *filename) FILE *f=fopen(filename,"w"); for(int nv=0; nvs_block[1]) + { double s=s_block[0]; s_block[0]=s_block[1]; s_block[1]=s; } + if (ns_prism==2 && s_prism[0]>s_prism[1]) + { double s=s_prism[0]; s_prism[0]=s_prism[1]; s_prism[1]=s; } + int match = (ns_block==ns_prism) ? 1 : 0; for(int n=0; match==1 && n 1.0e-7 ) @@ -268,8 +279,9 @@ int main(int argc, char *argv[]) int num_vertices=0; char *vertexfile=0, *pointfile=0; vector3 zhat={0,0,1}; - double height=1.0; + double height=1.5; vector3 test_point={0,0,0}; int have_test_point=0; + vector3 test_dir={0,0,0}; int have_test_dir=0; for(int narg=1; narg=argc) usage("too few arguments to --dir"); + sscanf(argv[narg+1],"%le",&(test_dir.x)); + sscanf(argv[narg+2],"%le",&(test_dir.y)); + sscanf(argv[narg+3],"%le",&(test_dir.z)); + have_test_dir=1; + narg+=3; + } else if (!strcmp(argv[narg],"--height")) sscanf(argv[++narg],"%le",&height); else if (!strcmp(argv[narg],"--pointfile")) @@ -402,9 +422,6 @@ int main(int argc, char *argv[]) if (++NumPoints == NUMPTS ) Done=1; } - // vector3 vp=prism_coordinate_c2p(prsm,v); - // vp.z=0; - // if (point_in_polygon(vp.x,vp.y,prsm->vertices.items,num_vertices)) if (point_in_objectp(v,the_prism)) { fprintf(gpfile1,"%e %e %e\n",v.x,v.y,v.z); fprintf(ppfile1,"SP(%e,%e,%e) {0};\n",v.x,v.y,v.z); @@ -415,6 +432,32 @@ int main(int argc, char *argv[]) fprintf(ppfile2,"SP(%e,%e,%e) {0};\n",v.x,v.y,v.z); if (have_test_point) printf("Point not in object.\n"); } + + if (have_test_dir) + { +#if 0 + double s_block[2]; + int ns_block=intersect_line_with_object(v, test_dir, the_block, s_block); + FILE *f=fopen("/tmp/test-prism.blocklines","w"); + for(int ns=0; ns Date: Wed, 28 Mar 2018 15:06:06 +0900 Subject: [PATCH 05/10] updates --- utils/ctlgeom.h | 3 +++ utils/geom.c | 46 ++++++++++++++++++++++++++-------------------- utils/prism.c | 44 +++++++++++++++++++++++++------------------- 3 files changed, 54 insertions(+), 39 deletions(-) diff --git a/utils/ctlgeom.h b/utils/ctlgeom.h index 8acc197..7052792 100644 --- a/utils/ctlgeom.h +++ b/utils/ctlgeom.h @@ -69,6 +69,9 @@ extern vector3 normal_to_object(vector3 p, GEOMETRIC_OBJECT o); extern vector3 normal_to_fixed_object(vector3 p, GEOMETRIC_OBJECT o); extern int intersect_line_with_object(vector3 p, vector3 d, GEOMETRIC_OBJECT o, double s[2]); +extern double intersect_line_segment_with_object(vector3 p, vector3 d, + GEOMETRIC_OBJECT o, + double a, double b); extern MATERIAL_TYPE material_of_point_inobject(vector3 p, boolean *inobject); extern MATERIAL_TYPE material_of_point_inobject0( GEOMETRIC_OBJECT_LIST geometry, vector3 p, boolean *inobject); diff --git a/utils/geom.c b/utils/geom.c index cce75aa..fce5362 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -62,6 +62,9 @@ using namespace ctlio; #define K_PI 3.14159265358979323846 #define CHECK(cond, s) if (!(cond)){fprintf(stderr,s "\n");exit(EXIT_FAILURE);} +#define MAX(a,b) ((a) > (b) ? (a) : (b)) +#define MIN(a,b) ((a) < (b) ? (a) : (b)) + /**************************************************************************/ /* If v is a vector in the lattice basis, normalize v so that @@ -894,21 +897,33 @@ int intersect_line_with_object(vector3 p, vector3 d, geometric_object o, } } } // case GEOM BLOCK - case GEOM PRISM: - { double *sprism=0; - int n=intersect_line_with_prism(o.subclass.prism_data, p0, d, &sprism); -if (n>2) - printf("--point %e %e %e --dir %e %e %e\n",p0.x,p0.y,p0.z,d.x,d.y,d.z); - CHECK( n<=2, "more than 2 intersections in intersect_line_with_object"); - memcpy(s, sprism, n*sizeof(double)); - if (n>0) free(sprism); - return n; - } // case GEOM PRISM + default: return 0; } } +/* Compute the intersections with o of a line along p+s*d in the interval s in [a,b], returning + the length of the s intersection in this interval. (Note: o must not be a compound object.) */ +double intersect_line_segment_with_object(vector3 p, vector3 d, geometric_object o, double a, double b) +{ + if (o.which_subclass==GEOM PRISM) + { + return intersect_line_segment_with_prism(o.subclass.prism_data, p, d, a, b); + } + else + { double s[2]; + if (2 == intersect_line_with_object(p, d, o, s)) + { double ds = (s[0] < s[1] ? MIN(s[1],b) - MAX(s[0],a) + : MIN(s[0],b) - MAX(s[1],a) + ); + return (ds > 0 ? ds : 0.0); + } + else + return 0.0; + } +} + /**************************************************************************/ /* Given a basis (matrix columns are the basis unit vectors) and the @@ -961,9 +976,6 @@ matrix3x3 CTLIO square_basis(matrix3x3 basis, vector3 size) /**************************************************************************/ /* geom_box utilities: */ - -#define MAX(a,b) ((a) > (b) ? (a) : (b)) -#define MIN(a,b) ((a) < (b) ? (a) : (b)) static void geom_box_union(geom_box *bu, const geom_box *b1, const geom_box *b2) { @@ -1264,13 +1276,7 @@ static double overlap_integrand(integer ndim, number *x, void *data_) a0 = data->c0 - w; b0 = data->c0 + w; } - if (2 == intersect_line_with_object(p, data->dir, data->o, s)) { - double ds = (s[0] < s[1] - ? MIN(s[1],b0) - MAX(s[0],a0) - : MIN(s[0],b0) - MAX(s[1],a0)); - return (ds > 0 ? ds * scale_result : 0.0); - } - return 0.0; + return intersect_line_segment_with_object(p, data->dir, data->o, a0, b0) * scale_result; } number overlap_with_object(geom_box b, int is_ellipsoid, geometric_object o, diff --git a/utils/prism.c b/utils/prism.c index 6d06324..89df513 100644 --- a/utils/prism.c +++ b/utils/prism.c @@ -176,12 +176,10 @@ boolean point_in_polygon(double px, double py, vector3 *vertices, int num_vertic } /***************************************************************/ -/* return the number of intersections of the line p + s*d with */ -/* surfaces of the prism. */ -/* if **sarray is non-null, on return it points to a newly */ -/* allocated array of length equal to the return value. */ +/* return the number of intersections of prism surfaces with */ +/* the line segment p + s*d , avertices.num_items; vector3 *vertices = prsm->vertices.items; @@ -196,9 +194,12 @@ int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double **sarray // lengths to be small or large double length_scale = vector3_norm(vector3_minus(vertices[1], vertices[0])); - // identify intersections of line with all prism side surfaces + // identify s-values s0, s1, ... of line-segment intersections + // with all prism side surfaces. + // measure[0] = (s_2-s_1) + (s_4-s_3) + ... + // measure[1] = (s_1-s_0) + (s_3-s_2) + ... int num_intersections=0; - double *svalues=0; + double last_s, measure[2]={0.0,0.0}; for(int nv=0; nvz_max) ) continue; + // Now we know the line p+s*d intersects the prism...does it do so + // within the range of the specified line segment? + if (sb) + continue; + + if (num_intersections>0) + measure[num_intersections%2] += s-last_s; num_intersections++; - if (sarray) - { svalues = (double *)realloc(svalues, num_intersections*sizeof(double)); - svalues[num_intersections-1] = s; - }; + last_s=s; }; // identify intersections of line with prism top and bottom surfaces @@ -235,18 +240,19 @@ int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double **sarray for(int LowerUpper=0; LowerUpper<2; LowerUpper++) { double z0 = LowerUpper ? height : 0.0; double s = (z0 - p_prsm.z)/dz; + if (sb) + continue; if (!point_in_polygon(p_prsm.x + s*d_prsm.x, p_prsm.y+s*d_prsm.y, vertices, num_vertices)) continue; - num_intersections++; - if (sarray) - { svalues = (double *)realloc(svalues, num_intersections*sizeof(double)); - svalues[num_intersections-1] = s; - }; - }; - if (sarray) *sarray=svalues; - return num_intersections; + if (num_intersections>0) + measure[num_intersections%2] += s-last_s; + num_intersections++; + last_s=s; + } + //if (num_intersections%2) // point lies inside object + return 0.0; } /***************************************************************/ From e06c6d32dd8ff11c06908704f79909df5e14d190 Mon Sep 17 00:00:00 2001 From: Homer Reid Date: Wed, 28 Mar 2018 16:27:09 -0400 Subject: [PATCH 06/10] new function intersect_line_segment_with_object --- utils/prism.c | 47 +++++++++++++++++------- utils/test-prism.c | 91 ++++++++++++++++++++-------------------------- 2 files changed, 74 insertions(+), 64 deletions(-) diff --git a/utils/prism.c b/utils/prism.c index 89df513..2b3a94e 100644 --- a/utils/prism.c +++ b/utils/prism.c @@ -175,12 +175,36 @@ boolean point_in_polygon(double px, double py, vector3 *vertices, int num_vertic return num_side_intersections%2; } +void add_to_list(double s, double *slist, int length) +{ + switch(length) + { case 0: + slist[0] = s; + break; + case 1: + if (s>=slist[0]) + slist[1]=s; + else + { slist[1]=slist[0]; slist[0]=s; } + break; + default: + if (svertices.num_items; vector3 *vertices = prsm->vertices.items; double height = prsm->height; @@ -196,10 +220,8 @@ double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, doub // identify s-values s0, s1, ... of line-segment intersections // with all prism side surfaces. - // measure[0] = (s_2-s_1) + (s_4-s_3) + ... - // measure[1] = (s_1-s_0) + (s_3-s_2) + ... + double slist[2]; int num_intersections=0; - double last_s, measure[2]={0.0,0.0}; for(int nv=0; nvb) continue; - if (num_intersections>0) - measure[num_intersections%2] += s-last_s; - num_intersections++; - last_s=s; + add_to_list(s, slist, num_intersections++); }; // identify intersections of line with prism top and bottom surfaces @@ -245,14 +264,16 @@ double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, doub if (!point_in_polygon(p_prsm.x + s*d_prsm.x, p_prsm.y+s*d_prsm.y, vertices, num_vertices)) continue; - if (num_intersections>0) - measure[num_intersections%2] += s-last_s; - num_intersections++; - last_s=s; + add_to_list(s, slist, num_intersections++); } - //if (num_intersections%2) // point lies inside object - return 0.0; + if (num_intersections==0) + return 0.0; + if (num_intersections%2) // point lies inside object + return slist[0]-a; + else + return slist[1]-slist[0]; + } /***************************************************************/ diff --git a/utils/test-prism.c b/utils/test-prism.c index 2ec7093..a09ad14 100644 --- a/utils/test-prism.c +++ b/utils/test-prism.c @@ -171,8 +171,8 @@ int unit_test1() /************************************************************************/ /* second unit test: create the same parallelepiped in two ways, as a */ /* block and as a prism, then check intersections of 1000 */ -/* randomly generated lines with the objects and confirm the */ -/* results agree. */ +/* randomly generated line segments with the objects and */ +/* confirm that the results agree. */ /************************************************************************/ int unit_test2() { @@ -208,24 +208,14 @@ int unit_test2() d.y = SinTheta*sin(Phi); d.z = CosTheta; - double s_block[2]; - int ns_block=intersect_line_with_object(p, d, the_block, s_block); + // randomly generated endpoints of line segment + double a = myurand(-1.0, 1.0); + double b = a + myurand(-1.0, 1.0); - double s_prism[2]; - int ns_prism=intersect_line_with_object(p, d, the_prism, s_prism); - - // if there are 2 s-values they may be reported in different orders - // so sort them before comparing - if (ns_block==2 && s_block[0]>s_block[1]) - { double s=s_block[0]; s_block[0]=s_block[1]; s_block[1]=s; } - if (ns_prism==2 && s_prism[0]>s_prism[1]) - { double s=s_prism[0]; s_prism[0]=s_prism[1]; s_prism[1]=s; } - - int match = (ns_block==ns_prism) ? 1 : 0; - for(int n=0; match==1 && n 1.0e-7 ) - match=0; + double l_block=intersect_line_segment_with_object(p, d, the_block, a, b); + double l_prism=intersect_line_segment_with_object(p, d, the_prism, a, b); + int match = fabs(l_block - l_prism) > 1.0e-7; if (match==0) num_failed++; } @@ -247,6 +237,7 @@ void print_usage(char *msg, int print_usage) printf("\n"); printf(" --point x y z\n"); printf(" --pointfile MyPointFile\n"); + printf(" --lineseg dx dy dz a b\n"); }; exit(1); } @@ -278,18 +269,19 @@ int main(int argc, char *argv[]) vector3 *vertices=0; int num_vertices=0; char *vertexfile=0, *pointfile=0; - vector3 zhat={0,0,1}; + vector3 axis={0,0,1}; double height=1.5; vector3 test_point={0,0,0}; int have_test_point=0; - vector3 test_dir={0,0,0}; int have_test_dir=0; + vector3 test_dir={0,0,0}; int have_test_lineseg=0; + double test_a, test_b; for(int narg=1; narg=argc) usage("too few arguments to --axis"); - sscanf(argv[narg+1],"%le",&(zhat.x)); - sscanf(argv[narg+2],"%le",&(zhat.y)); - sscanf(argv[narg+3],"%le",&(zhat.z)); + sscanf(argv[narg+1],"%le",&(axis.x)); + sscanf(argv[narg+2],"%le",&(axis.y)); + sscanf(argv[narg+3],"%le",&(axis.z)); narg+=3; } else if (!strcmp(argv[narg],"--point")) @@ -300,13 +292,15 @@ int main(int argc, char *argv[]) have_test_point=1; narg+=3; } - else if (!strcmp(argv[narg],"--dir")) - { if (narg+3>=argc) usage("too few arguments to --dir"); + else if (!strcmp(argv[narg],"--lineseg")) + { if (narg+5>=argc) usage("too few arguments to --lineseg"); sscanf(argv[narg+1],"%le",&(test_dir.x)); sscanf(argv[narg+2],"%le",&(test_dir.y)); sscanf(argv[narg+3],"%le",&(test_dir.z)); - have_test_dir=1; - narg+=3; + sscanf(argv[narg+4],"%le",&(test_a)); + sscanf(argv[narg+5],"%le",&(test_b)); + have_test_lineseg=1; + narg+=5; } else if (!strcmp(argv[narg],"--height")) sscanf(argv[++narg],"%le",&height); @@ -336,12 +330,11 @@ int main(int argc, char *argv[]) }; fclose(f); - geometric_object the_prism=make_prism(NULL, vertices, num_vertices, height, zhat); + geometric_object the_prism=make_prism(NULL, vertices, num_vertices, height, axis); prism *prsm=the_prism.subclass.prism_data; prism2gmsh(prsm, "test-prism.geo"); prism2gnuplot(prsm, "test-prism.gp"); - /***************************************************************/ /* read points from file or generate random points and check */ /* inclusion */ @@ -390,6 +383,14 @@ int main(int argc, char *argv[]) box.high = vector3_plus(box.high,delta); box.low = vector3_minus(box.low,delta); + void* m = NULL; + vector3 c = { 0, 0, 0 }; + vector3 xhat = make_vector3(1,0,0); + vector3 yhat = make_vector3(0,1,0); + vector3 zhat = make_vector3(0,0,1); + vector3 size = make_vector3(LX,LY,LZ); + geometric_object the_block = make_block(m, c, xhat, yhat, zhat, size); + f = (pointfile ? fopen(pointfile,"r") : 0); FILE *gpfile1 = fopen("test-prism.in.gp","w"); @@ -433,29 +434,17 @@ int main(int argc, char *argv[]) if (have_test_point) printf("Point not in object.\n"); } - if (have_test_dir) + if (have_test_lineseg) { -#if 0 - double s_block[2]; - int ns_block=intersect_line_with_object(v, test_dir, the_block, s_block); - FILE *f=fopen("/tmp/test-prism.blocklines","w"); - for(int ns=0; ns Date: Wed, 16 May 2018 16:24:18 -0400 Subject: [PATCH 07/10] finished implementing all geometry primitives for prisms; added passing unit test utils/test-prism.c --- utils/Makefile.am | 12 +- utils/ctlgeom.h | 7 +- utils/prism.c | 161 +++++++++-------- utils/test-prism.c | 430 ++++++++++++++++++++------------------------- 4 files changed, 287 insertions(+), 323 deletions(-) diff --git a/utils/Makefile.am b/utils/Makefile.am index 8fa2acf..fc4030d 100644 --- a/utils/Makefile.am +++ b/utils/Makefile.am @@ -13,6 +13,13 @@ libctlgeom_la_CPPFLAGS = -DLIBCTLGEOM geomtst_SOURCES = geomtst.c geomtst_LDADD = libctlgeom.la +check_PROGRAMS = test-prism + +test_prism_SOURCES = test-prism.c +test_prism_LDADD = libctlgeom.la + +TESTS = test-prism + dist_man_MANS = gen-ctl-io.1 BUILT_SOURCES = gen-ctl-io geom-ctl-io.c ctlgeom-types.h nlopt-constants.scm @@ -39,10 +46,5 @@ geom-ctl-io.c: ctl-io.c endif -check_PROGRAMS = test-prism - -test_prism_SOURCES = test-prism.c -test_prism_LDADD = libctlgeom.la - clean-local: rm -f ctl-io.[ch] nlopt-constants.scm diff --git a/utils/ctlgeom.h b/utils/ctlgeom.h index 7052792..f5f76d9 100644 --- a/utils/ctlgeom.h +++ b/utils/ctlgeom.h @@ -150,7 +150,7 @@ GEOMETRIC_OBJECT make_ellipsoid(MATERIAL_TYPE material, vector3 center, // routines in prism.c GEOMETRIC_OBJECT make_prism(MATERIAL_TYPE material, - vector3 *vertices, int num_vertices, + const vector3 *vertices, int num_vertices, double height, vector3 axis); vector3 prism_coordinate_p2c(prism *prsm, vector3 vp); @@ -158,11 +158,12 @@ vector3 prism_coordinate_c2p(prism *prsm, vector3 vc); vector3 prism_vector_p2c(prism *prsm, vector3 vp); vector3 prism_vector_c2p(prism *prsm, vector3 vc); void get_prism_bounding_box(prism *prsm, geom_box *box); -int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double **sarray); +double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b); +boolean point_in_polygon(double px, double py, vector3 *vertices, int num_vertices); boolean point_in_prism(prism *prsm, vector3 p); vector3 normal_to_prism(prism *prsm, vector3 p); void display_prism_info(int indentby, prism *prsm); - +int vector3_nearly_equal(vector3 v1, vector3 v2); /**************************************************************************/ diff --git a/utils/prism.c b/utils/prism.c index 89df513..39971b1 100644 --- a/utils/prism.c +++ b/utils/prism.c @@ -20,17 +20,24 @@ */ /* - * prism.c -- geometry routines for prisms - * a prism is a planar polygon, consisting of any - * number of user-specified vertices, extruded - * through a given thickness (the "height") in the - * direction of a given axis. + * prism.c -- geometry routines 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", + * in which the prism floor lies in the XY plane with centroid + * at the origin and the prism axis is the positive Z-axis. + * homer reid 4/2018 */ #include #include #include +/*--------------------------------------------------------------*/ +/* repeat some definitions from geom.c that are needed here too;*/ +/* can these be moved to ctlgeom.h? */ +/*--------------------------------------------------------------*/ #ifndef LIBCTLGEOM # include "ctl-io.h" #else @@ -40,18 +47,10 @@ #include #ifdef CXX_CTL_IO -using namespace ctlio; -# define CTLIO ctlio:: + using namespace ctlio; # define GEOM geometric_object:: -# define BLK block:: -# define CYL cylinder:: -# define MAT material_type:: #else -# define CTLIO # define GEOM -# define BLK -# define CYL -# define MAT #endif #ifdef __cplusplus @@ -66,7 +65,6 @@ using namespace ctlio; # define FREE1(p) free(p) #endif -#define K_PI 3.14159265358979323846 #define CHECK(cond, s) if (!(cond)){fprintf(stderr,s "\n");exit(EXIT_FAILURE);} /***************************************************************/ @@ -96,9 +94,9 @@ void get_prism_bounding_box(prism *prsm, geom_box *box) // set x,y coordinates of low, high to bounding box // of prism base polygon (in prism coordinate system) - vector3 low = vertices[0]; - vector3 high = vertices[0]; - for(int nv=1; nvvertices.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) + return 0; + return point_in_polygon(xp.x, xp.y, vertices, num_vertices); +} + +/***************************************************************/ +/* insert s into slist, which is sorted in increasing order */ +/* has maximum length 2 */ +/***************************************************************/ +int insert_s_in_list(double s, double slist[2], int length) +{ + if (length==0) + { slist[0]=s; + return 1; + } + if ( s < slist[0] ) + { slist[1]=slist[0]; + slist[0]=s; + } + else if (length==1 || s < slist[1] ) + slist[1]=s; + + return 2; +} + +/***************************************************************/ +/* compute all values of s at which the line p+s*d intersects */ +/* a prism face, retaining the first two entries in a list */ +/* of s values sorted in increasing order m */ +/***************************************************************/ +int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double slist[2]) { int num_vertices = prsm->vertices.num_items; vector3 *vertices = prsm->vertices.items; @@ -194,13 +228,10 @@ double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, doub // lengths to be small or large double length_scale = vector3_norm(vector3_minus(vertices[1], vertices[0])); - // identify s-values s0, s1, ... of line-segment intersections - // with all prism side surfaces. - // measure[0] = (s_2-s_1) + (s_4-s_3) + ... - // measure[1] = (s_1-s_0) + (s_3-s_2) + ... + // identify intersections with prism side faces int num_intersections=0; - double last_s, measure[2]={0.0,0.0}; - for(int nv=0; nvz_max) ) continue; - // Now we know the line p+s*d intersects the prism...does it do so - // within the range of the specified line segment? - if (sb) - continue; - - if (num_intersections>0) - measure[num_intersections%2] += s-last_s; - num_intersections++; - last_s=s; - }; + num_intersections=insert_s_in_list(s, slist, num_intersections); + } - // identify intersections of line with prism top and bottom surfaces + // 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)) - for(int LowerUpper=0; LowerUpper<2; LowerUpper++) + for(LowerUpper=0; LowerUpper<2; LowerUpper++) { double z0 = LowerUpper ? height : 0.0; double s = (z0 - p_prsm.z)/dz; - if (sb) - continue; if (!point_in_polygon(p_prsm.x + s*d_prsm.x, p_prsm.y+s*d_prsm.y, vertices, num_vertices)) continue; - - if (num_intersections>0) - measure[num_intersections%2] += s-last_s; - num_intersections++; - last_s=s; + num_intersections=insert_s_in_list(s,slist,num_intersections); } - //if (num_intersections%2) // point lies inside object - return 0.0; + return num_intersections; } /***************************************************************/ /***************************************************************/ /***************************************************************/ -boolean point_in_prism(prism *prsm, vector3 xc) -{ - 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) - return 0; - return point_in_polygon(xp.x, xp.y, vertices, num_vertices); +double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b) +{ + double slist[2]={0.0, 0.0}; + if (2!=intersect_line_with_prism(prsm, p, d, slist)) + return 0.0; + double ds = fmin(slist[1],b) - fmax(slist[0],a); + return ds > 0.0 ? ds : 0.0; } /***************************************************************/ @@ -319,8 +334,10 @@ vector3 normal_to_prism(prism *prsm, vector3 xc) vector3 retval; double min_distance; + // side walls - for(int nv=0; nvm_c2p); - for(int nv=0; nv=3, "fewer than 3 vertices in make_prism"); // compute centroid of vertices vector3 centroid = {0.0, 0.0, 0.0}; - for(int nv=0; nvvertices.num_items = num_vertices; prsm->vertices.items = (vector3 *)malloc(num_vertices*sizeof(vector3)); - for(int nv=0; nvvertices.items[nv] = prism_coordinate_c2p(prsm,vertices[nv]); // note the center and centroid are different! diff --git a/utils/test-prism.c b/utils/test-prism.c index 2ec7093..76ba38c 100644 --- a/utils/test-prism.c +++ b/utils/test-prism.c @@ -1,5 +1,27 @@ +/* libctl: flexible Guile-based control files for scientific software + * Copyright (C) 1998-2014 Massachusetts Institute of Technology and Steven G. Johnson + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + * + * Steven G. Johnson can be contacted at stevenj@alum.mit.edu. + */ + /************************************************************************/ /* test-prism.c: unit test for prisms in libctlgeom */ +/* homer reid 5/2018 */ /************************************************************************/ #include #include @@ -25,47 +47,38 @@ static vector3 make_vector3(double x, double y, double z) return v; } -boolean point_in_polygon(double px, double py, vector3 *vertices, int num_vertices); - /************************************************************************/ /* return a uniform random number in [a,b] */ /************************************************************************/ -static double myurand(double a, double b) -{ - return a + (b-a)*(rand()/((double)RAND_MAX)); -} +static double urand(double a, double b) +{ return a + (b-a)*(rand()/((double)RAND_MAX)); } -static double mydrand() -{ return myurand(0.0,1.0); } +static double drand() +{ return urand(0.0,1.0); } /************************************************************************/ -/* return a random unit vector, uniformly distributed over a parallelepiped */ +/* random point uniformly distributed over a parallelepiped */ /************************************************************************/ vector3 random_point(vector3 min_corner, vector3 max_corner) -{ return make_vector3( myurand(min_corner.x, max_corner.x), - myurand(min_corner.y, max_corner.y), - myurand(min_corner.z, max_corner.z) +{ return make_vector3( urand(min_corner.x, max_corner.x), + urand(min_corner.y, max_corner.y), + urand(min_corner.z, max_corner.z) ); } /************************************************************************/ -/* return a random unit vector, uniformly distributed over a sphere */ +/* random unit vector uniformly distributed over a sphere */ /************************************************************************/ vector3 random_unit_vector3(void) { - double z, t, r; - vector3 v; - - z = 2*mydrand() - 1; - t = 2*K_PI*mydrand(); - r = sqrt(1 - z*z); - v.x = r * cos(t); - v.y = r * sin(t); - v.z = z; - return v; + double z = urand(-1.0,1.0); + double rho = sqrt(1 - z*z); + double phi = urand(0.0, 2.0*K_PI); + return make_vector3( rho*cos(phi), rho*sin(phi), z ); } /***************************************************************/ +/* write prism vertices and edges to text file. */ /* after running this routine to produce a file named MyFile, */ /* the prism may be plotted in gnuplot like this: */ /* gnuplot> splot 'MyFile' u 1:2:3 w lp pt 7 ps 1 */ @@ -77,7 +90,8 @@ void prism2gnuplot(prism *prsm, char *filename) double height = prsm->height; FILE *f=fopen(filename,"w"); - for(int nv=0; nv=0.0 ? 1.0 : -1.0; } + +vector3 standardize(vector3 v) +{ vector3 sv=unit_vector3(v); + double sign = (sv.z!=0.0 ? sgn(sv.z) : sv.y!=0.0 ? sgn(sv.y) : sgn(sv.x)); + return vector3_scale(sign,sv); +} + /************************************************************************/ -/* first unit test: create the same parallelepiped in two ways, as a */ -/* block and as a prism, then check inclusion of 10000 */ -/* randomly generated points and confirm the results agree. */ +/* first unit test: check inclusion of randomly-generated points */ /************************************************************************/ -int unit_test1() +int test_point_inclusion(geometric_object the_block, geometric_object the_prism, + int num_tests, int write_log) { - void* m = NULL; - vector3 c = { 0, 0, 0 }; - vector3 xhat = make_vector3(1,0,0); - vector3 yhat = make_vector3(0,1,0); - vector3 zhat = make_vector3(0,0,1); - vector3 size = make_vector3(LX,LY,LZ); - geometric_object the_block = make_block(m, c, xhat, yhat, zhat, size); - - vector3 v[4]; - v[0].x=-0.5*LX; v[0].y=-0.5*LY; v[0].z=-0.5*LZ; - v[1].x=+0.5*LX; v[1].y=-0.5*LY; v[1].z=-0.5*LZ; - v[2].x=+0.5*LX; v[2].y=+0.5*LY; v[2].z=-0.5*LZ; - v[3].x=-0.5*LX; v[3].y=+0.5*LY; v[3].z=-0.5*LZ; - geometric_object the_prism=make_prism(m, v, 4, LZ, zhat); - - int num_failed=0; + vector3 size = the_block.subclass.block_data->size; vector3 min_corner = vector3_scale(-1.0, size); vector3 max_corner = vector3_scale(+1.0, size); - FILE *f=fopen("/tmp/test-prism.points","w"); - for(int n=0; nsize; + vector3 min_corner = vector3_scale(-1.0, size); + vector3 max_corner = vector3_scale(+1.0, size); + FILE *f = write_log ? fopen("/tmp/test-prism.normals","w") : 0; + + int num_failed=0; + int n; + for(n=0; nsize; + vector3 min_corner = vector3_scale(-1.0, size); + vector3 max_corner = vector3_scale(+1.0, size); + FILE *f = write_log ? fopen("/tmp/test-prism.segments","w") : 0; + + int num_failed=0; + int n; + for(n=0; n 1.0e-6*fmax(fabs(sblock),fabs(sprism)) ) + num_failed++; + + if (f) + fprintf(f," %e %e %s\n",sblock,sprism,fabs(sblock-sprism) > 1.0e-6*fmax(fabs(sblock),fabs(sprism)) ? "fail" : "success"); + } + if (f) fclose(f); + + printf("%i/%i segments failed\n",num_failed,num_tests); + return num_failed; +} + +/***************************************************************/ +/* unit tests: create the same parallelepiped two ways (as a */ +/* block and as a prism) and verify that geometric primitives */ +/* give identical results */ +/***************************************************************/ +int run_unit_tests() { void* m = NULL; vector3 c = { 0, 0, 0 }; @@ -184,57 +263,28 @@ int unit_test2() vector3 size = make_vector3(LX,LY,LZ); geometric_object the_block = make_block(m, c, xhat, yhat, zhat, size); - vector3 v[4]; + vector3 v[4]; v[0].x=-0.5*LX; v[0].y=-0.5*LY; v[0].z=-0.5*LZ; v[1].x=+0.5*LX; v[1].y=-0.5*LY; v[1].z=-0.5*LZ; v[2].x=+0.5*LX; v[2].y=+0.5*LY; v[2].z=-0.5*LZ; v[3].x=-0.5*LX; v[3].y=+0.5*LY; v[3].z=-0.5*LZ; geometric_object the_prism=make_prism(m, v, 4, LZ, zhat); - int num_failed=0; - vector3 min_corner = vector3_scale(-1.0, size); - vector3 max_corner = vector3_scale(+1.0, size); - FILE *f=fopen("/tmp/test-prism.lines","w"); - for(int n=0; ns_block[1]) - { double s=s_block[0]; s_block[0]=s_block[1]; s_block[1]=s; } - if (ns_prism==2 && s_prism[0]>s_prism[1]) - { double s=s_prism[0]; s_prism[0]=s_prism[1]; s_prism[1]=s; } - - int match = (ns_block==ns_prism) ? 1 : 0; - for(int n=0; match==1 && n 1.0e-7 ) - match=0; - - if (match==0) - num_failed++; - } - fclose(f); - - printf("%i/%i lines failed\n",num_failed,NUMLINES); - return num_failed; + int write_log=0; + + if (write_log) + prism2gnuplot(the_prism.subclass.prism_data, "/tmp/test-prism.prism"); + + int num_failed_1 = test_point_inclusion(the_block, the_prism, NUMPTS, write_log); + int num_failed_2 = test_normal_to_object(the_block, the_prism, NUMLINES, write_log); + int num_failed_3 = test_line_segment_intersection(the_block, the_prism, NUMLINES, write_log); + + return num_failed_1 + num_failed_2 + num_failed_3; } +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ void print_usage(char *msg, int print_usage) { if (!msg) @@ -246,7 +296,9 @@ void print_usage(char *msg, int print_usage) printf(" --axis x y z\n"); printf("\n"); printf(" --point x y z\n"); - printf(" --pointfile MyPointFile\n"); + printf(" --dir x y z\n"); + printf(" --a a\n"); + printf(" --b b\n"); }; exit(1); } @@ -266,23 +318,19 @@ int main(int argc, char *argv[]) geom_initialize(); if (argc<=1) // if no arguments, run unit tests - { - int num_failed_1 = unit_test1(); - int num_failed_2 = unit_test2(); - return num_failed_1 + num_failed_2; - } + return run_unit_tests(); /***************************************************************/ /* process arguments *******************************************/ /***************************************************************/ - vector3 *vertices=0; - int num_vertices=0; - char *vertexfile=0, *pointfile=0; + char *vertexfile=0; vector3 zhat={0,0,1}; double height=1.5; - vector3 test_point={0,0,0}; int have_test_point=0; - vector3 test_dir={0,0,0}; int have_test_dir=0; - for(int narg=1; narg Date: Wed, 16 May 2018 16:29:24 -0400 Subject: [PATCH 08/10] finished implementing all geometry primitives for prisms; added passing unit test utils/test-prism.c --- utils/prism.c | 2 -- utils/test-prism.c | 2 +- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/utils/prism.c b/utils/prism.c index 3462333..1aec5e0 100644 --- a/utils/prism.c +++ b/utils/prism.c @@ -237,8 +237,6 @@ int insert_s_in_list(double s, double slist[2], int length) /***************************************************************/ int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double slist[2]) { - if (bvertices.num_items; vector3 *vertices = prsm->vertices.items; double height = prsm->height; diff --git a/utils/test-prism.c b/utils/test-prism.c index 8be8964..296c6ad 100644 --- a/utils/test-prism.c +++ b/utils/test-prism.c @@ -324,7 +324,7 @@ int main(int argc, char *argv[]) /* process arguments *******************************************/ /***************************************************************/ char *vertexfile=0; - vector3 zhat={0,0,1}; + vector3 axis={0,0,1}; double height=1.5; vector3 test_point={0,0,0}; vector3 test_dir={0,0,1}; From c7bf25cdb4228fd57200372e1cf1ba3662bbbe69 Mon Sep 17 00:00:00 2001 From: Homer Reid Date: Wed, 16 May 2018 21:53:11 -0400 Subject: [PATCH 09/10] inserted content of prism.c into geom.c, made updates to enable mpb and meep builds --- examples/Makefile.am | 5 +- utils/Makefile.am | 2 +- utils/ctlgeom.h | 13 -- utils/geom.c | 459 ++++++++++++++++++++++++++++++++++++++- utils/prism.c | 500 ------------------------------------------- utils/test-prism.c | 5 +- 6 files changed, 464 insertions(+), 520 deletions(-) delete mode 100644 utils/prism.c diff --git a/examples/Makefile.am b/examples/Makefile.am index fbf28fd..7eddb07 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -57,7 +57,7 @@ CTL_DEFS = -DCTL_SCM='"'$(LIBCTL_DIR)/base/ctl.scm'"' \ -DVERSION_STRING='"'$(VERSION_STRING)'"' example_SOURCES = $(MY_SOURCES) -nodist_example_SOURCES = main.c geom.c prism.c ctl-io.h ctl-io.c +nodist_example_SOURCES = main.c geom.c ctl-io.h ctl-io.c BUILT_SOURCES = $(nodist_example_SOURCES) example_LDADD = $(MY_LIBS) $(LIBCTL) example_LDFLAGS = $(MY_LDFLAGS) @@ -69,9 +69,6 @@ main.c: $(LIBCTL_DIR)/base/main.c geom.c: $(LIBCTL_DIR)/utils/geom.c cp -f $(LIBCTL_DIR)/utils/geom.c $@ -prism.c: $(LIBCTL_DIR)/utils/prism.c - cp -f $(LIBCTL_DIR)/utils/prism.c $@ - ctl-io.c: $(SPECIFICATION_FILE) $(LIBCTL_DIR)/utils/geom.scm $(GEN_CTL_IO) $(GEN_CTL_IO) --code -o $@ $(SPECIFICATION_FILE) $(LIBCTL_DIR) diff --git a/utils/Makefile.am b/utils/Makefile.am index 5a88d53..ad2aef0 100644 --- a/utils/Makefile.am +++ b/utils/Makefile.am @@ -6,7 +6,7 @@ noinst_PROGRAMS = geomtst EXTRA_DIST = gen-ctl-io.in README geom.scm geom-ctl-io-defaults.c nlopt.c -libctlgeom_la_SOURCES = geom.c $(top_srcdir)/src/ctl-math.c $(top_srcdir)/src/integrator.c geom-ctl-io.c prism.c ctlgeom-types.h +libctlgeom_la_SOURCES = geom.c $(top_srcdir)/src/ctl-math.c $(top_srcdir)/src/integrator.c geom-ctl-io.c ctlgeom-types.h libctlgeom_la_LDFLAGS = -version-info @SHARED_VERSION_INFO@ libctlgeom_la_CPPFLAGS = -DLIBCTLGEOM -I$(top_srcdir)/src diff --git a/utils/ctlgeom.h b/utils/ctlgeom.h index f5f76d9..e204a4f 100644 --- a/utils/ctlgeom.h +++ b/utils/ctlgeom.h @@ -152,19 +152,6 @@ GEOMETRIC_OBJECT make_ellipsoid(MATERIAL_TYPE material, vector3 center, GEOMETRIC_OBJECT make_prism(MATERIAL_TYPE material, const vector3 *vertices, int num_vertices, double height, vector3 axis); - -vector3 prism_coordinate_p2c(prism *prsm, vector3 vp); -vector3 prism_coordinate_c2p(prism *prsm, vector3 vc); -vector3 prism_vector_p2c(prism *prsm, vector3 vp); -vector3 prism_vector_c2p(prism *prsm, vector3 vc); -void get_prism_bounding_box(prism *prsm, geom_box *box); -double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b); -boolean point_in_polygon(double px, double py, vector3 *vertices, int num_vertices); -boolean point_in_prism(prism *prsm, vector3 p); -vector3 normal_to_prism(prism *prsm, vector3 p); -void display_prism_info(int indentby, prism *prsm); -int vector3_nearly_equal(vector3 v1, vector3 v2); - /**************************************************************************/ #ifdef __cplusplus diff --git a/utils/geom.c b/utils/geom.c index e9d89b1..3b7ef0c 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -65,6 +65,17 @@ using namespace ctlio; #define MAX(a,b) ((a) > (b) ? (a) : (b)) #define MIN(a,b) ((a) < (b) ? (a) : (b)) +// forward declarations of prism-related routines, at the bottom of this file +vector3 prism_vector_p2c(prism *prsm, vector3 vp); +vector3 prism_vector_c2p(prism *prsm, vector3 vc); +void get_prism_bounding_box(prism *prsm, geom_box *box); +double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b); +boolean point_in_polygon(double px, double py, vector3 *vertices, int num_vertices); +boolean point_in_prism(prism *prsm, vector3 p); +vector3 normal_to_prism(prism *prsm, vector3 p); +void display_prism_info(int indentby, prism *prsm); + + /**************************************************************************/ /* If v is a vector in the lattice basis, normalize v so that @@ -130,7 +141,6 @@ void geom_fix_object(geometric_object o) } case GEOM PRISM: { - // FIXME not really sure what to do here } break; case GEOM COMPOUND_GEOMETRIC_OBJECT: @@ -2059,3 +2069,450 @@ geometric_object make_ellipsoid(material_type material, vector3 center, = 2.0 / size.z; return o; } + +/*************************************************************** + * 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 + * vertices, extruded through a given thickness (the "height") in the + * direction of a given 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. + * 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_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_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; + + // set x,y coordinates of low, high to bounding box + // of prism base polygon (in prism coordinate system) + vector3 low=vertices[0], high=vertices[0]; + int nv; + for(nv=1; nvlow = prism_coordinate_p2c(prsm, low); + box->high = prism_coordinate_p2c(prsm, high); +} + +/***************************************************************/ +/* find the value of s at which the line p+s*d intersects the */ +/* line segment connecting v1 to v2 (in 2 dimensions) */ +/* 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 */ +/***************************************************************/ +boolean intersect_line_with_segment(double px, double py, double dx, double dy, + vector3 v1, vector3 v2, 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; + double DetM = M00*M11 - M01*M10; + double L2 = bx*bx + by*by; // squared length of edge + if ( fabs(DetM) < 1.0e-10*L2 ) // d zero or nearly parallel to edge-->no intersection + return 0; + + double t = (M00*RHSy-M10*RHSx)/DetM; + if (s) *s = (M11*RHSx-M01*RHSy)/DetM; + + if (t<0.0 || t>1.0) // intersection of lines does not lie between vertices + return 0; + + return 1; +} + +// 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) +{ + double ss; + int status=intersect_line_with_segment(px,py,dx,dy,v1,v2,&ss); + if (status==0 || ss<0.0) + return 0; + if (s) *s=ss; + return 1; +} + +/***************************************************************/ +/* 2D point-in-polygon test: return 1 if the point 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; */ +/* point lies in polygon iff this is number is odd. */ +/***************************************************************/ +boolean point_in_polygon(double px, double py, vector3 *vertices, int num_vertices) +{ + double dx=0.0, dy=-1.0; + int num_side_intersections=0; + int nv; + for(nv=0; nv=slist[0]) + slist[1]=s; + else + { slist[1]=slist[0]; slist[0]=s; } + break; + default: + if (svertices.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) + return 0; + return point_in_polygon(xp.x, xp.y, vertices, num_vertices); +} + +/***************************************************************/ +/* insert s into slist, which is sorted in increasing order */ +/* has maximum length 2 */ +/***************************************************************/ +int insert_s_in_list(double s, double slist[2], int length) +{ + if (length==0) + { slist[0]=s; + return 1; + } + if ( s < slist[0] ) + { slist[1]=slist[0]; + slist[0]=s; + } + else if (length==1 || s < slist[1] ) + slist[1]=s; + + return 2; +} + +/***************************************************************/ +/* compute all values of s at which the line p+s*d intersects */ +/* a prism face, retaining the first two entries in a list */ +/* of s values sorted in increasing order m */ +/***************************************************************/ +int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double slist[2]) +{ + 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); + + // use length of first edge as a length scale for judging + // lengths to be small or large + double length_scale = vector3_norm(vector3_minus(vertices[1], vertices[0])); + + // identify intersections with prism side faces + int num_intersections=0; + int nv; + for(nv=0; nvz_max) ) + continue; + + num_intersections=insert_s_in_list(s, slist, num_intersections); + } + + // 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)) + 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)) + continue; + num_intersections=insert_s_in_list(s,slist,num_intersections); + } + return num_intersections; +} + +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ +double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b) +{ + double slist[2]={0.0, 0.0}; + if (2!=intersect_line_with_prism(prsm, p, d, slist)) + return 0.0; + double ds = fmin(slist[1],b) - fmax(slist[0],a); + return ds > 0.0 ? ds : 0.0; +} + +/***************************************************************/ +/* compute the minimum-length vector from p to the plane */ +/* containing o (origin) and spanned by basis vectors v1,v2 */ +/* algorithm: solve the 3x3 system p-s*v3 = o + t*v1 + u*v2 */ +/* where v3 = v1 x v2 and s,t,u are unknowns */ +/***************************************************************/ +vector3 normal_to_plane(vector3 o, vector3 v1, vector3 v2, vector3 p) +{ + vector3 RHS = vector3_minus(p,o); + + // handle the degenerate-to-2D case + if ( (fabs(v1.z) + fabs(v2.z)) < 2.0e-7 && fabs(RHS.z)<1.0e-7 ) + { vector3 zhat={0,0,1}; + vector3 v3 = vector3_cross(zhat, v1); + double M00 = v1.x; double M10 = v3.x; + double M01 = v1.y; double M11 = v3.y; + double DetM = M00*M11 - M01*M10; + if ( fabs(DetM) < 1.0e-10 ) + return vector3_scale(0.0, v3); + // double t = (M00*RHSy-M10*RHSx)/DetM; + double s= (M11*RHS.x-M01*RHS.y)/DetM; + return vector3_scale(-1.0*s,v3); + } + + vector3 v3 = vector3_cross(v1, v2); + matrix3x3 M; + M.c0 = v1; + M.c1 = v2; + M.c2 = vector3_scale(-1.0, v3); + vector3 tus = matrix3x3_vector3_mult(matrix3x3_inverse(M),RHS); + return vector3_scale(-1.0*tus.z, v3); +} + +/***************************************************************/ +/* find the face of the prism for which the normal distance */ +/* from x 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 centroid = prsm->centroid; + double height = prsm->height; + vector3 *vertices = prsm->vertices.items; + int num_vertices = prsm->vertices.num_items; + + vector3 xp = prism_coordinate_c2p(prsm, xc); + vector3 axis = {0,0,0}; axis.z=height; + + vector3 retval; + double min_distance; + + // side walls + int nv; + for(nv=0; nvvertices.items; + int num_vertices = prsm->vertices.num_items; + double height = prsm->height; + vector3 z0 = {0.0, 0.0, 1.0}; + vector3 axis = matrix3x3_vector3_mult(prsm->m_c2p,z0); + + 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); + int nv; + for(nv=0; nv=3, "fewer than 3 vertices in make_prism"); + + // compute centroid of vertices + vector3 centroid = {0.0, 0.0, 0.0}; + int nv; + for(nv=0; nvcentroid = centroid; + prsm->height = height; + prsm->m_p2c = m_p2c; + prsm->m_c2p = m_c2p; + + // 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)); + for(nv=0; nvvertices.items[nv] = prism_coordinate_c2p(prsm,vertices[nv]); + + // note the center and centroid are different! + vector3 center = vector3_plus(centroid, vector3_scale(0.5*height,zhat) ); + geometric_object o=make_geometric_object(material, center); + o.which_subclass=GEOM PRISM; + o.subclass.prism_data = prsm; + + return o; +} diff --git a/utils/prism.c b/utils/prism.c deleted file mode 100644 index 1aec5e0..0000000 --- a/utils/prism.c +++ /dev/null @@ -1,500 +0,0 @@ -/* libctl: flexible Guile-based control files for scientific software - * Copyright (C) 1998-2014 Massachusetts Institute of Technology and Steven G. Johnson - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, write to the - * Free Software Foundation, Inc., 59 Temple Place - Suite 330, - * Boston, MA 02111-1307, USA. - * - * Steven G. Johnson can be contacted at stevenj@alum.mit.edu. - */ - -/* - * prism.c -- geometry routines 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", - * in which the prism floor lies in the XY plane with centroid - * at the origin and the prism axis is the positive Z-axis. - * homer reid 4/2018 - */ - -#include -#include -#include - -/*--------------------------------------------------------------*/ -/* repeat some definitions from geom.c that are needed here too;*/ -/* can these be moved to ctlgeom.h? */ -/*--------------------------------------------------------------*/ -#ifndef LIBCTLGEOM -# include "ctl-io.h" -#else -# define material_type void* - static void material_type_copy(void **src, void **dest) { *dest = *src; } -#endif -#include - -#ifdef CXX_CTL_IO - using namespace ctlio; -# define GEOM geometric_object:: -#else -# define GEOM -#endif - -#ifdef __cplusplus -# define MALLOC(type,num) (new type[num]) -# define MALLOC1(type) (new type) -# define FREE(p) delete[] (p) -# define FREE1(p) delete (p) -#else -# define MALLOC(type,num) ((type *) malloc(sizeof(type) * (num))) -# define MALLOC1(type) MALLOC(type,1) -# define FREE(p) free(p) -# define FREE1(p) free(p) -#endif - -#define CHECK(cond, s) if (!(cond)){fprintf(stderr,s "\n");exit(EXIT_FAILURE);} - -/***************************************************************/ -/* 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_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_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; - - // set x,y coordinates of low, high to bounding box - // of prism base polygon (in prism coordinate system) - vector3 low=vertices[0], high=vertices[0]; - int nv; - for(nv=1; nvlow = prism_coordinate_p2c(prsm, low); - box->high = prism_coordinate_p2c(prsm, high); -} - -/***************************************************************/ -/* find the value of s at which the line p+s*d intersects the */ -/* line segment connecting v1 to v2 (in 2 dimensions) */ -/* 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 */ -/***************************************************************/ -boolean intersect_line_with_segment(double px, double py, double dx, double dy, - vector3 v1, vector3 v2, 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; - double DetM = M00*M11 - M01*M10; - double L2 = bx*bx + by*by; // squared length of edge - if ( fabs(DetM) < 1.0e-10*L2 ) // d zero or nearly parallel to edge-->no intersection - return 0; - - double t = (M00*RHSy-M10*RHSx)/DetM; - if (s) *s = (M11*RHSx-M01*RHSy)/DetM; - - if (t<0.0 || t>1.0) // intersection of lines does not lie between vertices - return 0; - - return 1; -} - -// 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) -{ - double ss; - int status=intersect_line_with_segment(px,py,dx,dy,v1,v2,&ss); - if (status==0 || ss<0.0) - return 0; - if (s) *s=ss; - return 1; -} - -/***************************************************************/ -/* 2D point-in-polygon test: return 1 if the point 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; */ -/* point lies in polygon iff this is number is odd. */ -/***************************************************************/ -boolean point_in_polygon(double px, double py, vector3 *vertices, int num_vertices) -{ - double dx=0.0, dy=-1.0; - int num_side_intersections=0; - int nv; - for(nv=0; nv=slist[0]) - slist[1]=s; - else - { slist[1]=slist[0]; slist[0]=s; } - break; - default: - if (svertices.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) - return 0; - return point_in_polygon(xp.x, xp.y, vertices, num_vertices); -} - -/***************************************************************/ -/* insert s into slist, which is sorted in increasing order */ -/* has maximum length 2 */ -/***************************************************************/ -int insert_s_in_list(double s, double slist[2], int length) -{ - if (length==0) - { slist[0]=s; - return 1; - } - if ( s < slist[0] ) - { slist[1]=slist[0]; - slist[0]=s; - } - else if (length==1 || s < slist[1] ) - slist[1]=s; - - return 2; -} - -/***************************************************************/ -/* compute all values of s at which the line p+s*d intersects */ -/* a prism face, retaining the first two entries in a list */ -/* of s values sorted in increasing order m */ -/***************************************************************/ -int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double slist[2]) -{ - 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); - - // use length of first edge as a length scale for judging - // lengths to be small or large - double length_scale = vector3_norm(vector3_minus(vertices[1], vertices[0])); - - // identify intersections with prism side faces - int num_intersections=0; - int nv; - for(nv=0; nvz_max) ) - continue; - - num_intersections=insert_s_in_list(s, slist, num_intersections); - } - - // 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)) - 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)) - continue; - num_intersections=insert_s_in_list(s,slist,num_intersections); - } - return num_intersections; -} - -/***************************************************************/ -/***************************************************************/ -/***************************************************************/ -double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b) -{ - double slist[2]={0.0, 0.0}; - if (2!=intersect_line_with_prism(prsm, p, d, slist)) - return 0.0; - double ds = fmin(slist[1],b) - fmax(slist[0],a); - return ds > 0.0 ? ds : 0.0; -} - -/***************************************************************/ -/* compute the minimum-length vector from p to the plane */ -/* containing o (origin) and spanned by basis vectors v1,v2 */ -/* algorithm: solve the 3x3 system p-s*v3 = o + t*v1 + u*v2 */ -/* where v3 = v1 x v2 and s,t,u are unknowns */ -/***************************************************************/ -vector3 normal_to_plane(vector3 o, vector3 v1, vector3 v2, vector3 p) -{ - vector3 RHS = vector3_minus(p,o); - - // handle the degenerate-to-2D case - if ( (fabs(v1.z) + fabs(v2.z)) < 2.0e-7 && fabs(RHS.z)<1.0e-7 ) - { vector3 zhat={0,0,1}; - vector3 v3 = vector3_cross(zhat, v1); - double M00 = v1.x; double M10 = v3.x; - double M01 = v1.y; double M11 = v3.y; - double DetM = M00*M11 - M01*M10; - if ( fabs(DetM) < 1.0e-10 ) - return vector3_scale(0.0, v3); - // double t = (M00*RHSy-M10*RHSx)/DetM; - double s= (M11*RHS.x-M01*RHS.y)/DetM; - return vector3_scale(-1.0*s,v3); - } - - vector3 v3 = vector3_cross(v1, v2); - matrix3x3 M; - M.c0 = v1; - M.c1 = v2; - M.c2 = vector3_scale(-1.0, v3); - vector3 tus = matrix3x3_vector3_mult(matrix3x3_inverse(M),RHS); - return vector3_scale(-1.0*tus.z, v3); -} - -/***************************************************************/ -/* find the face of the prism for which the normal distance */ -/* from x 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 centroid = prsm->centroid; - double height = prsm->height; - vector3 *vertices = prsm->vertices.items; - int num_vertices = prsm->vertices.num_items; - - vector3 xp = prism_coordinate_c2p(prsm, xc); - vector3 axis = {0,0,0}; axis.z=height; - - vector3 retval; - double min_distance; - - // side walls - int nv; - for(nv=0; nvvertices.items; - int num_vertices = prsm->vertices.num_items; - double height = prsm->height; - vector3 z0 = {0.0, 0.0, 1.0}; - vector3 axis = matrix3x3_vector3_mult(prsm->m_c2p,z0); - - 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); - int nv; - for(nv=0; nv=3, "fewer than 3 vertices in make_prism"); - - // compute centroid of vertices - vector3 centroid = {0.0, 0.0, 0.0}; - int nv; - for(nv=0; nvcentroid = centroid; - prsm->height = height; - prsm->m_p2c = m_p2c; - prsm->m_c2p = m_c2p; - - // 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)); - for(nv=0; nvvertices.items[nv] = prism_coordinate_c2p(prsm,vertices[nv]); - - // note the center and centroid are different! - vector3 center = vector3_plus(centroid, vector3_scale(0.5*height,zhat) ); - geometric_object o=make_geometric_object(material, center); - o.which_subclass=GEOM PRISM; - o.subclass.prism_data = prsm; - - return o; -} diff --git a/utils/test-prism.c b/utils/test-prism.c index 296c6ad..704102c 100644 --- a/utils/test-prism.c +++ b/utils/test-prism.c @@ -40,6 +40,9 @@ #define LY 1.0 #define LZ 1.5 +// routine defined in geom.c +int vector3_nearly_equal(vector3 v1, vector3 v2); + static vector3 make_vector3(double x, double y, double z) { vector3 v; @@ -404,7 +407,7 @@ int main(int argc, char *argv[]) printf("point {%e,%e,%e}: \n",test_point.x,test_point.y,test_point.z); printf(" %s prism\n", in_prism ? "in" : "not in"); printf(" normal to prism: {%e,%e,%e}\n",nhat.x,nhat.y,nhat.z); - printf(" intersection with line segment {%e,%e,%e} + (%e,%e)*{%e,%e,%e}: %s\n", + printf(" intersection with line segment {%e,%e,%e} + (%e,%e)*{%e,%e,%e}: %e\n", test_point.x, test_point.y, test_point.z, a,b,test_dir.x, test_dir.y, test_dir.z,s); } From e22d1fd937acbe807cb5346358e055eb936472f3 Mon Sep 17 00:00:00 2001 From: Homer Reid Date: Wed, 16 May 2018 22:10:44 -0400 Subject: [PATCH 10/10] fixed compilation issues for test-prism --- utils/ctlgeom.h | 5 +++-- utils/geom.c | 27 +++++++++++++-------------- utils/test-prism.c | 10 ++++++---- 3 files changed, 22 insertions(+), 20 deletions(-) diff --git a/utils/ctlgeom.h b/utils/ctlgeom.h index e204a4f..d963e45 100644 --- a/utils/ctlgeom.h +++ b/utils/ctlgeom.h @@ -147,11 +147,12 @@ GEOMETRIC_OBJECT make_block(MATERIAL_TYPE material, vector3 center, GEOMETRIC_OBJECT make_ellipsoid(MATERIAL_TYPE material, vector3 center, vector3 e1, vector3 e2, vector3 e3, vector3 size); - -// routines in prism.c GEOMETRIC_OBJECT make_prism(MATERIAL_TYPE material, const vector3 *vertices, int num_vertices, double height, vector3 axis); + +int vector3_nearly_equal(vector3 v1, vector3 v2, double tolerance); + /**************************************************************************/ #ifdef __cplusplus diff --git a/utils/geom.c b/utils/geom.c index 3b7ef0c..5fd6562 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -66,16 +66,14 @@ using namespace ctlio; #define MIN(a,b) ((a) < (b) ? (a) : (b)) // forward declarations of prism-related routines, at the bottom of this file -vector3 prism_vector_p2c(prism *prsm, vector3 vp); -vector3 prism_vector_c2p(prism *prsm, vector3 vc); -void get_prism_bounding_box(prism *prsm, geom_box *box); -double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b); -boolean point_in_polygon(double px, double py, vector3 *vertices, int num_vertices); -boolean point_in_prism(prism *prsm, vector3 p); -vector3 normal_to_prism(prism *prsm, vector3 p); -void display_prism_info(int indentby, prism *prsm); - - +static vector3 prism_vector_p2c(prism *prsm, vector3 vp); +static vector3 prism_vector_c2p(prism *prsm, vector3 vc); +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); /**************************************************************************/ /* If v is a vector in the lattice basis, normalize v so that @@ -2433,9 +2431,9 @@ void display_prism_info(int indentby, prism *prsm) /***************************************************************/ // like vector3_equal but tolerant of small floating-point discrepancies /***************************************************************/ -int vector3_nearly_equal(vector3 v1, vector3 v2) +int vector3_nearly_equal(vector3 v1, vector3 v2, double tolerance) { - return (vector3_norm( vector3_minus(v1,v2) ) < 1.0e-7*vector3_norm(v1)); + return (vector3_norm( vector3_minus(v1,v2) ) <= tolerance*vector3_norm(v1)); } /***************************************************************/ @@ -2468,12 +2466,13 @@ GEOMETRIC_OBJECT make_prism(MATERIAL_TYPE material, // make sure all vertices lie in a plane normal to the given axis vector3 zhat = unit_vector3(axis); + double tolerance=1.0e-6; for(nv=0; nv