diff --git a/AUTHORS b/AUTHORS index 996550d..6f57dbb 100644 --- a/AUTHORS +++ b/AUTHORS @@ -3,6 +3,7 @@ libctl was written by Steven G. Johnson (stevenj@alum.mit.edu) with contribution M.T. Homer Reid Christopher Hogan Ardavan Oskooi +Daniel W. Boyce The multidimensional integration routines in src/integrator.c were adapted from the HIntlib Library by Rudolf Schuerer and from the GNU Scientific diff --git a/utils/ctlgeom.h b/utils/ctlgeom.h index 33349f4..0bf19ee 100644 --- a/utils/ctlgeom.h +++ b/utils/ctlgeom.h @@ -164,6 +164,15 @@ GEOMETRIC_OBJECT make_prism_with_center(MATERIAL_TYPE material, vector3 center, const vector3 *vertices, int num_vertices, double height, vector3 axis); +// slanted prism with `center` field computed automatically from vertices, height, axis, sidewall_angle +GEOMETRIC_OBJECT make_slanted_prism(MATERIAL_TYPE material, const vector3 *vertices_bottom, int num_vertices, + double height, vector3 axis, double sidewall_angle); + +// as make_slanted_prism, but with a rigid translation so that the prism is centered at center +GEOMETRIC_OBJECT make_slanted_prism_with_center(MATERIAL_TYPE material, vector3 center, + const vector3 *vertices_bottom, int num_vertices, double height, + vector3 axis, double sidewall_angle); + int vector3_nearly_equal(vector3 v1, vector3 v2, double tolerance); /**************************************************************************/ diff --git a/utils/geom.c b/utils/geom.c index 579e77a..544039d 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -77,6 +77,7 @@ static boolean point_in_prism(prism *prsm, vector3 pc); static vector3 normal_to_prism(prism *prsm, vector3 pc); static double intersect_line_segment_with_prism(prism *prsm, vector3 pc, vector3 dc, double a, double b); +static double get_prism_volume(prism *prsm); static void get_prism_bounding_box(prism *prsm, geom_box *box); static void display_prism_info(int indentby, geometric_object *o); static void init_prism(geometric_object *o); @@ -913,15 +914,7 @@ double geom_object_volume(GEOMETRIC_OBJECT o) { return o.subclass.block_data->which_subclass == BLK BLOCK_SELF ? vol : vol * (K_PI / 6); } case GEOM PRISM: { - vector3_list vertices = o.subclass.prism_data->vertices_p; - double area = 0; - int i; - for (i = 0; i < vertices.num_items; ++i) { - int i1 = (i + 1) % vertices.num_items; - area += 0.5 * (vertices.items[i1].x - vertices.items[i].x) * - (vertices.items[i1].y + vertices.items[i].y); - } - return fabs(area) * o.subclass.prism_data->height; + return get_prism_volume(o.subclass.prism_data); } default: return 0; /* unsupported object types? */ } @@ -1903,8 +1896,10 @@ geometric_object make_ellipsoid(material_type material, vector3 center, vector3 /*************************************************************** * The remainder of this file implements geometric primitives for prisms. * A prism is a planar polygon, consisting of 3 or more user-specified - * vertices, extruded through a given thickness (the "height") in the - * direction of a given unit vector (the "axis.") + * vertices (the "bottom_vertices), extruded through a given thickness + * (the "height") in the direction of a given unit vector (the "axis") + * with the walls of the extrusion tapering at a given angle angle + * (the "sidewall_angle). * 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. @@ -1975,6 +1970,7 @@ int intersect_line_with_segment(vector3 q0, vector3 q1, vector3 q2, vector3 u, d double DetM = M00 * M11 - M01 * M10; double L2 = M01 * M01 + M11 * M11; // squared length of edge, used to set length scale if (fabs(DetM) < 1.0e-10 * L2) { // d zero or nearly parallel to edge + if (vector3_nearly_equal(q0, q1, 1e-12) || vector3_nearly_equal(q0, q2, 1e-12)) return IN_SEGMENT; double q01x = q0.x - q1.x, q01y = q0.y - q1.y, q01 = sqrt(q01x * q01x + q01y * q01y); double q02x = q0.x - q2.x, q02y = q0.y - q2.y, q02 = sqrt(q02x * q02x + q02y * q02y); double dot = q01x * q02x + q01y * q02y; @@ -2042,7 +2038,7 @@ boolean node_in_or_on_polygon(vector3 q0, vector3 *nodes, int num_nodes, // Is q0 on a vertex or edge? // Transform coordinate system of nodes such that q0 is at 0|0 for (nn = 0; nn < num_nodes; nn++) { - int status = intersect_ray_with_segment(q0, nodes[nn], nodes[(nn + 1) % num_nodes], xAxis, 0); + int status = intersect_ray_with_segment(q0, nodes[nn], nodes[(nn + 1) % num_nodes], unit_vector3(vector3_minus(nodes[(nn + 1) % num_nodes], nodes[nn])), 0); if (status == IN_SEGMENT) { return include_boundaries; } // Find start point which is not on the x axis (from q0) @@ -2114,8 +2110,12 @@ boolean point_in_or_on_prism(prism *prsm, vector3 pc, boolean include_boundaries double height = prsm->height; vector3 pp = prism_coordinate_c2p(prsm, pc); if (pp.z < 0.0 || pp.z > prsm->height) return 0; - vector3 *nodes = prsm->vertices_p.items; - int num_nodes = prsm->vertices_p.num_items; + int num_nodes = prsm->vertices_bottom_p.num_items; + vector3 nodes[num_nodes]; + int nv; + for (nv = 0; nv < num_nodes; nv++) { + nodes[nv] = vector3_plus(prsm->vertices_bottom_p.items[nv], vector3_scale(pp.z, prsm->top_polygon_diff_vectors_scaled_p.items[nv])); + } return node_in_or_on_polygon(pp, nodes, num_nodes, include_boundaries); } @@ -2150,38 +2150,38 @@ static int dcmp(const void *pd1, const void *pd2) { int intersect_line_with_prism(prism *prsm, vector3 pc, vector3 dc, double *slist) { vector3 pp = prism_coordinate_c2p(prsm, pc); vector3 dp = prism_vector_c2p(prsm, dc); - vector3 *vps = prsm->vertices_p.items; - int num_vertices = prsm->vertices_p.num_items; + vector3 *vps_bottom = prsm->vertices_bottom_p.items; + vector3 *vps_top = prsm->vertices_top_p.items; + int num_vertices = prsm->vertices_bottom_p.num_items; double height = prsm->height; - // use length of first polygon edge as a length scale for judging - // lengths to be small or large - double length_scale = vector3_norm(vector3_minus(vps[1], vps[0])); - // identify intersections with prism side faces + double tus_tolerance = 1e-8; int num_intersections = 0; int nv; for (nv = 0; nv < num_vertices; nv++) { int nvp1 = nv + 1; if (nvp1 == num_vertices) nvp1 = 0; - // first solve the in-plane portion of the problem: determine the - // intersection of the XY-plane projection of the line with the - // polygon edge between vertices (nv,nv+1). - double s; - int status = intersect_line_with_segment(pp, vps[nv], vps[nvp1], dp, &s); - if (status == NON_INTERSECTING || status == ON_RAY) continue; - - // OK, we know the XY-plane projection of the line intersects the polygon edge; - // now go back to 3D, ask for the z-coordinate at which the line intersects - // the z-axis extrusion of the edge, and determine whether this point lies - // inside or outside the height of the prism. - double zp_intersect = pp.z + s * dp.z; - double AbsTol = 1.0e-7 * length_scale; - double zp_min = 0.0 - AbsTol; - double zp_max = height + AbsTol; - if ((zp_intersect < zp_min) || (zp_intersect > zp_max)) continue; - + // checks if dp is parallel to the plane of the prism side face under consideration + vector3 v1 = vector3_minus(vps_bottom[nvp1], vps_bottom[nv]); + vector3 v2 = vector3_minus(vps_top[nv], vps_bottom[nv]); + double dot_tolerance = 1e-6; + if (fabs(vector3_dot(dp, vector3_cross(v1, v2))) <= dot_tolerance) continue; + + // to find the intersection point pp + s*dp between the line and the + // prism side face, we will solve the vector equation + // pp + s*dp = o + t*v1 + u*v2 + // where o is vps_bottom[nv], v1 is vps_bottom[nvp1]-vps_bottom[nv], + // v2 is vps_top[nv]-vps_bottom[nv], and 0 <= t <= 1, 0 <= u <= 1. + matrix3x3 M; + M.c0 = v1; + M.c1 = v2; + M.c2 = vector3_scale(-1, dp); + vector3 RHS = vector3_minus(pp, vps_bottom[nv]); + vector3 tus = matrix3x3_vector3_mult(matrix3x3_inverse(M), RHS); + if (tus.x < -tus_tolerance || tus.x > 1+tus_tolerance || tus.y < -tus_tolerance || tus.y > 1+tus_tolerance) continue; + double s = tus.z; slist[num_intersections++] = s; } @@ -2191,12 +2191,30 @@ int intersect_line_with_prism(prism *prsm, vector3 pc, vector3 dc, double *slist for (LowerUpper = 0; LowerUpper < 2; LowerUpper++) { double z0p = LowerUpper ? height : 0.0; double s = (z0p - pp.z) / dp.z; + vector3 *vps = LowerUpper ? vps_top : vps_bottom; if (!node_in_polygon(pp.x + s * dp.x, pp.y + s * dp.y, vps, num_vertices)) continue; slist[num_intersections++] = s; } qsort((void *)slist, num_intersections, sizeof(double), dcmp); - return num_intersections; + // if num_intersections is zero then just return that + if (num_intersections == 0) return num_intersections; + else { + // remove duplicates from slist + double duplicate_tolerance = 1e-3; + int num_unique_elements = 1; + double slist_unique[num_vertices+2]; + slist_unique[0] = slist[0]; + for (nv = 1; nv < num_intersections; nv++) { + if (fabs(slist[nv] - slist[nv-1]) > duplicate_tolerance*fabs(slist[nv])) { + slist_unique[num_unique_elements] = slist[nv]; + num_unique_elements++; + } + } + slist = slist_unique; + num_intersections = num_unique_elements; + return num_intersections; + } } /***************************************************************/ @@ -2296,10 +2314,19 @@ double min_distance_to_quadrilateral(vector3 p, vector3 o, vector3 v1, vector3 v // fc==0/1 for floor/ceiling double min_distance_to_prism_roof_or_ceiling(vector3 pp, prism *prsm, int fc) { - vector3 *vps = prsm->vertices_p.items; - int num_vertices = prsm->vertices_p.num_items; - vector3 op = {0.0, 0.0, 0.0}; - if (fc == 1) op.z = prsm->height; // origin of floor/ceiling + int num_vertices = prsm->vertices_bottom_p.num_items; + vector3 op = {0.0, 0.0, 0.0}; // origin of floor/ceiling + vector3 vps[num_vertices]; + if (fc == 1) { + memcpy(vps, prsm->vertices_top_p.items, num_vertices * sizeof(vector3)); + for (int i = 0; i < num_vertices; i++) { + vps[i].z = 0; + } + op.z = prsm->height; + } + else { + memcpy(vps, prsm->vertices_bottom_p.items, num_vertices * sizeof(vector3)); + } vector3 zhatp = {0, 0, 1.0}; double s = normal_distance_to_plane(pp, op, vps[0], vps[1], zhatp, 0); vector3 ppProj = @@ -2322,8 +2349,9 @@ vector3 normal_to_prism(prism *prsm, vector3 pc) { if (prsm->height == 0.0) return prsm->axis; double height = prsm->height; - vector3 *vps = prsm->vertices_p.items; - int num_vertices = prsm->vertices_p.num_items; + vector3 *vps_bottom = prsm->vertices_bottom_p.items; + vector3 *vps_diff_to_top = prsm->top_polygon_diff_vectors_p.items; + int num_vertices = prsm->vertices_bottom_p.num_items; vector3 zhatp = {0.0, 0.0, 1.0}; vector3 axisp = vector3_scale(height, zhatp); @@ -2335,9 +2363,9 @@ vector3 normal_to_prism(prism *prsm, vector3 pc) { // consider side walls for (nv = 0; nv < num_vertices; nv++) { int nvp1 = (nv == (num_vertices - 1) ? 0 : nv + 1); - vector3 v0p = vps[nv]; - vector3 v1p = vector3_minus(vps[nvp1], vps[nv]); - vector3 v2p = axisp; + vector3 v0p = vps_bottom[nv]; + vector3 v1p = vector3_minus(vps_bottom[nvp1], vps_bottom[nv]); + vector3 v2p = vps_diff_to_top[nv]; vector3 v3p = unit_vector3(vector3_cross(v1p, v2p)); double s = min_distance_to_quadrilateral(pp, v0p, v1p, v2p, v3p); if (fabs(s) < min_distance) { @@ -2357,19 +2385,125 @@ vector3 normal_to_prism(prism *prsm, vector3 pc) { return prism_vector_p2c(prsm, retval); } +/***************************************************************/ +/* Compute the area of a polygon using its vertices. */ +/***************************************************************/ +double get_area_of_polygon_from_nodes(vector3 *nodes, int num_nodes){ + double area = 0; + int i; + for (i = 0; i < num_nodes; ++i) { + int i1 = (i + 1) % num_nodes; + area += 0.5 * (nodes[i1].x - nodes[i].x) * + (nodes[i1].y + nodes[i].y); + } + return fabs(area); +} + +/***************************************************************/ +/* This computes the volume of an irregular triangular prism */ +/* following the scheme of http://darrenirvine.blogspot.com/2011/10/volume-of-irregular-triangular-prism.html */ +/* The two end triangles have points a, b, and c. Angle abc */ +/* is right, and angles bac and acb are acute, of course. The */ +/* primary constraint for this method is that the lines be- */ +/* tween 'a's between triangles, betweens 'b's between */ +/* triangles, and between 'c's between triangles must all be */ +/* parallel, though the end triangles need not be parallel. */ + +/***************************************************************/ +double get_volume_irregular_triangular_prism(vector3 a0, vector3 b0, vector3 c0, vector3 a1, vector3 b1, vector3 c1) { + vector3 side_a = vector3_minus(a1, a0); + vector3 side_b = vector3_minus(b1, b0); + vector3 side_c = vector3_minus(c1, c0); + if (vector3_norm(vector3_cross(side_a, side_b)) != 0 || + vector3_norm(vector3_cross(side_b, side_c)) != 0 || + vector3_norm(vector3_cross(side_c, side_a)) != 0) { + //throw an error + } + double length_side_a = vector3_norm(side_a); + double length_side_b = vector3_norm(side_b); + double length_side_c = vector3_norm(side_c); + double average_length = (length_side_a + length_side_b + length_side_c) / 3.0; + + vector3 point_on_plane = a0; + vector3 plane_normal_vector = unit_vector3(side_a); + vector3 plane_to_a1 = vector3_minus(a1, point_on_plane); + vector3 plane_to_b1 = vector3_minus(b1, point_on_plane); + vector3 plane_to_c1 = vector3_minus(c1, point_on_plane); + vector3 proj_plane_to_a1_on_normal_vector = vector3_scale(vector3_dot(plane_normal_vector, plane_to_a1), plane_normal_vector); + vector3 proj_plane_to_b1_on_normal_vector = vector3_scale(vector3_dot(plane_normal_vector, plane_to_b1), plane_normal_vector); + vector3 proj_plane_to_c1_on_normal_vector = vector3_scale(vector3_dot(plane_normal_vector, plane_to_c1), plane_normal_vector); + vector3 a1_on_plane = vector3_minus(a1, proj_plane_to_a1_on_normal_vector); + vector3 b1_on_plane = vector3_minus(b1, proj_plane_to_b1_on_normal_vector); + vector3 c1_on_plane = vector3_minus(c1, proj_plane_to_c1_on_normal_vector); + double base = vector3_norm(vector3_minus(c1_on_plane, b1_on_plane)); + double height = vector3_norm(vector3_minus(a1_on_plane, b1_on_plane)); + double cross_sectional_area = fabs(0.5 * base * height); + + return fabs(average_length * cross_sectional_area); +} + +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ +double get_prism_volume(prism *prsm) { + if (prsm->sidewall_angle == 0.0) { + return get_area_of_polygon_from_nodes(prsm->vertices_bottom_p.items, prsm->vertices_bottom_p.num_items) * fabs(prsm->height); + } + else { + int num_vertices = prsm->vertices_bottom_p.num_items; + double bottom_polygon_area = get_area_of_polygon_from_nodes(prsm->vertices_bottom_p.items, prsm->vertices_bottom_p.num_items); + double top_polygon_area = get_area_of_polygon_from_nodes(prsm->vertices_top_p.items, prsm->vertices_top_p.num_items); + double volume; + vector3 *wedges_a; + wedges_a = (vector3 *)malloc(num_vertices * sizeof(vector3)); + CHECK(wedges_a, "out of memory"); + vector3 *wedges_b; + wedges_b = (vector3 *)malloc(num_vertices * sizeof(vector3)); + CHECK(wedges_b, "out of memory"); + vector3 *wedges_c; + wedges_c = (vector3 *)malloc(num_vertices * sizeof(vector3)); + CHECK(wedges_c, "out of memory"); + if (bottom_polygon_area > top_polygon_area) { + volume = fabs(top_polygon_area * prsm->height); + memcpy(wedges_a, prsm->vertices_top_p.items, num_vertices * sizeof(vector3)); + memcpy(wedges_b, prsm->vertices_top_p.items, num_vertices * sizeof(vector3)); + for (int nv = 0; nv < num_vertices; nv++) { + wedges_b[nv].z = 0.0; + } + memcpy(wedges_c, prsm->vertices_bottom_p.items, num_vertices * sizeof(vector3)); + } + else { + volume = fabs(bottom_polygon_area * prsm->height); + memcpy(wedges_a, prsm->vertices_bottom_p.items, num_vertices * sizeof(vector3)); + memcpy(wedges_b, prsm->vertices_bottom_p.items, num_vertices * sizeof(vector3)); + for (int nv = 0; nv < num_vertices; nv++) { + wedges_b[nv].z = prsm->height; + } + memcpy(wedges_c, prsm->vertices_top_p.items, num_vertices * sizeof(vector3)); + } + for (int nv = 0; nv < num_vertices; nv++) { + int nvp1 = (nv + 1 == num_vertices ? 0 : nv + 1); + volume += get_volume_irregular_triangular_prism(wedges_a[nv], wedges_b[nv], wedges_c[nv], wedges_a[nvp1], wedges_b[nvp1], wedges_c[nvp1]); + } + return volume; + } +} + /***************************************************************/ /***************************************************************/ /***************************************************************/ void get_prism_bounding_box(prism *prsm, geom_box *box) { - vector3 *vertices = prsm->vertices.items; - int num_vertices = prsm->vertices.num_items; - box->low = box->high = vertices[0]; + vector3 *vertices_bottom = prsm->vertices_bottom.items; + vector3 *vertices_top = prsm->vertices_top.items; + int num_vertices = prsm->vertices_bottom.num_items; + box->low = box->high = vertices_bottom[0]; int nv, fc; for (nv = 0; nv < num_vertices; nv++) for (fc = 0; fc < 2; fc++) // 'floor,ceiling' { - vector3 v = vertices[nv]; - if (fc == 1) v = vector3_plus(v, vector3_scale(prsm->height, prsm->axis)); + vector3 v; + if (fc == 0) v = vertices_bottom[nv]; + if (fc == 1) v = vertices_top[nv]; box->low.x = fmin(box->low.x, v.x); box->low.y = fmin(box->low.y, v.y); @@ -2387,11 +2521,11 @@ void get_prism_bounding_box(prism *prsm, geom_box *box) { void display_prism_info(int indentby, geometric_object *o) { prism *prsm = o->subclass.prism_data; - vector3 *vs = prsm->vertices.items; - int num_vertices = prsm->vertices.num_items; + vector3 *vs = prsm->vertices_bottom.items; + int num_vertices = prsm->vertices_bottom.num_items; - ctl_printf("%*s height %g, axis (%g,%g,%g), %i vertices:\n", indentby, "", prsm->height, - prsm->axis.x, prsm->axis.y, prsm->axis.z, num_vertices); + ctl_printf("%*s height %g, axis (%g,%g,%g), sidewall angle: %g radians, %i vertices:\n", indentby, "", prsm->height, + prsm->axis.x, prsm->axis.y, prsm->axis.z, prsm->sidewall_angle, num_vertices); int nv; for (nv = 0; nv < num_vertices; nv++) ctl_printf("%*s (%g,%g,%g)\n", indentby, "", vs[nv].x, vs[nv].y, vs[nv].z); @@ -2403,6 +2537,7 @@ void display_prism_info(int indentby, geometric_object *o) { int vector3_nearly_equal(vector3 v1, vector3 v2, double tolerance) { return (vector3_norm(vector3_minus(v1, v2)) <= tolerance * vector3_norm(v1)); } +matrix3x3 sidewall_scaling_matrix; /***************************************************************/ /* return the unit normal vector to the triangle with the given*/ @@ -2417,28 +2552,29 @@ vector3 triangle_normal(vector3 v1, vector3 v2, vector3 v3) { /***************************************************************/ /* On entry, the only fields in o->prism that are assumed to */ -/* be initialized are: vertices, height, and (optionally) axis.*/ -/* If axis has not been initialized (i.e. it is set to its */ -/* default value, which is the zero vector) then the prism axis*/ -/* is automatically computed as the normal to the vertex plane.*/ -/* If o->center is equal to auto_center on entry, then it is */ -/* set to the prism center, as computed from the vertices, */ -/* axis, and height. Otherwise, the prism is rigidly translated*/ -/* to center it at the specified value of o->center. */ +/* be initialized are: vertices_bottom, height, (optionally) */ +/* axis, and sidewall_angle. If axis has not been initialized */ +/* (i.e. it is set to its default value, which is the zero */ +/* vector) then the prism axis is automatically computed as */ +/* the normal to the vertex plane. If o->center is equal to */ +/* auto_center on entry, then it is set to the prism center, */ +/* as computed from the vertices, axis, and height. Otherwise, */ +/* the prism is rigidly translated to center it at the */ +/* specified value of o->center. */ /***************************************************************/ // special vector3 that signifies 'no value specified' vector3 auto_center = {NAN, NAN, NAN}; void init_prism(geometric_object *o) { prism *prsm = o->subclass.prism_data; - vector3 *vertices = prsm->vertices.items; - int num_vertices = prsm->vertices.num_items; + vector3 *vertices_bottom = prsm->vertices_bottom.items; + int num_vertices = prsm->vertices_bottom.num_items; CHECK(num_vertices >= 3, "fewer than 3 vertices in init_prism"); // compute centroid of vertices vector3 centroid = {0.0, 0.0, 0.0}; int nv; for (nv = 0; nv < num_vertices; nv++) - centroid = vector3_plus(centroid, vertices[nv]); + centroid = vector3_plus(centroid, vertices_bottom[nv]); prsm->centroid = centroid = vector3_scale(1.0 / ((double)num_vertices), centroid); // make sure all vertices lie in a plane, i.e. that the normal @@ -2448,7 +2584,7 @@ void init_prism(geometric_object *o) { double tol = 1.0e-6; for (nv = 0; nv < num_vertices; nv++) { int nvp1 = (nv + 1) % num_vertices; - vector3 tri_normal = triangle_normal(centroid, vertices[nv], vertices[nvp1]); + vector3 tri_normal = triangle_normal(centroid, vertices_bottom[nv], vertices_bottom[nvp1]); if (vector3_norm(tri_normal) == 0.0) // vertices collinear with centroid continue; if (!plane_normal_set) { @@ -2487,7 +2623,7 @@ void init_prism(geometric_object *o) { else { vector3 shift = vector3_minus(o->center, current_center); for (nv = 0; nv < num_vertices; nv++) - vertices[nv] = vector3_plus(vertices[nv], shift); + vertices_bottom[nv] = vector3_plus(vertices_bottom[nv], shift); centroid = vector3_plus(centroid, shift); } @@ -2516,18 +2652,158 @@ void init_prism(geometric_object *o) { yhat = y0hat; } else { - xhat = unit_vector3(vector3_minus(vertices[1], vertices[0])); + xhat = unit_vector3(vector3_minus(vertices_bottom[1], vertices_bottom[0])); yhat = unit_vector3(vector3_cross(zhat, xhat)); } matrix3x3 m_p2c = {xhat, yhat, zhat}; prsm->m_p2c = m_p2c; prsm->m_c2p = matrix3x3_inverse(m_p2c); - // compute vertices in prism coordinate system - prsm->vertices_p.num_items = num_vertices; - prsm->vertices_p.items = (vector3 *)malloc(num_vertices * sizeof(vector3)); + // compute vertices_bottom in prism coordinate system + prsm->vertices_bottom_p.num_items = num_vertices; + prsm->vertices_bottom_p.items = (vector3 *)malloc(num_vertices * sizeof(vector3)); for (nv = 0; nv < num_vertices; nv++) - prsm->vertices_p.items[nv] = prism_coordinate_c2p(prsm, vertices[nv]); + prsm->vertices_bottom_p.items[nv] = prism_coordinate_c2p(prsm, vertices_bottom[nv]); + + // Calculate difference vertices of the top polygon and vectors between bottom + // polygon and the top polygon, where: + // * the bottom polygon is the one passed in to the the make_prism() function, + // stored in vertices_bottom and vertices_bottom_p, + // * the top polygon is the top surface (parallel to the bottom polygon) resulting + // from the extrusion of the bottom polygon. Whether or not the extrusion tapers + // is dependent on the value of sidewall_angle. + // + // The top polygon is calculated by first copying the values of vertices_bottom_p into + // vertices_top_p, except z=prsm->height for all top vertices. If prsm->sidewall_angle + // is equal to zero, then no further calculations are performed on the top vertices. + // If not, we know that all EDGES of the the top polygon will be offset so that in the + // xy plane they are parallel to the edges of the bottom polygon. The offset amount is + // determined by the sidewall angle and the height of the prism. To perform the + // calculation, each of the edges of the top polygon (without an offset) are stored in + // an array of edges (edge is a struct defined if prsm->sidewall_angle!=0 containing + // the endpoints a1 a2, with a third vector v defined a2-a1). Then the vector normal to + // v is calculated, and the offset vector. A test is performed to determine in which + // direction (the direction of +offset or -offset) from the edge we can find points + // inside the polygon by performing a node_in_or_on_polygon test at a finite distance + // away from the midpoint of the edge: + // edge.a1 + 0.5*edge.v + 1e-3*offset. + // This information is used to determine in which direction the offset of the edge is + // applied, in conjunction with whether prsm->sidewall_angle is positive or negative + // (if positive, the offset will be applied in towards the points where + // node_in_or_on_polygon is true, else the offset will be applied out away from those + // points). After the offsets are applied to the edges, the intersections between the + // new edges are calculated, which are the new values of vertices_top_p. + // + // Some side notes on the difference vectors: + // * The value of each of the top polygon vertices can be found + // vertices_bottom_p + top_polygon_diff_vectors_p + // vertices_bottom + top_polygon_diff_vectors + // * A linearly interpolated value of the polygon vertices between the bottom + // polygon and the top can be found + // vertices_bottom_p + top_polygon_diff_vectors_scaled_p * z + number theta = (K_PI/2) - fabs(prsm->sidewall_angle); + prsm->vertices_top_p.num_items = num_vertices; + prsm->vertices_top_p.items = (vector3 *)malloc(num_vertices * sizeof(vector3)); + CHECK(prsm->vertices_top_p.items, "out of memory"); + memcpy(prsm->vertices_top_p.items, prsm->vertices_bottom_p.items, num_vertices * sizeof(vector3)); + for (nv = 0; nv < num_vertices; nv++) { + prsm->vertices_top_p.items[nv].z = prsm->height; + } + + if (prsm->sidewall_angle != 0.0) { + typedef struct { + vector3 a1, a2, v; // v will be defined as a2 - a1 + } edge; + + // find the point at the bottom left corner of the polygon + double smallest_x = HUGE_VAL; + double smallest_y = HUGE_VAL; + int index_for_point_a = -1; + int index_for_point_b = -1; + int index_for_point_c = -1; + for (nv = 0; nv < num_vertices; nv++) { + double current_x = prsm->vertices_bottom_p.items[nv].x; + double current_y = prsm->vertices_bottom_p.items[nv].y; + if (current_x < smallest_x) { + smallest_x = current_x; + smallest_y = current_y; + index_for_point_b = nv; + } + else if (current_x == smallest_x && current_y < smallest_y) { + smallest_y = current_y; + index_for_point_b = nv; + } + } + if (index_for_point_b == -1) { + exit(EXIT_FAILURE); + } + else { + index_for_point_a = (index_for_point_b + 1 == num_vertices ? 0 : index_for_point_b + 1); + index_for_point_c = (index_for_point_b - 1 == -1 ? num_vertices - 1 : index_for_point_b - 1); + } + // find orientation of the polygon + vector3 A = prsm->vertices_bottom_p.items[index_for_point_a]; + vector3 B = prsm->vertices_bottom_p.items[index_for_point_b]; + vector3 C = prsm->vertices_bottom_p.items[index_for_point_c]; + double orientation_number = (B.x - A.x)*(C.y - A.y)-(C.x - A.x)*(B.y - A.y); + int orientation_positive_or_negative = (orientation_number < 0 ? 0 : 1); + + edge *top_polygon_edges; + top_polygon_edges = (edge *)malloc(num_vertices * sizeof(edge)); + number w = prsm->height / tan(theta); + + for (nv = 0; nv < num_vertices; nv++) { + top_polygon_edges[nv].a1 = prsm->vertices_top_p.items[(nv - 1 == -1 ? num_vertices - 1 : nv - 1)]; + top_polygon_edges[nv].a2 = prsm->vertices_top_p.items[nv]; + top_polygon_edges[nv].v = vector3_minus(top_polygon_edges[nv].a2, top_polygon_edges[nv].a1); + + vector3 normal_vector = (orientation_positive_or_negative ? unit_vector3(vector3_cross(top_polygon_edges[nv].v, zhat)) : unit_vector3(vector3_cross(top_polygon_edges[nv].v, vector3_scale(-1, zhat)))); + + // positive sidewall angles means the prism tapers in towards the rest of the prism body + // negative sidewall angles means the prism tapers out away from the rest of the prism body + vector3 offset = vector3_scale(prsm->sidewall_angle > 0 ? w : -w, normal_vector); + top_polygon_edges[nv].a1 = vector3_plus(top_polygon_edges[nv].a1, offset); + top_polygon_edges[nv].a2 = vector3_plus(top_polygon_edges[nv].a2, offset); + } + + for (nv = 0; nv < num_vertices; nv++) { + number x1 = top_polygon_edges[nv].a1.x; + number y1 = top_polygon_edges[nv].a1.y; + number x2 = top_polygon_edges[nv].a2.x; + number y2 = top_polygon_edges[nv].a2.y; + number x3 = top_polygon_edges[(nv + 1 == num_vertices ? 0 : nv + 1)].a1.x; + number y3 = top_polygon_edges[(nv + 1 == num_vertices ? 0 : nv + 1)].a1.y; + number x4 = top_polygon_edges[(nv + 1 == num_vertices ? 0 : nv + 1)].a2.x; + number y4 = top_polygon_edges[(nv + 1 == num_vertices ? 0 : nv + 1)].a2.y; + + // Intersection point calculated with https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection#Given_two_points_on_each_line + number px = ((x1*y2-y1*x2)*(x3-x4)-(x1-x2)*(x3*y4-y3*x4)) / ((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4)); + number py = ((x1*y2-y1*x2)*(y3-y4)-(y1-y2)*(x3*y4-y3*x4)) / ((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4)); + prsm->vertices_top_p.items[nv].x = px; + prsm->vertices_top_p.items[nv].y = py; + } + } + + prsm->top_polygon_diff_vectors_p.num_items = num_vertices; + prsm->top_polygon_diff_vectors_p.items = (vector3 *)malloc(num_vertices * sizeof(vector3)); + CHECK(prsm->top_polygon_diff_vectors_p.items, "out of memory"); + for (nv = 0; nv < num_vertices; nv++) { + prsm->top_polygon_diff_vectors_p.items[nv] = vector3_minus(prsm->vertices_top_p.items[nv], prsm->vertices_bottom_p.items[nv]); + } + + prsm->top_polygon_diff_vectors_scaled_p.num_items = num_vertices; + prsm->top_polygon_diff_vectors_scaled_p.items = (vector3 *)malloc(num_vertices * sizeof(vector3)); + CHECK(prsm->top_polygon_diff_vectors_scaled_p.items, "out of memory"); + for (nv = 0; nv < num_vertices; nv++) { + prsm->top_polygon_diff_vectors_scaled_p.items[nv] = vector3_scale(1/prsm->height, prsm->top_polygon_diff_vectors_p.items[nv]); + } + + prsm->vertices_top.num_items = num_vertices; + prsm->vertices_top.items = (vector3 *)malloc(num_vertices * sizeof(vector3)); + CHECK(prsm->vertices_top.items, "out of memory"); + for (nv = 0; nv < num_vertices; nv++) { + prsm->vertices_top.items[nv] = prism_coordinate_p2c(prsm, prsm->vertices_top_p.items[nv]); + } // workspace is an internally-stored double-valued array of length num_vertices+2 // that is used by some geometry routines @@ -2544,20 +2820,34 @@ geometric_object make_prism(material_type material, const vector3 *vertices, int return make_prism_with_center(material, auto_center, vertices, num_vertices, height, axis); } -// prism in which all vertices are translated to ensure that the prism is centered at center -geometric_object make_prism_with_center(material_type material, vector3 center, - const vector3 *vertices, int num_vertices, double height, - vector3 axis) { +// prism in which all vertices are translated to ensure that the prism is centered at center. +geometric_object make_prism_with_center(material_type material, vector3 center, const vector3 *vertices_bottom, + int num_vertices, double height, vector3 axis) { + return make_slanted_prism_with_center(material, center, vertices_bottom, num_vertices, height, axis, 0); +} + +// slanted prism with center determined automatically from vertices, height, axis, and sidewall_angle +geometric_object make_slanted_prism(material_type material, const vector3 *vertices_bottom, + int num_vertices, double height, vector3 axis, double sidewall_angle) { + return make_slanted_prism_with_center(material, auto_center, vertices_bottom, num_vertices, height, axis, sidewall_angle); +} + +// Have both make_prism_with_center and make_slanted_prism_with_center keep the same parameters to maintain ABI +// compatibility, though make_prism_with_center just calls make_slanted_prism_with_center with the sidewall angle equal +// to zero. To make a slanted prism, the user will have to call make_slanted_prism for now. +geometric_object make_slanted_prism_with_center(material_type material, vector3 center, const vector3 *vertices_bottom, + int num_vertices, double height, vector3 axis, double sidewall_angle) { geometric_object o = make_geometric_object(material, center); o.which_subclass = GEOM PRISM; prism *prsm = o.subclass.prism_data = MALLOC1(prism); CHECK(prsm, "out of memory"); - prsm->vertices.num_items = num_vertices; - prsm->vertices.items = (vector3 *)malloc(num_vertices * sizeof(vector3)); - CHECK(prsm->vertices.items, "out of memory"); - memcpy(prsm->vertices.items, vertices, num_vertices * sizeof(vector3)); + prsm->vertices_bottom.num_items = num_vertices; + prsm->vertices_bottom.items = (vector3 *)malloc(num_vertices * sizeof(vector3)); + CHECK(prsm->vertices_bottom.items, "out of memory"); + memcpy(prsm->vertices_bottom.items, vertices_bottom, num_vertices * sizeof(vector3)); prsm->height = height; prsm->axis = axis; + prsm->sidewall_angle = sidewall_angle; init_prism(&o); return o; } diff --git a/utils/geom.scm b/utils/geom.scm index ef9f6ba..24c6495 100644 --- a/utils/geom.scm +++ b/utils/geom.scm @@ -92,13 +92,13 @@ ; some notes regarding prisms: ; (a) When instantiating a prism, typically only the -; fields `vertices`, `height,` and (optionally) `axis` -; will be initialized by the user; all remaining fields are -; derived properties that are computed internally. (So, morally -; they should be thought of as having been declared using -; `define-derived-property` or `define-post-processed-property,` -; except here the code that does the derivation or -; post-processing is implemented in C, not scheme.) +; fields `vertices_bottom`, `height,` (optionally) `axis`, +; and `sidewall_angle` will be initialized by the user; all +; remaining fields are derived properties that are computed +; internally. (So, morally, they should be thought of as having +; been declared using `define-derived-property` or +; `define-post-processed-property,` except here the code that does +; the derivation or post-processing is implemented in C, not scheme.) ; (b) The suffix _p (for "prism") is used to identify variables ; that store coordinates of points or components of vectors ; in the prism coordinate system. (The prism coordinate system @@ -117,13 +117,23 @@ ; (center = centroid + 0.5*height*axis), so---in contrast to all other ; types of geometric-object---there is no need to specify the `center` ; field when instantiating a prism. +; (f) The sidwall angle determines an angle at which the prism is extruded. +; A positive sidewall angle determines a prism that extrudes inward at +; the given angle, and a negative sidewall angle determines a prisms +; that extrudes outward. This is useful for modeling a prism formed in +; a foundry that cannot grow objects with a perfectly normal sidewall. (define-class prism geometric-object ; fields to be filled in by users - (define-property vertices '() (make-list-type 'vector3)) + (define-property vertices_bottom '() (make-list-type 'vector3)) (define-property height 0 'number) (define-property axis (vector3 0 0 0) 'vector3) + (define-property sidewall_angle 0 'number) ; derived fields computed internally - (define-property vertices_p '() (make-list-type 'vector3)) + (define-property vertices_bottom_p '() (make-list-type 'vector3)) + (define-property top_polygon_diff_vectors_p '() (make-list-type 'vector3)) + (define-property top_polygon_diff_vectors_scaled_p '() (make-list-type 'vector3)) + (define-property vertices_top_p '() (make-list-type 'vector3)) + (define-property vertices_top '() (make-list-type 'vector3)) (define-property centroid (vector3 0 0 0) 'vector3) (define-property workspace '() (make-list-type 'number)) (define-property m_c2p identity_matrix 'matrix3x3) diff --git a/utils/test-prism.c b/utils/test-prism.c index 7188f9e..56c8b53 100644 --- a/utils/test-prism.c +++ b/utils/test-prism.c @@ -79,16 +79,16 @@ void GPQuad(FILE *f, vector3 v, vector3 l1, vector3 l2, prism *prsm) { /***************************************************************/ /***************************************************************/ void my_get_prism_bounding_box(prism *prsm, geom_box *box) { - vector3 *vertices = prsm->vertices_p.items; - int num_vertices = prsm->vertices_p.num_items; + vector3 *vertices_bottom = prsm->vertices_bottom_p.items; + int num_vertices = prsm->vertices_bottom_p.num_items; double height = prsm->height; - box->low = box->high = prism_coordinate_p2c(prsm, vertices[0]); + box->low = box->high = prism_coordinate_p2c(prsm, vertices_bottom[0]); int nv, fc; for (nv = 0; nv < num_vertices; nv++) for (fc = 0; fc < 2; fc++) // 'floor,ceiling' { - vector3 vp = vertices[nv]; + vector3 vp = vertices_bottom[nv]; if (fc == 1) vp.z = height; vector3 vc = prism_coordinate_p2c(prsm, vp); @@ -129,13 +129,15 @@ vector3 random_point_in_box(vector3 min_corner, vector3 max_corner) { /* random point uniformly distributed over a planar polygon */ /* (all z coordinates are 0) */ /************************************************************************/ -vector3 random_point_in_polygon(vector3 *vertices, int num_vertices) { +vector3 random_point_in_polygon(prism *prsm) { // randomly choose a vertex and generate random point within the triangle // formed by that vertex, the next vertex, and the centroid + vector3 *vertices_bottom = prsm->vertices_bottom.items; + int num_vertices = prsm->vertices_bottom.num_items; int which_vertex = rand() % num_vertices; vector3 v0 = {0, 0, 0}; - vector3 v1 = vertices[which_vertex]; - vector3 v2 = vertices[(which_vertex + 1) % num_vertices]; + vector3 v1 = vertices_bottom[which_vertex]; + vector3 v2 = vertices_bottom[(which_vertex + 1) % num_vertices]; double xi = urand(0.0, 1.0), eta = urand(0.0, 1.0 - xi); return vector3_plus(vector3_scale(xi, vector3_minus(v1, v0)), vector3_scale(eta, vector3_minus(v2, v0))); @@ -146,8 +148,8 @@ vector3 random_point_in_polygon(vector3 *vertices, int num_vertices) { /************************************************************************/ vector3 random_point_on_prism(geometric_object o) { prism *prsm = o.subclass.prism_data; - vector3 *vertices = prsm->vertices_p.items; - int num_vertices = prsm->vertices_p.num_items; + vector3 *vertices_bottom = prsm->vertices_bottom_p.items; + int num_vertices = prsm->vertices_bottom_p.num_items; double height = prsm->height; // choose a face @@ -155,15 +157,15 @@ vector3 random_point_on_prism(geometric_object o) { int which_face = rand() % num_faces; if (which_face < num_vertices) // side face { - vector3 min_corner = vertices[which_face]; - vector3 max_corner = vertices[(which_face + 1) % num_vertices]; + vector3 min_corner = vertices_bottom[which_face]; + vector3 max_corner = vertices_bottom[(which_face + 1) % num_vertices]; max_corner.z = height; return random_point_in_box(prism_coordinate_p2c(prsm, min_corner), prism_coordinate_p2c(prsm, max_corner)); } else // floor or ceiling { - vector3 p = random_point_in_polygon(vertices, num_vertices); + vector3 p = random_point_in_polygon(prsm); if (which_face == num_faces - 1) p.z = height; return prism_coordinate_p2c(prsm, p); } @@ -185,31 +187,40 @@ vector3 random_unit_vector3() { /* gnuplot> splot 'MyFile' u 1:2:3 w lp pt 7 ps 1 */ /***************************************************************/ void prism2gnuplot(prism *prsm, char *filename) { - vector3 *vertices = prsm->vertices_p.items; - int num_vertices = prsm->vertices_p.num_items; + int num_vertices = prsm->vertices_bottom_p.num_items; double height = prsm->height; + vector3_list vertices_bottom; + vertices_bottom.num_items = num_vertices; + vertices_bottom.items = (vector3 *)malloc(num_vertices * sizeof(vector3)); + memcpy(vertices_bottom.items, prsm->vertices_bottom_p.items, num_vertices * sizeof(vector3)); + vector3_list vertices_top; + vertices_top.num_items = num_vertices; + vertices_top.items = (vector3 *)malloc(num_vertices * sizeof(vector3)); + int nv; + for (nv = 0; nv < num_vertices; nv++) { + vertices_top.items[nv] = vector3_plus(prsm->vertices_bottom_p.items[nv], prsm->top_polygon_diff_vectors_p.items[nv]); + } FILE *f = fopen(filename, "w"); - int nv; for (nv = 0; nv < num_vertices; nv++) { - vector3 vap = vertices[nv]; + vector3 vap = vertices_bottom.items[nv]; vap.z = 0.0; - vector3 vbp = vertices[nv]; + vector3 vbp = vertices_top.items[nv]; vbp.z = height; - vector3 vcp = vertices[(nv + 1) % num_vertices]; + vector3 vcp = vertices_top.items[(nv + 1) % num_vertices]; vcp.z = height; - vector3 vdp = vertices[(nv + 1) % num_vertices]; + vector3 vdp = vertices_bottom.items[(nv + 1) % num_vertices]; vdp.z = 0.0; vector3 vac = prism_coordinate_p2c(prsm, vap); vector3 vbc = prism_coordinate_p2c(prsm, vbp); vector3 vcc = prism_coordinate_p2c(prsm, vcp); vector3 vdc = prism_coordinate_p2c(prsm, vdp); - fprintf(f, "%e %e %e %e %e %e \n", vac.x, vac.y, vac.z, vap.x, vap.y, vap.z); - fprintf(f, "%e %e %e %e %e %e \n", vbc.x, vbc.y, vbc.z, vbp.x, vbp.y, vbp.z); - fprintf(f, "%e %e %e %e %e %e \n", vcc.x, vcc.y, vcc.z, vcp.x, vcp.y, vcp.z); - fprintf(f, "%e %e %e %e %e %e \n", vdc.x, vdc.y, vdc.z, vdp.x, vdp.y, vdp.z); - fprintf(f, "%e %e %e %e %e %e \n", vac.x, vac.y, vac.z, vap.x, vap.y, vap.z); + fprintf(f, "%0.16e %0.16e %0.16e \n", vac.x, vac.y, vac.z); + fprintf(f, "%0.16e %0.16e %0.16e \n", vbc.x, vbc.y, vbc.z); + fprintf(f, "%0.16e %0.16e %0.16e \n", vcc.x, vcc.y, vcc.z); + fprintf(f, "%0.16e %0.16e %0.16e \n", vdc.x, vdc.y, vdc.z); + fprintf(f, "%0.16e %0.16e %0.16e \n", vac.x, vac.y, vac.z); fprintf(f, "\n\n"); } fclose(f); @@ -219,8 +230,8 @@ void prism2gnuplot(prism *prsm, char *filename) { /* write prism vertices and edges to GMSH geometry (.geo) file */ /***************************************************************/ void prism2gmsh(prism *prsm, char *filename) { - vector3 *vertices = prsm->vertices_p.items; - int num_vertices = prsm->vertices_p.num_items; + vector3 *vertices_bottom = prsm->vertices_bottom_p.items; + int num_vertices = prsm->vertices_bottom_p.num_items; double height = prsm->height; vector3 zhat = prsm->m_p2c.c2; vector3 axis = vector3_scale(height, zhat); @@ -228,7 +239,7 @@ void prism2gmsh(prism *prsm, char *filename) { FILE *f = fopen(filename, "w"); int nv; for (nv = 0; nv < num_vertices; nv++) { - vector3 vp = vertices[nv]; + vector3 vp = vertices_bottom[nv]; vector3 vc = prism_coordinate_p2c(prsm, vp); fprintf(f, "Point(%i)={%e, %e, %e};\n", nv, vc.x, vc.y, vc.z); } @@ -254,7 +265,7 @@ vector3 standardize(vector3 v) { } /************************************************************************/ -/* first unit test: check inclusion of randomly-generated points */ +/* 1st unit test: check inclusion of randomly-generated points */ /************************************************************************/ int test_point_inclusion(geometric_object the_block, geometric_object the_prism, int num_tests, int write_log) { @@ -289,7 +300,7 @@ int test_point_inclusion(geometric_object the_block, geometric_object the_prism, } /************************************************************************/ -/* second unit test: check calculation of normals to objects */ +/* 2nd unit test: check calculation of normals to objects */ /************************************************************************/ #define PFACE 0.1 int test_normal_to_object(geometric_object the_block, geometric_object the_prism, int num_tests, @@ -327,7 +338,7 @@ int test_normal_to_object(geometric_object the_block, geometric_object the_prism } /************************************************************************/ -/* third unit test: check-line segment intersections */ +/* 3rd unit test: check-line segment intersections */ /************************************************************************/ int test_line_segment_intersection(geometric_object the_block, geometric_object the_prism, int num_tests, int write_log) { @@ -367,7 +378,7 @@ int test_line_segment_intersection(geometric_object the_block, geometric_object } /************************************************************************/ -/* fourth unit test: check of point in polygon test with slanted H */ +/* 4th unit test: check of point in polygon test with slanted H */ /************************************************************************/ int test_point_in_polygon(int write_log) { // make array of test points that should always pass @@ -449,6 +460,630 @@ int test_point_in_polygon(int write_log) { return num_failed; } +/************************************************************************/ +/* 5th unit test: saves a prism with a square base with a normal */ +/* sidewall angle and a prism with the same base polygon with non- */ +/* normal sidewall angle to separate GNU plot files. */ +/************************************************************************/ +int test_square_base_sidewall_prisms_to_gnuplot() { + void *m = NULL; + + int num_nodes_square = 4; + vector3 nodes_square[num_nodes_square]; + nodes_square[0] = make_vector3(-10.0, -10.0, 0.0); + nodes_square[1] = make_vector3(-10.0, 10.0, 0.0); + nodes_square[2] = make_vector3(10.0, 10.0, 0.0); + nodes_square[3] = make_vector3(10.0, -10.0, 0.0); + + double height_square = 100; + vector3 zhat = make_vector3(0, 0, 1); + + double normal_sidewall = 0; + geometric_object square_normal_sidewall_geom_object = make_prism(m, nodes_square, num_nodes_square, height_square, zhat); + prism *square_normal_sidewall_prism = square_normal_sidewall_geom_object.subclass.prism_data; + + double one_degree_sidewall = 1.0 * 2 * K_PI / 360.0; + geometric_object square_one_degree_sidewall_geom_object = make_slanted_prism(m, nodes_square, num_nodes_square, height_square, zhat, one_degree_sidewall); + prism *square_one_degree_sidewall_prism = square_one_degree_sidewall_geom_object.subclass.prism_data; + + prism2gnuplot(square_normal_sidewall_prism, "square_normal_sidewall_gnu_plot.dat"); + prism2gnuplot(square_one_degree_sidewall_prism, "square_one_degree_sidewall_gnu_plot.dat"); + + return 0; +} + +/************************************************************************/ +/* 6th unit test: saves a prism with a concave octagonal c-shaped */ +/* base with a normal sidewall angle and a prism with the same base */ +/* polygon with non-normal sidewall angle to separate GNU plot files. */ +/************************************************************************/ +int test_octagon_c_base_sidewall_prisms_to_gnuplot() { + void *m = NULL; + + int num_nodes_octagon_c = 16; + vector3 nodes_octagon_c[num_nodes_octagon_c]; + nodes_octagon_c[0] = make_vector3(114.905, 88.7434, 0.0); + nodes_octagon_c[1] = make_vector3(88.7434, 114.905, 0.0); + nodes_octagon_c[2] = make_vector3(51.7447, 114.905, 0.0); + nodes_octagon_c[3] = make_vector3(25.5827, 88.7434, 0.0); + nodes_octagon_c[4] = make_vector3(25.5827, 51.7447, 0.0); + nodes_octagon_c[5] = make_vector3(51.7447, 25.5827, 0.0); + nodes_octagon_c[6] = make_vector3(88.7434, 25.5827, 0.0); + nodes_octagon_c[7] = make_vector3(114.905, 51.7447, 0.0); + nodes_octagon_c[8] = make_vector3(140.488, 41.1477, 0.0); + nodes_octagon_c[9] = make_vector3(99.3401, 0.0, 0.0); + nodes_octagon_c[10] = make_vector3(41.1477, 0.0, 0.0); + nodes_octagon_c[11] = make_vector3(0.0, 41.1477, 0.0); + nodes_octagon_c[12] = make_vector3(0.0, 99.3401, 0.0); + nodes_octagon_c[13] = make_vector3(41.1477, 140.488, 0.0); + nodes_octagon_c[14] = make_vector3(99.3401, 140.488, 0.0); + nodes_octagon_c[15] = make_vector3(140.488, 99.3401, 0.0); + + double height_octagon_c = 127; + vector3 zhat = make_vector3(0, 0, 1); + + double normal_sidewall = 0; + geometric_object octagon_c_normal_sidewall_geom_object = make_prism(m, nodes_octagon_c, num_nodes_octagon_c, height_octagon_c, zhat); + prism *octagon_c_normal_sidewall_prism = octagon_c_normal_sidewall_geom_object.subclass.prism_data; + + double two_half_degree_sidewall = 2.5 * 2 * K_PI / 360.0; + geometric_object octagon_c_two_half_degree_sidewall_geom_object = make_slanted_prism(m, nodes_octagon_c, num_nodes_octagon_c, height_octagon_c, zhat, two_half_degree_sidewall); + prism *octagon_c_two_half_degree_sidewall_prism = octagon_c_two_half_degree_sidewall_geom_object.subclass.prism_data; + + prism2gnuplot(octagon_c_normal_sidewall_prism, "octagon_c_normal_sidewall_gnu_plot.dat"); + prism2gnuplot(octagon_c_two_half_degree_sidewall_prism, "octagon_c_two_half_degree_sidewall_gnu_plot.dat"); + + return 0; +} + +/************************************************************************/ +/* 7th unit test: test all of geom.c's prism helper functions on a */ +/* prism with a concave octagonal c-shaped base with both a normal */ +/* sidewall angle a 2.5-degree sidewall angle. */ +/************************************************************************/ +double relative_error(double actual, double expected) { + return fabs((actual-expected)/actual); +} + +int test_helper_functions_on_octagonal_c_prism() { + double tolerance = 5.0e-5; + + void *m = NULL; + + int num_nodes_octagon_c = 16; + vector3 nodes_octagon_c[num_nodes_octagon_c]; + nodes_octagon_c[0] = make_vector3(114.905, 88.7434, 0.0); + nodes_octagon_c[1] = make_vector3(88.7434, 114.905, 0.0); + nodes_octagon_c[2] = make_vector3(51.7447, 114.905, 0.0); + nodes_octagon_c[3] = make_vector3(25.5827, 88.7434, 0.0); + nodes_octagon_c[4] = make_vector3(25.5827, 51.7447, 0.0); + nodes_octagon_c[5] = make_vector3(51.7447, 25.5827, 0.0); + nodes_octagon_c[6] = make_vector3(88.7434, 25.5827, 0.0); + nodes_octagon_c[7] = make_vector3(114.905, 51.7447, 0.0); + nodes_octagon_c[8] = make_vector3(140.488, 41.1477, 0.0); + nodes_octagon_c[9] = make_vector3(99.3401, 0.0, 0.0); + nodes_octagon_c[10] = make_vector3(41.1477, 0.0, 0.0); + nodes_octagon_c[11] = make_vector3(0.0, 41.1477, 0.0); + nodes_octagon_c[12] = make_vector3(0.0, 99.3401, 0.0); + nodes_octagon_c[13] = make_vector3(41.1477, 140.488, 0.0); + nodes_octagon_c[14] = make_vector3(99.3401, 140.488, 0.0); + nodes_octagon_c[15] = make_vector3(140.488, 99.3401, 0.0); + + double height_octagon_c = 127; + vector3 zhat = make_vector3(0, 0, 1); + + double normal_sidewall = 0; + geometric_object octagon_c_normal_sidewall_geom_object = make_prism(m, nodes_octagon_c, num_nodes_octagon_c, height_octagon_c, zhat); + prism *octagon_c_normal_sidewall_prism = octagon_c_normal_sidewall_geom_object.subclass.prism_data; + + double two_half_degree_sidewall = 2.5 * 2 * K_PI / 360.0; + geometric_object octagon_c_two_half_degree_sidewall_geom_object = make_slanted_prism(m, nodes_octagon_c, num_nodes_octagon_c, height_octagon_c, zhat, two_half_degree_sidewall); + prism *octagon_c_two_half_degree_sidewall_prism = octagon_c_two_half_degree_sidewall_geom_object.subclass.prism_data; + + int num_tests_normal = 0; + int num_failed_normal = 0; + int num_tests_tapered = 0; + int num_failed_tapered = 0; + + printf("prism helper function testing:\n"); + + // test geom_object_volume + double volume_normal_sidewall_freecad = 1082462.27453587; + double volume_normal_sidewall_calculated = geom_object_volume(octagon_c_normal_sidewall_geom_object); + num_tests_normal++; + if (relative_error(volume_normal_sidewall_calculated, volume_normal_sidewall_freecad) > tolerance) { + num_failed_normal++; + } + + double volume_tapered_sidewall_freecad = 833978.754046812; + double volume_tapered_sidewall_calculated = geom_object_volume(octagon_c_two_half_degree_sidewall_geom_object); + num_tests_tapered++; + if (relative_error(volume_tapered_sidewall_calculated, volume_tapered_sidewall_freecad) > tolerance) { + num_failed_tapered++; + } + + // test point_in_prism + vector3_list point_in_prism_test_points_normal_sidewall; + point_in_prism_test_points_normal_sidewall.num_items = 25; + point_in_prism_test_points_normal_sidewall.items = (vector3 *)malloc(point_in_prism_test_points_normal_sidewall.num_items * sizeof(vector3)); + point_in_prism_test_points_normal_sidewall.items[0] = make_vector3(46.4462, 12.7914, 63.5000); // interior point + point_in_prism_test_points_normal_sidewall.items[1] = make_vector3(127.697, 46.4462, 95.2500); // interior point + point_in_prism_test_points_normal_sidewall.items[2] = make_vector3(70.2439, 0.00000, 31.7500); // point on external side face + point_in_prism_test_points_normal_sidewall.items[3] = make_vector3(101.824, 38.6637, 95.2500); // point on internal side face + point_in_prism_test_points_normal_sidewall.items[4] = make_vector3(19.1870, 49.0955, 127.000); // point on top face + point_in_prism_test_points_normal_sidewall.items[5] = make_vector3(134.092, 96.6909, 0.00000); // point on bottom face + point_in_prism_test_points_normal_sidewall.items[6] = make_vector3(127.6965, 94.04175, 127.0); // edge on top + point_in_prism_test_points_normal_sidewall.items[7] = make_vector3(70.24405, 114.905, 0.0000); // edge on bottom + point_in_prism_test_points_normal_sidewall.items[8] = make_vector3(41.1477, 0.00000, 100.000); // edge on side + point_in_prism_test_points_normal_sidewall.items[9] = make_vector3(140.488, 99.3401, 127.000); // vertex -> corner on top at edge of c + point_in_prism_test_points_normal_sidewall.items[10] = make_vector3(140.488, 99.3401, 129.000); // continuation of edge from vertex + point_in_prism_test_points_normal_sidewall.items[11] = make_vector3(141.902, 97.9259, 127.000); // continuation of edge from vertex + point_in_prism_test_points_normal_sidewall.items[12] = make_vector3(142.336, 100.105, 127.000); // continuation of edge from vertex + point_in_prism_test_points_normal_sidewall.items[13] = make_vector3(137.226, 99.9890, 125.000); // continuation of edge from vertex + point_in_prism_test_points_normal_sidewall.items[14] = make_vector3(25.5827, 88.7434, 127.000); // vertex -> corner on top inside c + point_in_prism_test_points_normal_sidewall.items[15] = make_vector3(25.5827, 88.7434, 129.000); // continuation of edge from vertex + point_in_prism_test_points_normal_sidewall.items[16] = make_vector3(24.1685, 87.3292, 127.000); // continuation of edge from vertex + point_in_prism_test_points_normal_sidewall.items[17] = make_vector3(25.5827, 90.7434, 127.000); // continuation of edge from vertex + point_in_prism_test_points_normal_sidewall.items[18] = make_vector3(26.9969, 88.1576, 125.000); // continuation of edge from vertex + point_in_prism_test_points_normal_sidewall.items[19] = make_vector3(41.1477, 0.00000, 127.000); // vertex -> corner on top outside c + point_in_prism_test_points_normal_sidewall.items[20] = make_vector3(114.905, 51.7447, 0.00000); // vertex -> corner on bottom at edge of c + point_in_prism_test_points_normal_sidewall.items[21] = make_vector3(51.7447, 114.905, 0.00000); // vertex -> corner on bottom inside c + point_in_prism_test_points_normal_sidewall.items[22] = make_vector3(0.00000, 99.3401, 0.00000); // vertex -> corner on bottom outside c + point_in_prism_test_points_normal_sidewall.items[23] = make_vector3(0.00000, 0.00000, 0.00000); // origin + point_in_prism_test_points_normal_sidewall.items[24] = make_vector3(70.2440, 70.2440, 63.5000); // center of the c + + int point_in_prism_expected_normal_sidewall[point_in_prism_test_points_normal_sidewall.num_items]; + point_in_prism_expected_normal_sidewall[0] = 1; // interior point + point_in_prism_expected_normal_sidewall[1] = 1; // interior point + point_in_prism_expected_normal_sidewall[2] = 1; // point on external side face + point_in_prism_expected_normal_sidewall[3] = 1; // point on internal side face + point_in_prism_expected_normal_sidewall[4] = 1; // point on top face + point_in_prism_expected_normal_sidewall[5] = 1; // point on bottom face + point_in_prism_expected_normal_sidewall[6] = 1; // edge on top + point_in_prism_expected_normal_sidewall[7] = 1; // edge on bottom + point_in_prism_expected_normal_sidewall[8] = 1; // edge on side + point_in_prism_expected_normal_sidewall[9] = 1; // vertex -> corner on top at edge of c + point_in_prism_expected_normal_sidewall[10] = 0; // continuation of edge from vertex + point_in_prism_expected_normal_sidewall[11] = 0; // continuation of edge from vertex + point_in_prism_expected_normal_sidewall[12] = 0; // continuation of edge from vertex + point_in_prism_expected_normal_sidewall[13] = 1; // continuation of edge from vertex + point_in_prism_expected_normal_sidewall[14] = 1; // vertex -> corner on top inside c + point_in_prism_expected_normal_sidewall[15] = 0; // continuation of edge from vertex + point_in_prism_expected_normal_sidewall[16] = 1; // continuation of edge from vertex + point_in_prism_expected_normal_sidewall[17] = 1; // continuation of edge from vertex + point_in_prism_expected_normal_sidewall[18] = 0; // continuation of edge from vertex + point_in_prism_expected_normal_sidewall[19] = 1; // vertex -> corner on top outside c + point_in_prism_expected_normal_sidewall[20] = 1; // vertex -> corner on bottom at edge of c + point_in_prism_expected_normal_sidewall[21] = 1; // vertex -> corner on bottom inside c + point_in_prism_expected_normal_sidewall[22] = 1; // vertex -> corner on bottom outside c + point_in_prism_expected_normal_sidewall[23] = 0; // origin + point_in_prism_expected_normal_sidewall[24] = 0; // center of the c + + int point_in_prism_actual_normal_sidewall[point_in_prism_test_points_normal_sidewall.num_items]; + for (int i = 0; i < point_in_prism_test_points_normal_sidewall.num_items; i++) { + num_tests_normal++; + point_in_prism_actual_normal_sidewall[i] = point_in_fixed_pobjectp(point_in_prism_test_points_normal_sidewall.items[i], &octagon_c_normal_sidewall_geom_object); + } + + for (int i = 0; i < point_in_prism_test_points_normal_sidewall.num_items; i++) { + if (point_in_prism_actual_normal_sidewall[i] != point_in_prism_expected_normal_sidewall[i]) { + ctl_printf("\tAt (%f, %f, %f) we expected point_in_fixed_pobjectp on the normal sidewall prism to return %i, but instead it returned %i\n", point_in_prism_test_points_normal_sidewall.items[i].x, point_in_prism_test_points_normal_sidewall.items[i].y, point_in_prism_test_points_normal_sidewall.items[i].z, point_in_prism_expected_normal_sidewall[i], point_in_prism_actual_normal_sidewall[i]); + num_failed_normal++; + } + } + + vector3_list point_in_prism_test_points_tapered_sidewall; + point_in_prism_test_points_tapered_sidewall.num_items = 25; + point_in_prism_test_points_tapered_sidewall.items = (vector3 *)malloc(point_in_prism_test_points_tapered_sidewall.num_items * sizeof(vector3)); + point_in_prism_test_points_tapered_sidewall.items[0] = make_vector3(46.446200000000005, 12.791350000000001, 63.500000000000000); // interior point + point_in_prism_test_points_tapered_sidewall.items[1] = make_vector3(123.45257948434455, 98.285670515655440, 63.500000000000000); // interior point + point_in_prism_test_points_tapered_sidewall.items[2] = make_vector3(102.72366312425248, 35.642282404302410, 63.500000000000000); // point on external side face + point_in_prism_test_points_tapered_sidewall.items[3] = make_vector3(21.423995187964220, 70.244056207821420, 95.250000000000000); // point on internal side face + point_in_prism_test_points_tapered_sidewall.items[4] = make_vector3(29.618783181268217, 110.86913318128589, 127.00000000000000); // point on top face + point_in_prism_test_points_tapered_sidewall.items[5] = make_vector3(134.09200000000000, 96.690900000000000, 0.0000000000000000); // point on bottom face + point_in_prism_test_points_tapered_sidewall.items[6] = make_vector3(20.037760250618970, 70.244062415642820, 127.00000000000000); // edge on top + point_in_prism_test_points_tapered_sidewall.items[7] = make_vector3(70.244050000000000, 114.90500000000000, 0.0000000000000000); // edge on bottom + point_in_prism_test_points_tapered_sidewall.items[8] = make_vector3(50.596305376632360, 22.810230125309488, 63.500000000000000); // edge on side + point_in_prism_test_points_tapered_sidewall.items[9] = make_vector3(130.69912049055401, 101.28725051332967, 127.00000000000000); // vertex -> corner on top at edge of c // failing + point_in_prism_test_points_tapered_sidewall.items[10] = make_vector3(130.62227962053330, 101.30253528002449, 127.99692620419043); // continuation of edge from vertex + point_in_prism_test_points_tapered_sidewall.items[11] = make_vector3(131.40622727174056, 100.58014373214311, 127.00000000000000); // continuation of edge from vertex + point_in_prism_test_points_tapered_sidewall.items[12] = make_vector3(131.62300162627434, 101.66993007518276, 127.00000000000000); // continuation of edge from vertex + point_in_prism_test_points_tapered_sidewall.items[13] = make_vector3(129.14497344366782, 101.59639296596830, 126.00307379580957); // continuation of edge from vertex + point_in_prism_test_points_tapered_sidewall.items[14] = make_vector3(20.037760250618966, 91.040214078020940, 127.00000000000000); // vertex -> corner on top inside c + point_in_prism_test_points_tapered_sidewall.items[15] = make_vector3(19.994147981293120, 91.058279066765540, 127.99888519167415); // continuation of edge from vertex + point_in_prism_test_points_tapered_sidewall.items[16] = make_vector3(19.330648063809882, 90.333112702498250, 127.00000000000000); // continuation of edge from vertex + point_in_prism_test_points_tapered_sidewall.items[17] = make_vector3(20.037760250618966, 92.040214078020940, 127.00000000000000); // continuation of edge from vertex + point_in_prism_test_points_tapered_sidewall.items[18] = make_vector3(20.788484706753895, 90.729250464799020, 126.00111480832585); // continuation of edge from vertex + point_in_prism_test_points_tapered_sidewall.items[19] = make_vector3(43.444489246735300, 5.5449397493810295, 127.00000000000000); // vertex -> corner on top outside c + point_in_prism_test_points_tapered_sidewall.items[20] = make_vector3(114.90500000000000, 51.744700000000000, 0.0000000000000000); // vertex -> corner on bottom at edge of c + point_in_prism_test_points_tapered_sidewall.items[21] = make_vector3(51.744700000000000, 114.90500000000000, 0.0000000000000000); // vertex -> corner on bottom inside c + point_in_prism_test_points_tapered_sidewall.items[22] = make_vector3(0.0000000000000000, 99.340100000000000, 0.0000000000000000); // vertex -> corner on bottom outside c + point_in_prism_test_points_tapered_sidewall.items[23] = make_vector3(0.0000000000000000, 0.0000000000000000, 0.0000000000000000); // origin + point_in_prism_test_points_tapered_sidewall.items[24] = make_vector3(70.244000000000000, 70.244000000000000, 63.500000000000000); // center of the c + + int point_in_prism_expected_tapered_sidewall[point_in_prism_test_points_tapered_sidewall.num_items]; + point_in_prism_expected_tapered_sidewall[0] = 1; // interior point + point_in_prism_expected_tapered_sidewall[1] = 1; // interior point + point_in_prism_expected_tapered_sidewall[2] = 1; // point on external side face + point_in_prism_expected_tapered_sidewall[3] = 1; // point on internal side face + point_in_prism_expected_tapered_sidewall[4] = 1; // point on top face + point_in_prism_expected_tapered_sidewall[5] = 1; // point on bottom face + point_in_prism_expected_tapered_sidewall[6] = 1; // edge on top + point_in_prism_expected_tapered_sidewall[7] = 1; // edge on bottom + point_in_prism_expected_tapered_sidewall[8] = 1; // edge on side + point_in_prism_expected_tapered_sidewall[9] = 1; // vertex -> corner on top at edge of c + point_in_prism_expected_tapered_sidewall[10] = 0; // continuation of edge from vertex + point_in_prism_expected_tapered_sidewall[11] = 0; // continuation of edge from vertex + point_in_prism_expected_tapered_sidewall[12] = 0; // continuation of edge from vertex + point_in_prism_expected_tapered_sidewall[13] = 1; // continuation of edge from vertex + point_in_prism_expected_tapered_sidewall[14] = 1; // vertex -> corner on top inside c + point_in_prism_expected_tapered_sidewall[15] = 0; // continuation of edge from vertex + point_in_prism_expected_tapered_sidewall[16] = 1; // continuation of edge from vertex + point_in_prism_expected_tapered_sidewall[17] = 1; // continuation of edge from vertex + point_in_prism_expected_tapered_sidewall[18] = 0; // continuation of edge from vertex + point_in_prism_expected_tapered_sidewall[19] = 1; // vertex -> corner on top outside c + point_in_prism_expected_tapered_sidewall[20] = 1; // vertex -> corner on bottom at edge of c + point_in_prism_expected_tapered_sidewall[21] = 1; // vertex -> corner on bottom inside c + point_in_prism_expected_tapered_sidewall[22] = 1; // vertex -> corner on bottom outside c + point_in_prism_expected_tapered_sidewall[23] = 0; // origin + point_in_prism_expected_tapered_sidewall[24] = 0; // center of the c + + int point_in_prism_actual_tapered_sidewall[point_in_prism_test_points_tapered_sidewall.num_items]; + for (int i = 0; i < point_in_prism_test_points_tapered_sidewall.num_items; i++) { + num_tests_tapered++; + point_in_prism_actual_tapered_sidewall[i] = point_in_fixed_pobjectp(point_in_prism_test_points_tapered_sidewall.items[i], &octagon_c_two_half_degree_sidewall_geom_object); + } + + for (int i = 0; i < point_in_prism_test_points_tapered_sidewall.num_items; i++) { + if (point_in_prism_actual_tapered_sidewall[i] != point_in_prism_expected_tapered_sidewall[i]) { + ctl_printf("\tAt (%f, %f, %f) we expected point_in_fixed_pobjectp on the tapered sidewall prism to return %i, but instead it returned %i\n", point_in_prism_test_points_tapered_sidewall.items[i].x, point_in_prism_test_points_tapered_sidewall.items[i].y, point_in_prism_test_points_tapered_sidewall.items[i].z, point_in_prism_expected_tapered_sidewall[i], point_in_prism_actual_tapered_sidewall[i]); + num_failed_tapered++; + } + } + + // test normal_to_prism + vector3_list normal_to_prism_test_points_normal_sidewall; + normal_to_prism_test_points_normal_sidewall.num_items = 30; + normal_to_prism_test_points_normal_sidewall.items = (vector3 *)malloc(normal_to_prism_test_points_normal_sidewall.num_items * sizeof(vector3)); + normal_to_prism_test_points_normal_sidewall.items[0] = make_vector3(98.2887, 98.2887, 63.5000); // points around sidewalls + normal_to_prism_test_points_normal_sidewall.items[1] = make_vector3(70.2441, 109.905, 63.5000); + normal_to_prism_test_points_normal_sidewall.items[2] = make_vector3(42.1992, 98.2886, 63.5000); + normal_to_prism_test_points_normal_sidewall.items[3] = make_vector3(30.5827, 70.2441, 63.5000); + normal_to_prism_test_points_normal_sidewall.items[4] = make_vector3(42.1992, 42.1992, 63.5000); + normal_to_prism_test_points_normal_sidewall.items[5] = make_vector3(70.2441, 30.5827, 63.5000); + normal_to_prism_test_points_normal_sidewall.items[6] = make_vector3(98.2886, 42.1992, 63.5000); + normal_to_prism_test_points_normal_sidewall.items[7] = make_vector3(129.610, 51.0656, 63.5000); + normal_to_prism_test_points_normal_sidewall.items[8] = make_vector3(123.450, 17.0383, 63.5000); + normal_to_prism_test_points_normal_sidewall.items[9] = make_vector3(70.2439, -5.0000, 63.5000); + normal_to_prism_test_points_normal_sidewall.items[10] = make_vector3(17.0383, 17.0383, 63.5000); + normal_to_prism_test_points_normal_sidewall.items[11] = make_vector3(-5.0000, 70.2439, 63.5000); + normal_to_prism_test_points_normal_sidewall.items[12] = make_vector3(17.0383, 123.450, 63.5000); + normal_to_prism_test_points_normal_sidewall.items[13] = make_vector3(70.2439, 145.488, 63.5000); + normal_to_prism_test_points_normal_sidewall.items[14] = make_vector3(123.450, 123.450, 63.5000); + normal_to_prism_test_points_normal_sidewall.items[15] = make_vector3(129.610, 89.4223, 63.5000); + normal_to_prism_test_points_normal_sidewall.items[16] = make_vector3(110.869, 110.869, -5.0000); // points off bottom face + normal_to_prism_test_points_normal_sidewall.items[17] = make_vector3(70.2440, 127.697, -5.0000); + normal_to_prism_test_points_normal_sidewall.items[18] = make_vector3(29.6188, 110.869, -5.0000); + normal_to_prism_test_points_normal_sidewall.items[19] = make_vector3(12.7914, 70.2440, -5.0000); + normal_to_prism_test_points_normal_sidewall.items[20] = make_vector3(29.6188, 29.6188, -5.0000); + normal_to_prism_test_points_normal_sidewall.items[21] = make_vector3(70.2440, 12.7914, -5.0000); + normal_to_prism_test_points_normal_sidewall.items[22] = make_vector3(110.869, 29.6188, -5.0000); + normal_to_prism_test_points_normal_sidewall.items[23] = make_vector3(108.747, 112.991, 132.000); // points off top face + normal_to_prism_test_points_normal_sidewall.items[24] = make_vector3(70.2440, 127.697, 132.000); + normal_to_prism_test_points_normal_sidewall.items[25] = make_vector3(29.6188, 110.869, 132.000); + normal_to_prism_test_points_normal_sidewall.items[26] = make_vector3(12.7914, 70.2440, 132.000); + normal_to_prism_test_points_normal_sidewall.items[27] = make_vector3(29.6188, 29.6188, 132.000); + normal_to_prism_test_points_normal_sidewall.items[28] = make_vector3(70.2440, 12.7914, 132.000); + normal_to_prism_test_points_normal_sidewall.items[29] = make_vector3(108.747, 27.4968, 132.000); + + vector3 normal_to_prism_expected_normal_sidewall[normal_to_prism_test_points_normal_sidewall.num_items]; + normal_to_prism_expected_normal_sidewall[0] = make_vector3(-0.707107, -0.707107, 0.0000000); // points around sidewalls + normal_to_prism_expected_normal_sidewall[1] = make_vector3(0.0000000, -1.000000, 0.0000000); + normal_to_prism_expected_normal_sidewall[2] = make_vector3(0.7071010, -0.707112, 0.0000000); + normal_to_prism_expected_normal_sidewall[3] = make_vector3(1.0000000, 0.0000000, 0.0000000); + normal_to_prism_expected_normal_sidewall[4] = make_vector3(0.7071070, 0.7071070, 0.0000000); + normal_to_prism_expected_normal_sidewall[5] = make_vector3(0.0000000, 1.0000000, 0.0000000); + normal_to_prism_expected_normal_sidewall[6] = make_vector3(-0.707112, 0.7071010, 0.0000000); + normal_to_prism_expected_normal_sidewall[7] = make_vector3(0.3826890, 0.9238770, 0.0000000); + normal_to_prism_expected_normal_sidewall[8] = make_vector3(0.7071050, -0.707108, 0.0000000); + normal_to_prism_expected_normal_sidewall[9] = make_vector3(0.0000000, -1.000000, 0.0000000); + normal_to_prism_expected_normal_sidewall[10] = make_vector3(-0.707107, -0.707107, 0.0000000); + normal_to_prism_expected_normal_sidewall[11] = make_vector3(-1.000000, 0.0000000, 0.0000000); + normal_to_prism_expected_normal_sidewall[12] = make_vector3(-0.707108, 0.7071050, 0.0000000); + normal_to_prism_expected_normal_sidewall[13] = make_vector3(0.0000000, 1.0000000, 0.0000000); + normal_to_prism_expected_normal_sidewall[14] = make_vector3(0.7071070, 0.7071070, 0.0000000); + normal_to_prism_expected_normal_sidewall[15] = make_vector3(0.3826800, -0.923881, 0.0000000); + normal_to_prism_expected_normal_sidewall[16] = make_vector3(0.0000000, 0.0000000, -1.000000); // points off bottom face + normal_to_prism_expected_normal_sidewall[17] = make_vector3(0.0000000, 0.0000000, -1.000000); + normal_to_prism_expected_normal_sidewall[18] = make_vector3(0.0000000, 0.0000000, -1.000000); + normal_to_prism_expected_normal_sidewall[19] = make_vector3(0.0000000, 0.0000000, -1.000000); + normal_to_prism_expected_normal_sidewall[20] = make_vector3(0.0000000, 0.0000000, -1.000000); + normal_to_prism_expected_normal_sidewall[21] = make_vector3(0.0000000, 0.0000000, -1.000000); + normal_to_prism_expected_normal_sidewall[22] = make_vector3(0.0000000, 0.0000000, -1.000000); + normal_to_prism_expected_normal_sidewall[23] = make_vector3(0.0000000, 0.0000000, 1.0000000); // points off top face + normal_to_prism_expected_normal_sidewall[24] = make_vector3(0.0000000, 0.0000000, 1.0000000); + normal_to_prism_expected_normal_sidewall[25] = make_vector3(0.0000000, 0.0000000, 1.0000000); + normal_to_prism_expected_normal_sidewall[26] = make_vector3(0.0000000, 0.0000000, 1.0000000); + normal_to_prism_expected_normal_sidewall[27] = make_vector3(0.0000000, 0.0000000, 1.0000000); + normal_to_prism_expected_normal_sidewall[28] = make_vector3(0.0000000, 0.0000000, 1.0000000); + normal_to_prism_expected_normal_sidewall[29] = make_vector3(0.0000000, 0.0000000, 1.0000000); + + vector3 normal_to_prism_actual_normal_sidewall[normal_to_prism_test_points_normal_sidewall.num_items]; + for (int i = 0; i < normal_to_prism_test_points_normal_sidewall.num_items; i++) { + num_tests_normal++; + normal_to_prism_actual_normal_sidewall[i] = unit_vector3(normal_to_object(normal_to_prism_test_points_normal_sidewall.items[i], octagon_c_normal_sidewall_geom_object)); + } + + for (int i = 0; i < normal_to_prism_test_points_normal_sidewall.num_items; i++) { + if (!vector3_nearly_equal(normal_to_prism_expected_normal_sidewall[i], normal_to_prism_actual_normal_sidewall[i], tolerance) + && !vector3_nearly_equal(normal_to_prism_expected_normal_sidewall[i], vector3_scale(-1, normal_to_prism_actual_normal_sidewall[i]), tolerance)) { + num_failed_normal++; + vector3 test_point = normal_to_prism_test_points_normal_sidewall.items[i]; + vector3 expected = normal_to_prism_expected_normal_sidewall[i]; + vector3 actual = normal_to_prism_actual_normal_sidewall[i]; + ctl_printf("\tAt (%f, %f, %f) the expected normal vector was (%f, %f, %f), but the actual\n\t\tnormal vector was (%f, %f, %f\n", test_point.x, test_point.y, test_point.z, expected.x, expected.y, expected.z, actual.x, actual.y, actual.z); + } + } + + vector3_list normal_to_prism_test_points_tapered_sidewall; + normal_to_prism_test_points_tapered_sidewall.num_items = 30; + normal_to_prism_test_points_tapered_sidewall.items = (vector3 *)malloc(normal_to_prism_test_points_tapered_sidewall.num_items * sizeof(vector3)); + normal_to_prism_test_points_tapered_sidewall.items[0] = make_vector3(106.256, 108.378, 63.2819); // points around sidewalls + normal_to_prism_test_points_tapered_sidewall.items[1] = make_vector3(70.2441, 122.673, 63.2819); + normal_to_prism_test_points_tapered_sidewall.items[2] = make_vector3(33.1711, 107.317, 63.2819); + normal_to_prism_test_points_tapered_sidewall.items[3] = make_vector3(17.8150, 70.2441, 63.2819); + normal_to_prism_test_points_tapered_sidewall.items[4] = make_vector3(33.1711, 33.1711, 63.2819); + normal_to_prism_test_points_tapered_sidewall.items[5] = make_vector3(70.2441, 17.8150, 63.2819); + normal_to_prism_test_points_tapered_sidewall.items[6] = make_vector3(106.256, 32.1101, 63.2819); + normal_to_prism_test_points_tapered_sidewall.items[7] = make_vector3(123.663, 39.7093, 63.2819); + normal_to_prism_test_points_tapered_sidewall.items[8] = make_vector3(113.360, 25.0055, 63.2819); + normal_to_prism_test_points_tapered_sidewall.items[9] = make_vector3(70.2439, 7.76771, 63.2819); + normal_to_prism_test_points_tapered_sidewall.items[10] = make_vector3(26.0665, 26.0665, 63.2819); + normal_to_prism_test_points_tapered_sidewall.items[11] = make_vector3(7.76771, 70.2439, 63.2819); + normal_to_prism_test_points_tapered_sidewall.items[12] = make_vector3(26.0665, 114.421, 63.2819); + normal_to_prism_test_points_tapered_sidewall.items[13] = make_vector3(70.2439, 132.720, 63.2819); + normal_to_prism_test_points_tapered_sidewall.items[14] = make_vector3(113.360, 115.482, 63.2819); + normal_to_prism_test_points_tapered_sidewall.items[15] = make_vector3(123.663, 100.779, 63.2819); + normal_to_prism_test_points_tapered_sidewall.items[16] = make_vector3(110.869, 110.869, -5.0000); // points off bottom face + normal_to_prism_test_points_tapered_sidewall.items[17] = make_vector3(70.2440, 127.697, -5.0000); + normal_to_prism_test_points_tapered_sidewall.items[18] = make_vector3(29.6188, 110.869, -5.0000); + normal_to_prism_test_points_tapered_sidewall.items[19] = make_vector3(12.7914, 70.2440, -5.0000); + normal_to_prism_test_points_tapered_sidewall.items[20] = make_vector3(29.6188, 29.6188, -5.0000); + normal_to_prism_test_points_tapered_sidewall.items[21] = make_vector3(70.2440, 12.7914, -5.0000); + normal_to_prism_test_points_tapered_sidewall.items[22] = make_vector3(110.869, 29.6188, -5.0000); + normal_to_prism_test_points_tapered_sidewall.items[23] = make_vector3(108.747, 112.991, 132.000); // points off top face + normal_to_prism_test_points_tapered_sidewall.items[24] = make_vector3(70.2440, 127.697, 132.000); + normal_to_prism_test_points_tapered_sidewall.items[25] = make_vector3(29.6188, 110.869, 132.000); + normal_to_prism_test_points_tapered_sidewall.items[26] = make_vector3(12.7914, 70.2440, 132.000); + normal_to_prism_test_points_tapered_sidewall.items[27] = make_vector3(29.6188, 29.6188, 132.000); + normal_to_prism_test_points_tapered_sidewall.items[28] = make_vector3(70.2440, 12.7914, 132.000); + normal_to_prism_test_points_tapered_sidewall.items[29] = make_vector3(108.747, 27.4968, 132.000); + + vector3 normal_to_prism_expected_tapered_sidewall[normal_to_prism_test_points_tapered_sidewall.num_items]; + normal_to_prism_expected_tapered_sidewall[0] = make_vector3(0.7064340, 0.7064340, -0.0436194); // points off top face + normal_to_prism_expected_tapered_sidewall[1] = make_vector3(0.0000000, 0.9990480, -0.0436194); + normal_to_prism_expected_tapered_sidewall[2] = make_vector3(-0.706428, 0.7064390, -0.0436194); + normal_to_prism_expected_tapered_sidewall[3] = make_vector3(-0.999048, 0.0000000, -0.0436194); + normal_to_prism_expected_tapered_sidewall[4] = make_vector3(-0.706434, -0.706434, -0.0436194); + normal_to_prism_expected_tapered_sidewall[5] = make_vector3(0.0000000, -0.999048, -0.0436194); + normal_to_prism_expected_tapered_sidewall[6] = make_vector3(0.7064390, -0.706428, -0.0436194); + normal_to_prism_expected_tapered_sidewall[7] = make_vector3(-0.382325, -0.922998, -0.0436194); + normal_to_prism_expected_tapered_sidewall[8] = make_vector3(-0.706432, 0.7064350, -0.0436194); + normal_to_prism_expected_tapered_sidewall[9] = make_vector3(0.0000000, 0.9990480, -0.0436194); + normal_to_prism_expected_tapered_sidewall[10] = make_vector3(0.7064340, 0.7064340, -0.0436194); + normal_to_prism_expected_tapered_sidewall[11] = make_vector3(0.9990480, 0.0000000, -0.0436194); + normal_to_prism_expected_tapered_sidewall[12] = make_vector3(0.7064350, -0.706432, -0.0436194); + normal_to_prism_expected_tapered_sidewall[13] = make_vector3(0.0000000, -0.999048, -0.0436194); + normal_to_prism_expected_tapered_sidewall[14] = make_vector3(-0.706434, -0.706434, -0.0436194); + normal_to_prism_expected_tapered_sidewall[15] = make_vector3(-0.382315, 0.9230020, -0.0436194); + normal_to_prism_expected_tapered_sidewall[16] = make_vector3(0.0000000, 0.0000000, -1.000000); // points off bottom face + normal_to_prism_expected_tapered_sidewall[17] = make_vector3(0.0000000, 0.0000000, -1.000000); + normal_to_prism_expected_tapered_sidewall[18] = make_vector3(0.0000000, 0.0000000, -1.000000); + normal_to_prism_expected_tapered_sidewall[19] = make_vector3(0.0000000, 0.0000000, -1.000000); + normal_to_prism_expected_tapered_sidewall[20] = make_vector3(0.0000000, 0.0000000, -1.000000); + normal_to_prism_expected_tapered_sidewall[21] = make_vector3(0.0000000, 0.0000000, -1.000000); + normal_to_prism_expected_tapered_sidewall[22] = make_vector3(0.0000000, 0.0000000, -1.000000); + normal_to_prism_expected_tapered_sidewall[23] = make_vector3(0.0000000, 0.0000000, 1.0000000); // points off top face + normal_to_prism_expected_tapered_sidewall[24] = make_vector3(0.0000000, 0.0000000, 1.0000000); + normal_to_prism_expected_tapered_sidewall[25] = make_vector3(0.0000000, 0.0000000, 1.0000000); + normal_to_prism_expected_tapered_sidewall[26] = make_vector3(0.0000000, 0.0000000, 1.0000000); + normal_to_prism_expected_tapered_sidewall[27] = make_vector3(0.0000000, 0.0000000, 1.0000000); + normal_to_prism_expected_tapered_sidewall[28] = make_vector3(0.0000000, 0.0000000, 1.0000000); + normal_to_prism_expected_tapered_sidewall[29] = make_vector3(0.0000000, 0.0000000, 1.0000000); + + vector3 normal_to_prism_actual_tapered_sidewall[normal_to_prism_test_points_tapered_sidewall.num_items]; + for (int i = 0; i < normal_to_prism_test_points_tapered_sidewall.num_items; i++) { + num_tests_tapered++; + normal_to_prism_actual_tapered_sidewall[i] = unit_vector3(normal_to_object(normal_to_prism_test_points_tapered_sidewall.items[i], octagon_c_two_half_degree_sidewall_geom_object)); + } + + for (int i = 0; i < normal_to_prism_test_points_tapered_sidewall.num_items; i++) { + if (!vector3_nearly_equal(normal_to_prism_expected_tapered_sidewall[i], normal_to_prism_actual_tapered_sidewall[i], tolerance) + && !vector3_nearly_equal(normal_to_prism_expected_tapered_sidewall[i], vector3_scale(-1, normal_to_prism_actual_tapered_sidewall[i]), tolerance)) { + num_failed_tapered++; + vector3 test_point = normal_to_prism_test_points_tapered_sidewall.items[i]; + vector3 expected = normal_to_prism_expected_tapered_sidewall[i]; + vector3 actual = normal_to_prism_actual_tapered_sidewall[i]; + ctl_printf("\tAt (%f, %f, %f) the expected normal vector was (%f, %f, %f), but the actual\n\t\tnormal vector was (%f, %f, %f\n", test_point.x, test_point.y, test_point.z, expected.x, expected.y, expected.z, actual.x, actual.y, actual.z); + } + } + + // test intersect_line_segment_with_prism + vector3_list intersect_line_with_prism_test_points_normal_sidewall; + intersect_line_with_prism_test_points_normal_sidewall.num_items = 9; + intersect_line_with_prism_test_points_normal_sidewall.items = (vector3 *)malloc(intersect_line_with_prism_test_points_normal_sidewall.num_items * sizeof(vector3)); + intersect_line_with_prism_test_points_normal_sidewall.items[0] = make_vector3(100.809, 144.033, 130.205); // line crossing top[15] to bottom[11] + intersect_line_with_prism_test_points_normal_sidewall.items[1] = make_vector3(17.0383, 123.450, 63.5000); // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_test_points_normal_sidewall.items[2] = make_vector3(17.0383, 123.450, 63.5000); // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_test_points_normal_sidewall.items[3] = make_vector3(17.0383, 123.450, 63.5000); // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_test_points_normal_sidewall.items[4] = make_vector3(12.7914, 70.2440, 63.5000); // interior point between 4, 5, 12, 13 in positive xhat + intersect_line_with_prism_test_points_normal_sidewall.items[5] = make_vector3(12.7914, 70.2440, 63.5000); // interior point between 4, 5, 12, 13 in negative xhat + intersect_line_with_prism_test_points_normal_sidewall.items[6] = make_vector3(40.1648, 142.861, 131.290); // between top point 14 and center of c on bottom + intersect_line_with_prism_test_points_normal_sidewall.items[7] = make_vector3(41.1477, 0.00000, 127.000); // between top point 11 and origin + intersect_line_with_prism_test_points_normal_sidewall.items[8] = make_vector3(51.7447, 114.905, 127.000); // between top point 3 and center of c on bottom + + vector3 intersect_line_with_prism_test_vectors_normal_sidewall[intersect_line_with_prism_test_points_normal_sidewall.num_items]; + intersect_line_with_prism_test_vectors_normal_sidewall[0] = make_vector3(-0.29372, -0.709099, -0.64102); // line crossing top[15] to bottom[11] + intersect_line_with_prism_test_vectors_normal_sidewall[1] = make_vector3(0.707107, -0.707107, 0.000000); // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_test_vectors_normal_sidewall[2] = make_vector3(0.707107, -0.707107, 0.000000); // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_test_vectors_normal_sidewall[3] = make_vector3(0.707107, -0.707107, 0.000000); // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_test_vectors_normal_sidewall[4] = make_vector3(1.000000, 0.0000000, 0.000000); // interior point between 4, 5, 12, 13 in positive xhat + intersect_line_with_prism_test_vectors_normal_sidewall[5] = make_vector3(-1.00000, 0.0000000, 0.000000); // interior point between 4, 5, 12, 13 in negative xhat + intersect_line_with_prism_test_vectors_normal_sidewall[6] = make_vector3(0.196571, -0.474559,-0.857994); // between top point 14 and center of c on bottom + intersect_line_with_prism_test_vectors_normal_sidewall[7] = make_vector3(-0.308223, 0.000000,-0.951314); // between top point 11 and origin + intersect_line_with_prism_test_vectors_normal_sidewall[8] = make_vector3(0.136135, -0.328658,-0.934586); // between top point 3 and center of c on bottom + + double intersect_line_with_prism_expected_normal_sidewall[intersect_line_with_prism_test_points_normal_sidewall.num_items]; + intersect_line_with_prism_expected_normal_sidewall[0] = 36.07816398; // line crossing top[15] to bottom[11] + intersect_line_with_prism_expected_normal_sidewall[1] = 25.58291121; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_expected_normal_sidewall[2] = 25.58291120; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_expected_normal_sidewall[3] = 51.16582241; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_expected_normal_sidewall[4] = 12.79135000; // interior point between 4, 5, 12, 13 in positive xhat + intersect_line_with_prism_expected_normal_sidewall[5] = 12.79135000; // interior point between 4, 5, 12, 13 in negative xhat + intersect_line_with_prism_expected_normal_sidewall[6] = 53.90914485; // between top point 14 and center of c on bottom + intersect_line_with_prism_expected_normal_sidewall[7] = 0.000000000; // between top point 11 and origin + intersect_line_with_prism_expected_normal_sidewall[8] = 0.000000000; // between top point 3 and center of c on bottom + + double intersect_line_with_prism_a_normal_sidewall[intersect_line_with_prism_test_points_normal_sidewall.num_items]; + intersect_line_with_prism_a_normal_sidewall[0] = 0; // line crossing top[15] to bottom[11] + intersect_line_with_prism_a_normal_sidewall[1] = 0; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_a_normal_sidewall[2] = 100; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_a_normal_sidewall[3] = 0; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_a_normal_sidewall[4] = 0; // interior point between 4, 5, 12, 13 in positive xhat + intersect_line_with_prism_a_normal_sidewall[5] = 0; // interior point between 4, 5, 12, 13 in negative xhat + intersect_line_with_prism_a_normal_sidewall[6] = 0; // between top point 14 and center of c on bottom + intersect_line_with_prism_a_normal_sidewall[7] = 0; // between top point 11 and origin + intersect_line_with_prism_a_normal_sidewall[8] = 0; // between top point 3 and center of c on bottom + + double intersect_line_with_prism_b_normal_sidewall[intersect_line_with_prism_test_points_normal_sidewall.num_items]; + intersect_line_with_prism_b_normal_sidewall[0] = 150; // line crossing top[15] to bottom[11] + intersect_line_with_prism_b_normal_sidewall[1] = 100; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_b_normal_sidewall[2] = 150; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_b_normal_sidewall[3] = 150; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_b_normal_sidewall[4] = 300; // interior point between 4, 5, 12, 13 in positive xhat + intersect_line_with_prism_b_normal_sidewall[5] = 300; // interior point between 4, 5, 12, 13 in negative xhat + intersect_line_with_prism_b_normal_sidewall[6] = 300; // between top point 14 and center of c on bottom + intersect_line_with_prism_b_normal_sidewall[7] = 300; // between top point 11 and origin + intersect_line_with_prism_b_normal_sidewall[8] = 300; // between top point 3 and center of c on bottom + + double intersect_line_with_prism_actual_normal_sidewall[intersect_line_with_prism_test_points_normal_sidewall.num_items]; + for (int i = 0; i < intersect_line_with_prism_test_points_normal_sidewall.num_items; i++) { + num_tests_normal++; + vector3 p = intersect_line_with_prism_test_points_normal_sidewall.items[i]; + vector3 d = intersect_line_with_prism_test_vectors_normal_sidewall[i]; + geometric_object o = octagon_c_normal_sidewall_geom_object; + double a = intersect_line_with_prism_a_normal_sidewall[i]; + double b = intersect_line_with_prism_b_normal_sidewall[i]; + intersect_line_with_prism_actual_normal_sidewall[i] = intersect_line_segment_with_object(p, d, o, a, b); + } + + for (int i = 0; i < intersect_line_with_prism_test_points_normal_sidewall.num_items; i++) { + double actual = intersect_line_with_prism_actual_normal_sidewall[i]; + double expected = intersect_line_with_prism_expected_normal_sidewall[i]; + if (fabs(fabs(actual)-fabs(expected)) > tolerance * fmax(fabs(actual), fabs(expected))) { + double px = intersect_line_with_prism_test_points_normal_sidewall.items[i].x; + double py = intersect_line_with_prism_test_points_normal_sidewall.items[i].y; + double pz = intersect_line_with_prism_test_points_normal_sidewall.items[i].z; + double dx = intersect_line_with_prism_test_vectors_normal_sidewall[i].x; + double dy = intersect_line_with_prism_test_vectors_normal_sidewall[i].y; + double dz = intersect_line_with_prism_test_vectors_normal_sidewall[i].z; + ctl_printf( + "\tThe line segment emanating from (%f, %f, %f) along s*d,\n\t\twith 0 <= s <= 300, d = (%f, %f, %f), was expected\n\t\tto have intersection length %f but instead had %f.\n", + px, py, pz, dx, dy, dz, expected, actual); + num_failed_normal++; + } + } + + vector3_list intersect_line_with_prism_test_points_tapered_sidewall; + intersect_line_with_prism_test_points_tapered_sidewall.num_items = 9; + intersect_line_with_prism_test_points_tapered_sidewall.items = (vector3 *)malloc(intersect_line_with_prism_test_points_tapered_sidewall.num_items * sizeof(vector3)); + intersect_line_with_prism_test_points_tapered_sidewall.items[0] = make_vector3(98.4872, 138.429, 130.281); // line crossing top[15] to bottom[11] + intersect_line_with_prism_test_points_tapered_sidewall.items[1] = make_vector3(19.0383, 121.528, 63.5000); // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_test_points_tapered_sidewall.items[2] = make_vector3(19.0383, 121.528, 63.5000); // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_test_points_tapered_sidewall.items[3] = make_vector3(19.0383, 121.528, 63.5000); // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_test_points_tapered_sidewall.items[4] = make_vector3(12.7914, 70.2440, 63.5000); // interior point between 4, 5, 12, 13 in positive xhat + intersect_line_with_prism_test_points_tapered_sidewall.items[5] = make_vector3(12.7914, 70.2440, 63.5000); // interior point between 4, 5, 12, 13 in negative xhat + intersect_line_with_prism_test_points_tapered_sidewall.items[6] = make_vector3(43.4445, 134.943, 127.000); // between top point 14 and center of c on bottom + intersect_line_with_prism_test_points_tapered_sidewall.items[7] = make_vector3(43.4445, 5.54494, 127.000); // between top point 11 and origin + intersect_line_with_prism_test_points_tapered_sidewall.items[8] = make_vector3(34.7428, 105.745, 127.000); // between top point 3 and center of c on bottom + + vector3 intersect_line_with_prism_test_vectors_tapered_sidewall[intersect_line_with_prism_test_points_tapered_sidewall.num_items]; + intersect_line_with_prism_test_vectors_tapered_sidewall[0] = make_vector3(-0.288786, -0.697187, -0.656149); // line crossing top[15] to bottom[11] + intersect_line_with_prism_test_vectors_tapered_sidewall[1] = make_vector3(0.6992010, -0.714925, 0.0000000); // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_test_vectors_tapered_sidewall[2] = make_vector3(0.6992010, -0.714925, 0.0000000); // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_test_vectors_tapered_sidewall[3] = make_vector3(0.6992010, -0.714925, 0.0000000); // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_test_vectors_tapered_sidewall[4] = make_vector3(1.0000000, 0.0000000, 0.0000000); // interior point between 4, 5, 12, 13 in positive xhat + intersect_line_with_prism_test_vectors_tapered_sidewall[5] = make_vector3(-1.000000, 0.0000000, 0.0000000); // interior point between 4, 5, 12, 13 in negative xhat + intersect_line_with_prism_test_vectors_tapered_sidewall[6] = make_vector3(0.1847880,-0.4461140, -0.875692); // between top point 14 and center of c on bottom + intersect_line_with_prism_test_vectors_tapered_sidewall[7] = make_vector3(-0.323393,-0.0412755, -0.945364); // between top point 11 and origin + intersect_line_with_prism_test_vectors_tapered_sidewall[8] = make_vector3(0.2599600,-0.2599600, -0.929969); // between top point 3 and center of c on bottom + + double intersect_line_with_prism_expected_tapered_sidewall[intersect_line_with_prism_test_points_tapered_sidewall.num_items]; + intersect_line_with_prism_expected_tapered_sidewall[0] = 21.67860775; // line crossing top[15] to bottom[11] + intersect_line_with_prism_expected_tapered_sidewall[1] = 20.03920840; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_expected_tapered_sidewall[2] = 20.03919670; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_expected_tapered_sidewall[3] = 40.07840510; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_expected_tapered_sidewall[4] = 10.01888013; // interior point between 4, 5, 12, 13 in positive xhat + intersect_line_with_prism_expected_tapered_sidewall[5] = 10.01888013; // interior point between 4, 5, 12, 13 in negative xhat + intersect_line_with_prism_expected_tapered_sidewall[6] = 35.53300107; // between top point 14 and center of c on bottom + intersect_line_with_prism_expected_tapered_sidewall[7] = 0.000000000; // between top point 11 and origin + intersect_line_with_prism_expected_tapered_sidewall[8] = 0.000000000; // between top point 3 and center of c on bottom + + double intersect_line_with_prism_a_tapered_sidewall[intersect_line_with_prism_test_points_tapered_sidewall.num_items]; + intersect_line_with_prism_a_tapered_sidewall[0] = 0; // line crossing top[15] to bottom[11] + intersect_line_with_prism_a_tapered_sidewall[1] = 0; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_a_tapered_sidewall[2] = 50; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_a_tapered_sidewall[3] = 0; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_a_tapered_sidewall[4] = 0; // interior point between 4, 5, 12, 13 in positive xhat + intersect_line_with_prism_a_tapered_sidewall[5] = 0; // interior point between 4, 5, 12, 13 in negative xhat + intersect_line_with_prism_a_tapered_sidewall[6] = 0; // between top point 14 and center of c on bottom + intersect_line_with_prism_a_tapered_sidewall[7] = 0; // between top point 11 and origin + intersect_line_with_prism_a_tapered_sidewall[8] = 0; // between top point 3 and center of c on bottom + + double intersect_line_with_prism_b_tapered_sidewall[intersect_line_with_prism_test_points_tapered_sidewall.num_items]; + intersect_line_with_prism_b_tapered_sidewall[0] = 150; // line crossing top[15] to bottom[11] + intersect_line_with_prism_b_tapered_sidewall[1] = 50; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_b_tapered_sidewall[2] = 150; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_b_tapered_sidewall[3] = 150; // line crossing midpoint 13-14 face to midpoint 9-10 face + intersect_line_with_prism_b_tapered_sidewall[4] = 300; // interior point between 4, 5, 12, 13 in positive xhat + intersect_line_with_prism_b_tapered_sidewall[5] = 300; // interior point between 4, 5, 12, 13 in negative xhat + intersect_line_with_prism_b_tapered_sidewall[6] = 300; // between top point 14 and center of c on bottom + intersect_line_with_prism_b_tapered_sidewall[7] = 300; // between top point 11 and origin + intersect_line_with_prism_b_tapered_sidewall[8] = 300; // between top point 3 and center of c on bottom + + double intersect_line_with_prism_actual_tapered_sidewall[intersect_line_with_prism_test_points_tapered_sidewall.num_items]; + for (int i = 0; i < intersect_line_with_prism_test_points_tapered_sidewall.num_items; i++) { + num_tests_tapered++; + vector3 p = intersect_line_with_prism_test_points_tapered_sidewall.items[i]; + vector3 d = intersect_line_with_prism_test_vectors_tapered_sidewall[i]; + geometric_object o = octagon_c_two_half_degree_sidewall_geom_object; + double a = intersect_line_with_prism_a_tapered_sidewall[i]; + double b = intersect_line_with_prism_b_tapered_sidewall[i]; + intersect_line_with_prism_actual_tapered_sidewall[i] = intersect_line_segment_with_object(p, d, o, a, b); + } + + for (int i = 0; i < intersect_line_with_prism_test_points_tapered_sidewall.num_items; i++) { + double actual = intersect_line_with_prism_actual_tapered_sidewall[i]; + double expected = intersect_line_with_prism_expected_tapered_sidewall[i]; + if (fabs(fabs(actual)-fabs(expected)) > tolerance * fmax(fabs(actual), fabs(expected))) { + double px = intersect_line_with_prism_test_points_tapered_sidewall.items[i].x; + double py = intersect_line_with_prism_test_points_tapered_sidewall.items[i].y; + double pz = intersect_line_with_prism_test_points_tapered_sidewall.items[i].z; + double dx = intersect_line_with_prism_test_vectors_tapered_sidewall[i].x; + double dy = intersect_line_with_prism_test_vectors_tapered_sidewall[i].y; + double dz = intersect_line_with_prism_test_vectors_tapered_sidewall[i].z; + ctl_printf( + "\tThe line segment emanating from (%f, %f, %f) along s*d,\n\t\twith 0 <= s <= 300, d = (%f, %f, %f), was expected\n\t\tto have intersection length %f but instead had %f.\n", + px, py, pz, dx, dy, dz, expected, actual); + num_failed_tapered++; + } + } + + printf("\tprism helper function testing summary: \n\t\t%i/%i tests failed with normal sidewall\n\t\t%i/%i tests failed with tapered sidewall\n", num_failed_normal, num_tests_normal, num_failed_tapered, num_tests_tapered); + + return num_failed_normal + num_failed_tapered; +} + /***************************************************************/ /* unit tests: create the same parallelepiped two ways (as a */ /* block and as a prism) and verify that geometric primitives */ @@ -506,11 +1141,14 @@ int run_unit_tests() { // 20180712 disabling this test because the new implementation of normal_to_object // for prisms is actually more accurate than the implementation for blocks, // although the distinction is only significant in cases where it is irrelevant - int num_failed_2 = 0; // test_normal_to_object(the_block, the_prism, NUMLINES, 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); int num_failed_4 = test_point_in_polygon(write_log); + int num_failed_5 = test_square_base_sidewall_prisms_to_gnuplot(); + int num_failed_6 = test_octagon_c_base_sidewall_prisms_to_gnuplot(); + int num_failed_7 = test_helper_functions_on_octagonal_c_prism(); - return num_failed_1 + num_failed_2 + num_failed_3 + num_failed_4; + return num_failed_1 + num_failed_3 + num_failed_4 + num_failed_5 + num_failed_6 + num_failed_7; } /***************************************************************/ @@ -606,7 +1244,7 @@ int main(int argc, char *argv[]) { /***************************************************************/ /* read vertices from vertex file and create prism *************/ /***************************************************************/ - vector3 *vertices = 0; + vector3 *vertices_bottom = 0; int num_vertices = 0; FILE *f = fopen(vertexfile, "r"); if (!f) usage("could not open vertexfile"); @@ -620,12 +1258,12 @@ int main(int argc, char *argv[]) { 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; + vertices_bottom = (vector3 *)realloc(vertices_bottom, num_vertices * sizeof(vector3)); + vertices_bottom[num_vertices - 1] = v; } fclose(f); - geometric_object the_prism = make_prism(NULL, vertices, num_vertices, height, axis); + geometric_object the_prism = make_prism(NULL, vertices_bottom, num_vertices, height, axis); prism *prsm = the_prism.subclass.prism_data; prism2gmsh(prsm, "test-prism.pp"); prism2gnuplot(prsm, "test-prism.gp");