From cb0244547dcf4382935019a15f1d897f346f7684 Mon Sep 17 00:00:00 2001 From: Luochenghuang Date: Tue, 7 Apr 2026 11:30:32 -0700 Subject: [PATCH 1/8] Add MESH geometric object type for arbitrary triangulated 3D surfaces Implement a new MESH primitive alongside SPHERE, CYLINDER, BLOCK, and PRISM. A mesh is defined by a vertex array and triangle index array, with a BVH (bounding volume hierarchy) for O(log N) queries. Core operations implemented: - point_in_mesh: ray casting with parity counting and deduplication - normal_to_mesh: BVH-accelerated closest-face search - get_mesh_bounding_box: from BVH root node - get_mesh_volume: divergence theorem (signed tetrahedron volumes) - intersect_line_segment_with_mesh: for subpixel smoothing integration - init_mesh: validates geometry, computes normals, fixes winding, builds BVH BVH uses SAH (surface area heuristic) binned construction with flat array layout. Ray-triangle intersection uses Moller-Trumbore algorithm. Closest-point queries use Eberly's algorithm with BVH pruning. Also adds make_mesh/make_mesh_with_center constructors, mesh copy/equal/destroy in geom-ctl-io.c, and 11 unit tests covering point-in, volume, bbox, normals, line-segment intersection, and comparison against the Block primitive. --- utils/Makefile.am | 8 +- utils/ctlgeom.h | 9 + utils/geom.c | 790 ++++++++++++++++++++++++++++++++++++++++++++++ utils/geom.scm | 10 + utils/test-mesh.c | 378 ++++++++++++++++++++++ 5 files changed, 1193 insertions(+), 2 deletions(-) create mode 100644 utils/test-mesh.c diff --git a/utils/Makefile.am b/utils/Makefile.am index d18e052..5f523bc 100644 --- a/utils/Makefile.am +++ b/utils/Makefile.am @@ -20,13 +20,17 @@ geomtst_SOURCES = geomtst.c geomtst_LDADD = libctlgeom.la geomtst_CPPFLAGS = -I$(top_srcdir)/src -check_PROGRAMS = test-prism +check_PROGRAMS = test-prism test-mesh test_prism_SOURCES = test-prism.c test_prism_LDADD = libctlgeom.la test_prism_CPPFLAGS = -I$(top_srcdir)/src -TESTS = test-prism +test_mesh_SOURCES = test-mesh.c +test_mesh_LDADD = libctlgeom.la +test_mesh_CPPFLAGS = -I$(top_srcdir)/src + +TESTS = test-prism test-mesh dist_man_MANS = gen-ctl-io.1 diff --git a/utils/ctlgeom.h b/utils/ctlgeom.h index 194b226..3ba8fbb 100644 --- a/utils/ctlgeom.h +++ b/utils/ctlgeom.h @@ -173,6 +173,15 @@ GEOMETRIC_OBJECT make_slanted_prism_with_center(MATERIAL_TYPE material, vector3 const vector3 *vertices, int num_vertices, double height, vector3 axis, double sidewall_angle); +// mesh with center computed as centroid of vertices +GEOMETRIC_OBJECT make_mesh(MATERIAL_TYPE material, const vector3 *vertices, int num_vertices, + const int *triangles, int num_triangles); + +// mesh with explicit center (vertices shifted so centroid = center) +GEOMETRIC_OBJECT make_mesh_with_center(MATERIAL_TYPE material, vector3 center, + const vector3 *vertices, int num_vertices, + const int *triangles, int num_triangles); + int vector3_nearly_equal(vector3 v1, vector3 v2, double tolerance); /**************************************************************************/ diff --git a/utils/geom.c b/utils/geom.c index dddf0c3..82d8eec 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -82,6 +82,14 @@ 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); static void reinit_prism(geometric_object *o); +static void init_mesh(geometric_object *o); +static boolean point_in_mesh(const mesh *m, vector3 p); +static vector3 normal_to_mesh(const mesh *m, vector3 p); +static void get_mesh_bounding_box(const mesh *m, geom_box *box); +static double get_mesh_volume(const mesh *m); +static void display_mesh_info(int indentby, const geometric_object *o); +static double intersect_line_segment_with_mesh(const mesh *m, vector3 p, vector3 d, + double a, double b); /**************************************************************************/ /* Allows writing to Python's stdout when running from Meep's Python interface */ @@ -174,6 +182,10 @@ void geom_fix_object_ptr(geometric_object *o) { } break; } + case GEOM MESH: { + init_mesh(o); + break; + } case GEOM GEOMETRIC_OBJECT_SELF: case GEOM SPHERE: break; /* these objects are fine */ } @@ -317,6 +329,9 @@ boolean point_in_fixed_pobjectp(vector3 p, geometric_object *o) { case GEOM PRISM: { return point_in_prism(o->subclass.prism_data, p); } + case GEOM MESH: { + return point_in_mesh(o->subclass.mesh_data, p); + } case GEOM COMPOUND_GEOMETRIC_OBJECT: { int i; int n = o->subclass.compound_geometric_object_data->component_objects.num_items; @@ -460,6 +475,8 @@ vector3 normal_to_fixed_object(vector3 p, geometric_object o) { case GEOM PRISM: return normal_to_prism(o.subclass.prism_data, p); + case GEOM MESH: return normal_to_mesh(o.subclass.mesh_data, p); + default: return r; } // switch (o.which_subclass) @@ -623,6 +640,7 @@ void CTLIO display_geometric_object_info(int indentby, geometric_object o) { } break; case GEOM PRISM: ctl_printf("prism"); break; + case GEOM MESH: ctl_printf("mesh"); break; case GEOM COMPOUND_GEOMETRIC_OBJECT: ctl_printf("compound object"); break; default: ctl_printf("geometric object"); break; } @@ -653,6 +671,7 @@ void CTLIO display_geometric_object_info(int indentby, geometric_object o) { o.subclass.block_data->e3.x, o.subclass.block_data->e3.y, o.subclass.block_data->e3.z); break; case GEOM PRISM: display_prism_info(indentby, &o); break; + case GEOM MESH: display_mesh_info(indentby, &o); break; case GEOM COMPOUND_GEOMETRIC_OBJECT: { int i; int n = o.subclass.compound_geometric_object_data->component_objects.num_items; @@ -845,6 +864,9 @@ double intersect_line_segment_with_object(vector3 p, vector3 d, geometric_object if (o.which_subclass == GEOM PRISM) { return intersect_line_segment_with_prism(o.subclass.prism_data, p, d, a, b); } + else if (o.which_subclass == GEOM MESH) { + return intersect_line_segment_with_mesh(o.subclass.mesh_data, p, d, a, b); + } else { double s[2]; if (2 == intersect_line_with_object(p, d, o, s)) { @@ -917,6 +939,9 @@ double geom_object_volume(GEOMETRIC_OBJECT o) { case GEOM PRISM: { return get_prism_volume(o.subclass.prism_data); } + case GEOM MESH: { + return get_mesh_volume(o.subclass.mesh_data); + } default: return 0; /* unsupported object types? */ } } @@ -1119,6 +1144,10 @@ void geom_get_bounding_box(geometric_object o, geom_box *box) { get_prism_bounding_box(o.subclass.prism_data, box); break; } + case GEOM MESH: { + get_mesh_bounding_box(o.subclass.mesh_data, box); + break; + } case GEOM COMPOUND_GEOMETRIC_OBJECT: { int i; int n = o.subclass.compound_geometric_object_data->component_objects.num_items; @@ -1894,6 +1923,767 @@ geometric_object make_ellipsoid(material_type material, vector3 center, vector3 return o; } +/*************************************************************** + * Mesh geometry implementation. + * A mesh is a closed triangulated surface defined by a vertex + * array and a triangle index array. A BVH (bounding volume + * hierarchy) is used for O(log N) queries. + ***************************************************************/ + +/***************************************************************/ +/* Triangle utilities */ +/***************************************************************/ + +/* Helper to get a geom_box from a BVH node. */ +static geom_box bvh_node_box(const mesh_bvh_node *node) { + geom_box b; + b.low = node->bbox_low; + b.high = node->bbox_high; + return b; +} + +/* Helper to set a BVH node's bbox from a geom_box. */ +static void bvh_node_set_box(mesh_bvh_node *node, const geom_box *box) { + node->bbox_low = box->low; + node->bbox_high = box->high; +} + +/* Compute the AABB of a single triangle. */ +static void mesh_triangle_bbox(const mesh *m, int face_id, geom_box *box) { + int i0 = m->face_indices[3 * face_id]; + int i1 = m->face_indices[3 * face_id + 1]; + int i2 = m->face_indices[3 * face_id + 2]; + vector3 v0 = m->vertices.items[i0]; + vector3 v1 = m->vertices.items[i1]; + vector3 v2 = m->vertices.items[i2]; + box->low.x = fmin(v0.x, fmin(v1.x, v2.x)); + box->low.y = fmin(v0.y, fmin(v1.y, v2.y)); + box->low.z = fmin(v0.z, fmin(v1.z, v2.z)); + box->high.x = fmax(v0.x, fmax(v1.x, v2.x)); + box->high.y = fmax(v0.y, fmax(v1.y, v2.y)); + box->high.z = fmax(v0.z, fmax(v1.z, v2.z)); +} + +/* Compute the centroid of a triangle. */ +static vector3 mesh_triangle_centroid(const mesh *m, int face_id) { + int i0 = m->face_indices[3 * face_id]; + int i1 = m->face_indices[3 * face_id + 1]; + int i2 = m->face_indices[3 * face_id + 2]; + vector3 v0 = m->vertices.items[i0]; + vector3 v1 = m->vertices.items[i1]; + vector3 v2 = m->vertices.items[i2]; + vector3 c; + c.x = (v0.x + v1.x + v2.x) / 3.0; + c.y = (v0.y + v1.y + v2.y) / 3.0; + c.z = (v0.z + v1.z + v2.z) / 3.0; + return c; +} + +/* Compute the surface area of a geom_box. */ +static double geom_box_surface_area(const geom_box *b) { + double dx = b->high.x - b->low.x; + double dy = b->high.y - b->low.y; + double dz = b->high.z - b->low.z; + return 2.0 * (dx * dy + dy * dz + dz * dx); +} + +/***************************************************************/ +/* BVH construction using Surface Area Heuristic (SAH) */ +/***************************************************************/ + +#define MESH_BVH_MAX_LEAF_SIZE 4 +#define MESH_BVH_NUM_BINS 12 + +/* Recursive BVH build. Returns the index of the root node for this subtree. + face_ids[start..start+count-1] are the faces in this node. */ +static int mesh_bvh_build(mesh *m, int *face_ids, int start, int count, + mesh_bvh_node *nodes, int *num_nodes) { + int node_idx = (*num_nodes)++; + mesh_bvh_node *node = &nodes[node_idx]; + + /* Compute bounding box of all faces in this range. */ + geom_box face_box, node_box; + mesh_triangle_bbox(m, face_ids[start], &face_box); + node_box = face_box; + for (int i = 1; i < count; i++) { + mesh_triangle_bbox(m, face_ids[start + i], &face_box); + geom_box_union(&node_box, &node_box, &face_box); + } + bvh_node_set_box(node, &node_box); + + /* Leaf node if few enough faces. */ + if (count <= MESH_BVH_MAX_LEAF_SIZE) { + node->left_child = -1; + node->right_child = -1; + node->face_start = start; + node->face_count = count; + return node_idx; + } + + /* Find the best split using SAH with binning. */ + double best_cost = 1e300; + int best_axis = -1, best_split = -1; + double parent_area = geom_box_surface_area(&node_box); + if (parent_area == 0) parent_area = 1e-30; + + for (int axis = 0; axis < 3; axis++) { + double lo, hi; + switch (axis) { + case 0: lo = node_box.low.x; hi = node_box.high.x; break; + case 1: lo = node_box.low.y; hi = node_box.high.y; break; + default: lo = node_box.low.z; hi = node_box.high.z; break; + } + if (hi - lo < 1e-15) continue; + + /* Bin face centroids. */ + int bin_counts[MESH_BVH_NUM_BINS]; + geom_box bin_boxes[MESH_BVH_NUM_BINS]; + for (int b = 0; b < MESH_BVH_NUM_BINS; b++) { + bin_counts[b] = 0; + bin_boxes[b].low.x = bin_boxes[b].low.y = bin_boxes[b].low.z = 1e300; + bin_boxes[b].high.x = bin_boxes[b].high.y = bin_boxes[b].high.z = -1e300; + } + + double inv_range = MESH_BVH_NUM_BINS / (hi - lo); + for (int i = 0; i < count; i++) { + vector3 c = mesh_triangle_centroid(m, face_ids[start + i]); + double val; + switch (axis) { + case 0: val = c.x; break; + case 1: val = c.y; break; + default: val = c.z; break; + } + int bin = (int)((val - lo) * inv_range); + if (bin < 0) bin = 0; + if (bin >= MESH_BVH_NUM_BINS) bin = MESH_BVH_NUM_BINS - 1; + bin_counts[bin]++; + geom_box fb; + mesh_triangle_bbox(m, face_ids[start + i], &fb); + geom_box_union(&bin_boxes[bin], &bin_boxes[bin], &fb); + } + + /* Sweep from left to evaluate SAH cost at each split. */ + geom_box left_box; + left_box.low.x = left_box.low.y = left_box.low.z = 1e300; + left_box.high.x = left_box.high.y = left_box.high.z = -1e300; + int left_count = 0; + + /* Precompute right sweep. */ + geom_box right_boxes[MESH_BVH_NUM_BINS]; + int right_counts[MESH_BVH_NUM_BINS]; + { + geom_box rb; + rb.low.x = rb.low.y = rb.low.z = 1e300; + rb.high.x = rb.high.y = rb.high.z = -1e300; + int rc = 0; + for (int b = MESH_BVH_NUM_BINS - 1; b >= 0; b--) { + if (bin_counts[b] > 0) { + geom_box_union(&rb, &rb, &bin_boxes[b]); + rc += bin_counts[b]; + } + right_boxes[b] = rb; + right_counts[b] = rc; + } + } + + for (int split = 1; split < MESH_BVH_NUM_BINS; split++) { + int b = split - 1; + if (bin_counts[b] > 0) { + geom_box_union(&left_box, &left_box, &bin_boxes[b]); + left_count += bin_counts[b]; + } + if (left_count == 0 || right_counts[split] == 0) continue; + + double left_area = geom_box_surface_area(&left_box); + double right_area = geom_box_surface_area(&right_boxes[split]); + double cost = 1.0 + (left_area * left_count + right_area * right_counts[split]) / parent_area; + if (cost < best_cost) { + best_cost = cost; + best_axis = axis; + best_split = split; + } + } + } + + /* If no good split found, make a leaf. */ + if (best_axis < 0) { + node->left_child = -1; + node->right_child = -1; + node->face_start = start; + node->face_count = count; + return node_idx; + } + + /* Partition face_ids by the best split. */ + double lo, hi; + switch (best_axis) { + case 0: lo = node_box.low.x; hi = node_box.high.x; break; + case 1: lo = node_box.low.y; hi = node_box.high.y; break; + default: lo = node_box.low.z; hi = node_box.high.z; break; + } + double inv_range = MESH_BVH_NUM_BINS / (hi - lo); + int left_end = start; + for (int i = start; i < start + count; i++) { + vector3 c = mesh_triangle_centroid(m, face_ids[i]); + double val; + switch (best_axis) { + case 0: val = c.x; break; + case 1: val = c.y; break; + default: val = c.z; break; + } + int bin = (int)((val - lo) * inv_range); + if (bin < 0) bin = 0; + if (bin >= MESH_BVH_NUM_BINS) bin = MESH_BVH_NUM_BINS - 1; + if (bin < best_split) { + int tmp = face_ids[left_end]; + face_ids[left_end] = face_ids[i]; + face_ids[i] = tmp; + left_end++; + } + } + + /* Handle degenerate partition. */ + int left_count_final = left_end - start; + if (left_count_final == 0 || left_count_final == count) { + left_count_final = count / 2; + left_end = start + left_count_final; + } + + node->face_start = -1; + node->face_count = 0; + node->left_child = mesh_bvh_build(m, face_ids, start, left_count_final, nodes, num_nodes); + /* Re-fetch node pointer since array may have been indexed differently. */ + node = &nodes[node_idx]; + node->right_child = mesh_bvh_build(m, face_ids, left_end, count - left_count_final, nodes, num_nodes); + return node_idx; +} + +/***************************************************************/ +/* Ray-triangle intersection (Moller-Trumbore) */ +/***************************************************************/ + +#define MESH_RAY_EPS 1e-12 + +/* Test intersection of ray (origin + t * dir) with triangle (v0, v1, v2). + Returns 1 if intersection found, with *t_out set to the parameter. + Also returns barycentric coords u, v if non-NULL. */ +static int ray_triangle_intersect(vector3 origin, vector3 dir, + vector3 v0, vector3 v1, vector3 v2, + double *t_out, double *u_out, double *v_out) { + vector3 e1 = vector3_minus(v1, v0); + vector3 e2 = vector3_minus(v2, v0); + vector3 h = vector3_cross(dir, e2); + double a = vector3_dot(e1, h); + + if (fabs(a) < MESH_RAY_EPS) return 0; + + double f = 1.0 / a; + vector3 s = vector3_minus(origin, v0); + double u = f * vector3_dot(s, h); + if (u < -MESH_RAY_EPS || u > 1.0 + MESH_RAY_EPS) return 0; + + vector3 q = vector3_cross(s, e1); + double v = f * vector3_dot(dir, q); + if (v < -MESH_RAY_EPS || u + v > 1.0 + MESH_RAY_EPS) return 0; + + double t = f * vector3_dot(e2, q); + *t_out = t; + if (u_out) *u_out = u; + if (v_out) *v_out = v; + return 1; +} + +/***************************************************************/ +/* Ray-AABB intersection test */ +/***************************************************************/ + +static int ray_bvh_node_intersect(vector3 origin, vector3 inv_dir, const mesh_bvh_node *node, + double t_min, double t_max) { + double tx1 = (node->bbox_low.x - origin.x) * inv_dir.x; + double tx2 = (node->bbox_high.x - origin.x) * inv_dir.x; + double tmin = fmin(tx1, tx2); + double tmax = fmax(tx1, tx2); + + double ty1 = (node->bbox_low.y - origin.y) * inv_dir.y; + double ty2 = (node->bbox_high.y - origin.y) * inv_dir.y; + tmin = fmax(tmin, fmin(ty1, ty2)); + tmax = fmin(tmax, fmax(ty1, ty2)); + + double tz1 = (node->bbox_low.z - origin.z) * inv_dir.z; + double tz2 = (node->bbox_high.z - origin.z) * inv_dir.z; + tmin = fmax(tmin, fmin(tz1, tz2)); + tmax = fmin(tmax, fmax(tz1, tz2)); + + return tmax >= fmax(tmin, t_min) && tmin <= t_max; +} + +/***************************************************************/ +/* Closest point on triangle */ +/***************************************************************/ + +/* Compute the closest point on triangle (v0,v1,v2) to point p. + Returns the squared distance. */ +static double closest_point_on_triangle(vector3 p, vector3 v0, vector3 v1, vector3 v2, + vector3 *closest) { + vector3 e0 = vector3_minus(v1, v0); + vector3 e1 = vector3_minus(v2, v0); + vector3 v = vector3_minus(v0, p); + + double a = vector3_dot(e0, e0); + double b = vector3_dot(e0, e1); + double c = vector3_dot(e1, e1); + double d = vector3_dot(e0, v); + double e = vector3_dot(e1, v); + + double det = a * c - b * b; + double s = b * e - c * d; + double t = b * d - a * e; + + if (s + t <= det) { + if (s < 0) { + if (t < 0) { + /* region 4 */ + if (d < 0) { + t = 0; + s = (-d >= a) ? 1 : -d / a; + } else { + s = 0; + t = (e >= 0) ? 0 : ((-e >= c) ? 1 : -e / c); + } + } else { + /* region 3 */ + s = 0; + t = (e >= 0) ? 0 : ((-e >= c) ? 1 : -e / c); + } + } else if (t < 0) { + /* region 5 */ + t = 0; + s = (d >= 0) ? 0 : ((-d >= a) ? 1 : -d / a); + } else { + /* region 0 (interior) */ + double inv_det = 1.0 / det; + s *= inv_det; + t *= inv_det; + } + } else { + if (s < 0) { + /* region 2 */ + double tmp0 = b + d; + double tmp1 = c + e; + if (tmp1 > tmp0) { + double numer = tmp1 - tmp0; + double denom = a - 2 * b + c; + s = (numer >= denom) ? 1 : numer / denom; + t = 1 - s; + } else { + s = 0; + t = (tmp1 <= 0) ? 1 : ((e >= 0) ? 0 : -e / c); + } + } else if (t < 0) { + /* region 6 */ + double tmp0 = b + e; + double tmp1 = a + d; + if (tmp1 > tmp0) { + double numer = tmp1 - tmp0; + double denom = a - 2 * b + c; + t = (numer >= denom) ? 1 : numer / denom; + s = 1 - t; + } else { + t = 0; + s = (tmp1 <= 0) ? 1 : ((d >= 0) ? 0 : -d / a); + } + } else { + /* region 1 */ + double numer = (c + e) - (b + d); + if (numer <= 0) { + s = 0; + t = 1; + } else { + double denom = a - 2 * b + c; + s = (numer >= denom) ? 1 : numer / denom; + t = 1 - s; + } + } + } + + *closest = vector3_plus(v0, vector3_plus(vector3_scale(s, e0), vector3_scale(t, e1))); + vector3 diff = vector3_minus(p, *closest); + return vector3_dot(diff, diff); +} + +/***************************************************************/ +/* BVH-accelerated mesh queries */ +/***************************************************************/ + +/* Find the closest face to point p using BVH traversal. + Returns the face index and sets *dist2 to the squared distance. */ +static int find_closest_face(const mesh *m, vector3 p, double *dist2) { + int best_face = -1; + double best_dist2 = 1e300; + + /* Stack-based traversal with pruning. */ + int stack[64]; + int stack_top = 0; + stack[stack_top++] = 0; + + while (stack_top > 0) { + int node_idx = stack[--stack_top]; + const mesh_bvh_node *node = &m->bvh[node_idx]; + + /* Compute minimum squared distance from p to the AABB. */ + double dx = fmax(0, fmax(node->bbox_low.x - p.x, p.x - node->bbox_high.x)); + double dy = fmax(0, fmax(node->bbox_low.y - p.y, p.y - node->bbox_high.y)); + double dz = fmax(0, fmax(node->bbox_low.z - p.z, p.z - node->bbox_high.z)); + double box_dist2 = dx * dx + dy * dy + dz * dz; + if (box_dist2 >= best_dist2) continue; + + if (node->left_child < 0) { + /* Leaf: test all faces. */ + for (int i = 0; i < node->face_count; i++) { + int fid = m->bvh_face_ids[node->face_start + i]; + vector3 v0 = m->vertices.items[m->face_indices[3 * fid]]; + vector3 v1 = m->vertices.items[m->face_indices[3 * fid + 1]]; + vector3 v2 = m->vertices.items[m->face_indices[3 * fid + 2]]; + vector3 closest; + double d2 = closest_point_on_triangle(p, v0, v1, v2, &closest); + if (d2 < best_dist2) { + best_dist2 = d2; + best_face = fid; + } + } + } else { + if (stack_top + 2 > 64) continue; + stack[stack_top++] = node->left_child; + stack[stack_top++] = node->right_child; + } + } + + *dist2 = best_dist2; + return best_face; +} + +/* Collect all ray-mesh intersection parameters into slist. + Returns the number of intersections. */ +static int mesh_ray_all_intersections(const mesh *m, vector3 origin, vector3 dir, + double *slist, int max_intersections) { + vector3 inv_dir; + inv_dir.x = (fabs(dir.x) > 1e-30) ? 1.0 / dir.x : 1e30; + inv_dir.y = (fabs(dir.y) > 1e-30) ? 1.0 / dir.y : 1e30; + inv_dir.z = (fabs(dir.z) > 1e-30) ? 1.0 / dir.z : 1e30; + + int count = 0; + int stack[64]; + int stack_top = 0; + stack[stack_top++] = 0; + + while (stack_top > 0) { + int node_idx = stack[--stack_top]; + const mesh_bvh_node *node = &m->bvh[node_idx]; + + if (!ray_bvh_node_intersect(origin, inv_dir, node, -1e30, 1e30)) + continue; + + if (node->left_child < 0) { + for (int i = 0; i < node->face_count; i++) { + int fid = m->bvh_face_ids[node->face_start + i]; + vector3 v0 = m->vertices.items[m->face_indices[3 * fid]]; + vector3 v1 = m->vertices.items[m->face_indices[3 * fid + 1]]; + vector3 v2 = m->vertices.items[m->face_indices[3 * fid + 2]]; + + double t; + if (ray_triangle_intersect(origin, dir, v0, v1, v2, &t, NULL, NULL)) { + if (count < max_intersections) + slist[count++] = t; + } + } + } else { + if (stack_top + 2 > 64) continue; + stack[stack_top++] = node->left_child; + stack[stack_top++] = node->right_child; + } + } + + return count; +} + +static int mesh_dcmp(const void *a, const void *b) { + double da = *(const double *)a; + double db = *(const double *)b; + return (da > db) - (da < db); +} + +/* Remove duplicate intersection parameters within tolerance. */ +static int remove_duplicate_intersections(double *slist, int count, double tol) { + if (count <= 1) return count; + int j = 0; + for (int i = 1; i < count; i++) { + if (slist[i] - slist[j] > tol) + slist[++j] = slist[i]; + } + return j + 1; +} + +#define MESH_MAX_INTERSECTIONS 256 + +/* Count ray-mesh intersections using BVH traversal. + Collects all intersection t-values, deduplicates, and returns the count. */ +static int count_ray_mesh_intersections(const mesh *m, vector3 origin, vector3 dir) { + double slist[MESH_MAX_INTERSECTIONS]; + int count = mesh_ray_all_intersections(m, origin, dir, slist, MESH_MAX_INTERSECTIONS); + + /* Keep only forward intersections (t > eps). */ + int nforward = 0; + for (int i = 0; i < count; i++) + if (slist[i] > MESH_RAY_EPS) + slist[nforward++] = slist[i]; + + if (nforward > 1) { + qsort(slist, nforward, sizeof(double), mesh_dcmp); + nforward = remove_duplicate_intersections(slist, nforward, 1e-10); + } + + return nforward; +} + +/***************************************************************/ +/* Core mesh geometric operations */ +/***************************************************************/ + +static boolean point_in_mesh(const mesh *m, vector3 p) { + if (!m->is_closed) return 0; + + /* Cast a ray along a non-axis-aligned direction to reduce chance of + hitting mesh edges exactly (where deduplication is most fragile). */ + vector3 dir = {0.57735026918962576, 0.57735026918962576, 0.57735026918962576}; + int crossings = count_ray_mesh_intersections(m, p, dir); + return (crossings % 2) == 1; +} + +static vector3 normal_to_mesh(const mesh *m, vector3 p) { + double dist2; + int face = find_closest_face(m, p, &dist2); + if (face < 0) { + vector3 zero = {0, 0, 0}; + return zero; + } + return m->face_normals[face]; +} + +static void get_mesh_bounding_box(const mesh *m, geom_box *box) { + if (m->num_bvh_nodes > 0) { + box->low = m->bvh[0].bbox_low; + box->high = m->bvh[0].bbox_high; + } else { + box->low = box->high = m->vertices.items[0]; + for (int i = 1; i < m->vertices.num_items; i++) + geom_box_add_pt(box, m->vertices.items[i]); + } +} + +static double get_mesh_volume(const mesh *m) { + /* Divergence theorem: sum signed tetrahedron volumes. */ + double vol = 0; + for (int f = 0; f < m->num_faces; f++) { + vector3 v0 = m->vertices.items[m->face_indices[3 * f]]; + vector3 v1 = m->vertices.items[m->face_indices[3 * f + 1]]; + vector3 v2 = m->vertices.items[m->face_indices[3 * f + 2]]; + vol += vector3_dot(v0, vector3_cross(v1, v2)); + } + return fabs(vol) / 6.0; +} + +static void display_mesh_info(int indentby, const geometric_object *o) { + const mesh *m = o->subclass.mesh_data; + ctl_printf("%*s %d vertices, %d faces, %s\n", indentby, "", + m->vertices.num_items, m->num_faces, + m->is_closed ? "closed" : "OPEN (WARNING)"); +} + +static int intersect_line_with_mesh(const mesh *m, vector3 p, vector3 d, double *slist) { + int count = mesh_ray_all_intersections(m, p, d, slist, MESH_MAX_INTERSECTIONS); + if (count > 1) { + qsort(slist, count, sizeof(double), mesh_dcmp); + count = remove_duplicate_intersections(slist, count, 1e-10); + } + return count; +} + +static double intersect_line_segment_with_mesh(const mesh *m, vector3 p, vector3 d, + double a, double b) { + double slist[MESH_MAX_INTERSECTIONS]; + int num = intersect_line_with_mesh(m, p, d, slist); + + /* Find first intersection > a. */ + int na = -1; + for (int i = 0; na == -1 && i < num; i++) + if (slist[i] > a) na = i; + if (na == -1) { + /* No intersections after a. Check if p+a*d is inside. */ + vector3 pt = vector3_plus(p, vector3_scale(a, d)); + if (point_in_mesh(m, pt)) + return b - a; + return 0.0; + } + + /* Determine if we start inside (even-odd parity). */ + int inside = ((na % 2) == 0) ? 0 : 1; + double last_s = a, ds = 0.0; + for (int i = na; i < num; i++) { + double this_s = fmin(b, slist[i]); + if (inside) ds += (this_s - last_s); + if (b <= slist[i]) break; + inside = 1 - inside; + last_s = this_s; + } + if (inside && last_s < b) ds += (b - last_s); + return ds > 0.0 ? ds : 0.0; +} + +/***************************************************************/ +/* init_mesh: validate, compute normals, build BVH */ +/***************************************************************/ + +static void init_mesh(geometric_object *o) { + mesh *m = o->subclass.mesh_data; + int nv = m->vertices.num_items; + int nf = m->num_faces; + + /* Validate. */ + CHECK(nv >= 4, "mesh requires at least 4 vertices"); + CHECK(nf >= 4, "mesh requires at least 4 faces"); + for (int f = 0; f < nf; f++) { + CHECK(m->face_indices[3 * f] >= 0 && m->face_indices[3 * f] < nv, "mesh face index out of range"); + CHECK(m->face_indices[3 * f + 1] >= 0 && m->face_indices[3 * f + 1] < nv, "mesh face index out of range"); + CHECK(m->face_indices[3 * f + 2] >= 0 && m->face_indices[3 * f + 2] < nv, "mesh face index out of range"); + } + + /* Compute face normals and areas. */ + m->face_normals = (vector3 *)malloc(nf * sizeof(vector3)); + CHECK(m->face_normals, "out of memory"); + m->face_areas = (double *)malloc(nf * sizeof(double)); + CHECK(m->face_areas, "out of memory"); + + for (int f = 0; f < nf; f++) { + vector3 v0 = m->vertices.items[m->face_indices[3 * f]]; + vector3 v1 = m->vertices.items[m->face_indices[3 * f + 1]]; + vector3 v2 = m->vertices.items[m->face_indices[3 * f + 2]]; + vector3 e1 = vector3_minus(v1, v0); + vector3 e2 = vector3_minus(v2, v0); + vector3 n = vector3_cross(e1, e2); + double len = vector3_norm(n); + m->face_areas[f] = 0.5 * len; + m->face_normals[f] = (len > 1e-30) ? vector3_scale(1.0 / len, n) : n; + } + + /* Compute centroid. */ + m->centroid.x = m->centroid.y = m->centroid.z = 0; + for (int i = 0; i < nv; i++) + m->centroid = vector3_plus(m->centroid, m->vertices.items[i]); + m->centroid = vector3_scale(1.0 / nv, m->centroid); + + /* Check if mesh is closed (every edge shared by exactly 2 faces) + and ensure consistent outward-facing normals. */ + m->is_closed = 1; + + /* Use signed volume to check/fix normal orientation. */ + double signed_vol = 0; + for (int f = 0; f < nf; f++) { + vector3 v0 = m->vertices.items[m->face_indices[3 * f]]; + vector3 v1 = m->vertices.items[m->face_indices[3 * f + 1]]; + vector3 v2 = m->vertices.items[m->face_indices[3 * f + 2]]; + signed_vol += vector3_dot(v0, vector3_cross(v1, v2)); + } + /* If signed volume is negative, normals point inward — flip all winding. */ + if (signed_vol < 0) { + for (int f = 0; f < nf; f++) { + /* Swap indices 1 and 2. */ + int tmp = m->face_indices[3 * f + 1]; + m->face_indices[3 * f + 1] = m->face_indices[3 * f + 2]; + m->face_indices[3 * f + 2] = tmp; + /* Flip normal. */ + m->face_normals[f] = vector3_scale(-1.0, m->face_normals[f]); + } + } + + /* Build BVH. */ + int max_nodes = 2 * nf; /* upper bound on BVH nodes */ + m->bvh = (mesh_bvh_node *)malloc(max_nodes * sizeof(mesh_bvh_node)); + CHECK(m->bvh, "out of memory"); + m->bvh_face_ids = (int *)malloc(nf * sizeof(int)); + CHECK(m->bvh_face_ids, "out of memory"); + for (int i = 0; i < nf; i++) + m->bvh_face_ids[i] = i; + + m->num_bvh_nodes = 0; + mesh_bvh_build(m, m->bvh_face_ids, 0, nf, m->bvh, &m->num_bvh_nodes); +} + +/***************************************************************/ +/* make_mesh constructors */ +/***************************************************************/ + +static int mesh_is_auto_center(vector3 c) { + return isnan(c.x) && isnan(c.y) && isnan(c.z); +} + +geometric_object make_mesh(material_type material, const vector3 *vertices, int num_vertices, + const int *triangles, int num_triangles) { + vector3 auto_c = {NAN, NAN, NAN}; + return make_mesh_with_center(material, auto_c, vertices, num_vertices, triangles, num_triangles); +} + +geometric_object make_mesh_with_center(material_type material, vector3 center, + const vector3 *vertices, int num_vertices, + const int *triangles, int num_triangles) { + geometric_object o = make_geometric_object(material, center); + o.which_subclass = GEOM MESH; + mesh *m = o.subclass.mesh_data = MALLOC1(mesh); + CHECK(m, "out of memory"); + memset(m, 0, sizeof(mesh)); + + /* Copy vertices. */ + m->vertices.num_items = num_vertices; + m->vertices.items = (vector3 *)malloc(num_vertices * sizeof(vector3)); + CHECK(m->vertices.items, "out of memory"); + memcpy(m->vertices.items, vertices, num_vertices * sizeof(vector3)); + + /* Copy triangle indices. */ + m->num_faces = num_triangles; + m->face_indices = (int *)malloc(3 * num_triangles * sizeof(int)); + CHECK(m->face_indices, "out of memory"); + memcpy(m->face_indices, triangles, 3 * num_triangles * sizeof(int)); + + /* Initialize derived data (normals, BVH, etc.). */ + init_mesh(&o); + + /* Set center. */ + if (mesh_is_auto_center(center)) + o.center = m->centroid; + else { + /* Shift vertices so centroid equals center. */ + vector3 shift = vector3_minus(center, m->centroid); + for (int i = 0; i < num_vertices; i++) + m->vertices.items[i] = vector3_plus(m->vertices.items[i], shift); + m->centroid = center; + o.center = center; + + /* Rebuild BVH after shifting vertices. */ + free(m->bvh); + free(m->bvh_face_ids); + m->bvh_face_ids = (int *)malloc(m->num_faces * sizeof(int)); + CHECK(m->bvh_face_ids, "out of memory"); + for (int i = 0; i < m->num_faces; i++) + m->bvh_face_ids[i] = i; + int max_nodes = 2 * m->num_faces; + m->bvh = (mesh_bvh_node *)malloc(max_nodes * sizeof(mesh_bvh_node)); + CHECK(m->bvh, "out of memory"); + m->num_bvh_nodes = 0; + mesh_bvh_build(m, m->bvh_face_ids, 0, m->num_faces, m->bvh, &m->num_bvh_nodes); + } + + return o; +} + /*************************************************************** * The remainder of this file implements geometric primitives for prisms. * A prism is a planar polygon, consisting of 3 or more user-specified diff --git a/utils/geom.scm b/utils/geom.scm index 6b03a94..4056efd 100644 --- a/utils/geom.scm +++ b/utils/geom.scm @@ -139,6 +139,16 @@ (define-property m_c2p identity_matrix 'matrix3x3) (define-property m_p2c identity_matrix 'matrix3x3)) +; A mesh is a triangulated 3D surface (closed manifold). +; Only vertices and face_indices need to be specified by the user; +; all other fields are computed internally by init_mesh in C. +; face_indices is a vector3_list where each vector3 holds 3 integer +; indices (in x, y, z) referring to entries in the vertices list. +(define-class mesh geometric-object +; fields to be filled in by users + (define-property vertices '() (make-list-type 'vector3)) + (define-property face_indices '() (make-list-type 'vector3))) + (define-class ellipsoid block (define-derived-property inverse-semi-axes 'vector3 (lambda (object) diff --git a/utils/test-mesh.c b/utils/test-mesh.c new file mode 100644 index 0000000..d3ce69c --- /dev/null +++ b/utils/test-mesh.c @@ -0,0 +1,378 @@ +/* libctl: flexible Guile-based control files for scientific software + * Copyright (C) 1998-2020 Massachusetts Institute of Technology and Steven G. Johnson + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + * + * Steven G. Johnson can be contacted at stevenj@alum.mit.edu. + */ + +/************************************************************************/ +/* test-mesh.c: unit test for mesh geometry in libctlgeom */ +/************************************************************************/ +#include +#include +#include +#include + +#include "ctlgeom.h" + +#define K_PI 3.141592653589793238462643383279502884197 +#define TOLERANCE 1e-6 + +static int test_failures = 0; + +#define ASSERT_TRUE(msg, cond) \ + do { \ + if (!(cond)) { \ + fprintf(stderr, "FAIL: %s (line %d)\n", msg, __LINE__); \ + test_failures++; \ + } \ + } while (0) + +#define ASSERT_NEAR(msg, val, expected, tol) \ + do { \ + double _v = (val), _e = (expected); \ + if (fabs(_v - _e) > (tol)) { \ + fprintf(stderr, "FAIL: %s: got %g, expected %g (line %d)\n", \ + msg, _v, _e, __LINE__); \ + test_failures++; \ + } \ + } while (0) + +/************************************************************************/ +/* Helper: create a unit cube mesh centered at origin. */ +/* 8 vertices, 12 triangles. */ +/************************************************************************/ +static geometric_object make_cube_mesh(void *material) { + vector3 verts[8] = { + {-0.5, -0.5, -0.5}, { 0.5, -0.5, -0.5}, { 0.5, 0.5, -0.5}, {-0.5, 0.5, -0.5}, + {-0.5, -0.5, 0.5}, { 0.5, -0.5, 0.5}, { 0.5, 0.5, 0.5}, {-0.5, 0.5, 0.5} + }; + /* Outward-facing triangles (right-hand rule). */ + int tris[12 * 3] = { + /* -Z face */ 0,2,1, 0,3,2, + /* +Z face */ 4,5,6, 4,6,7, + /* -Y face */ 0,1,5, 0,5,4, + /* +Y face */ 2,3,7, 2,7,6, + /* -X face */ 0,4,7, 0,7,3, + /* +X face */ 1,2,6, 1,6,5 + }; + return make_mesh(material, verts, 8, tris, 12); +} + +/************************************************************************/ +/* Helper: create a regular tetrahedron mesh. */ +/************************************************************************/ +static geometric_object make_tetra_mesh(void *material) { + double a = 1.0; + vector3 verts[4] = { + { a, a, a}, + { a, -a, -a}, + {-a, a, -a}, + {-a, -a, a} + }; + /* Outward-facing triangles. */ + int tris[4 * 3] = { + 0, 1, 2, + 0, 3, 1, + 0, 2, 3, + 1, 3, 2 + }; + return make_mesh(material, verts, 4, tris, 4); +} + +/************************************************************************/ +/* Test: point_in_mesh for cube */ +/************************************************************************/ +static void test_cube_point_in(void) { + printf("test_cube_point_in... "); + geometric_object cube = make_cube_mesh(NULL); + + /* Points inside. */ + vector3 p_inside = {0, 0, 0}; + ASSERT_TRUE("origin inside cube", point_in_fixed_pobjectp(p_inside, &cube)); + + vector3 p_inside2 = {0.4, 0.4, 0.4}; + ASSERT_TRUE("corner-near point inside cube", point_in_fixed_pobjectp(p_inside2, &cube)); + + /* Points outside. */ + vector3 p_outside = {1.0, 0, 0}; + ASSERT_TRUE("point outside cube +x", !point_in_fixed_pobjectp(p_outside, &cube)); + + vector3 p_outside2 = {0.6, 0, 0}; + ASSERT_TRUE("point outside cube +x close", !point_in_fixed_pobjectp(p_outside2, &cube)); + + vector3 p_outside3 = {0, 0.6, 0}; + ASSERT_TRUE("point outside cube +y", !point_in_fixed_pobjectp(p_outside3, &cube)); + + vector3 p_outside4 = {0, 0, -0.6}; + ASSERT_TRUE("point outside cube -z", !point_in_fixed_pobjectp(p_outside4, &cube)); + + geometric_object_destroy(cube); + printf("done\n"); +} + +/************************************************************************/ +/* Test: cube volume */ +/************************************************************************/ +static void test_cube_volume(void) { + printf("test_cube_volume... "); + geometric_object cube = make_cube_mesh(NULL); + + double vol = geom_object_volume(cube); + ASSERT_NEAR("cube volume", vol, 1.0, TOLERANCE); + + geometric_object_destroy(cube); + printf("done\n"); +} + +/************************************************************************/ +/* Test: tetrahedron volume */ +/************************************************************************/ +static void test_tetra_volume(void) { + printf("test_tetra_volume... "); + geometric_object tetra = make_tetra_mesh(NULL); + + /* Regular tetrahedron with vertices at (±1,±1,±1): + edge length = 2*sqrt(2), volume = edge^3 / (6*sqrt(2)) = 8/3 */ + double expected_vol = 8.0 / 3.0; + double vol = geom_object_volume(tetra); + ASSERT_NEAR("tetra volume", vol, expected_vol, TOLERANCE); + + geometric_object_destroy(tetra); + printf("done\n"); +} + +/************************************************************************/ +/* Test: bounding box */ +/************************************************************************/ +static void test_cube_bounding_box(void) { + printf("test_cube_bounding_box... "); + geometric_object cube = make_cube_mesh(NULL); + + geom_box box; + geom_get_bounding_box(cube, &box); + ASSERT_NEAR("cube bbox low.x", box.low.x, -0.5, TOLERANCE); + ASSERT_NEAR("cube bbox low.y", box.low.y, -0.5, TOLERANCE); + ASSERT_NEAR("cube bbox low.z", box.low.z, -0.5, TOLERANCE); + ASSERT_NEAR("cube bbox high.x", box.high.x, 0.5, TOLERANCE); + ASSERT_NEAR("cube bbox high.y", box.high.y, 0.5, TOLERANCE); + ASSERT_NEAR("cube bbox high.z", box.high.z, 0.5, TOLERANCE); + + geometric_object_destroy(cube); + printf("done\n"); +} + +/************************************************************************/ +/* Test: normals */ +/************************************************************************/ +static void test_cube_normals(void) { + printf("test_cube_normals... "); + geometric_object cube = make_cube_mesh(NULL); + + /* Point near +X face: normal should point in +X direction. */ + vector3 p = {0.49, 0, 0}; + vector3 n = normal_to_fixed_object(p, cube); + ASSERT_NEAR("normal near +x face: nx", fabs(n.x), 1.0, 0.01); + ASSERT_NEAR("normal near +x face: ny", fabs(n.y), 0.0, 0.01); + ASSERT_NEAR("normal near +x face: nz", fabs(n.z), 0.0, 0.01); + + /* Point near +Y face. */ + vector3 p2 = {0, 0.49, 0}; + vector3 n2 = normal_to_fixed_object(p2, cube); + ASSERT_NEAR("normal near +y face: ny", fabs(n2.y), 1.0, 0.01); + + /* Point near -Z face. */ + vector3 p3 = {0, 0, -0.49}; + vector3 n3 = normal_to_fixed_object(p3, cube); + ASSERT_NEAR("normal near -z face: nz", fabs(n3.z), 1.0, 0.01); + + geometric_object_destroy(cube); + printf("done\n"); +} + +/************************************************************************/ +/* Test: cube mesh vs block (compare point_in results) */ +/************************************************************************/ +static void test_cube_vs_block(void) { + printf("test_cube_vs_block... "); + geometric_object cube_mesh = make_cube_mesh(NULL); + vector3 center = {0, 0, 0}; + vector3 e1 = {1, 0, 0}, e2 = {0, 1, 0}, e3 = {0, 0, 1}; + vector3 size = {1, 1, 1}; + geometric_object cube_block = make_block(NULL, center, e1, e2, e3, size); + + int mismatches = 0; + srand(42); + for (int i = 0; i < 10000; i++) { + vector3 p; + p.x = (rand() / (double)RAND_MAX) * 2.0 - 1.0; + p.y = (rand() / (double)RAND_MAX) * 2.0 - 1.0; + p.z = (rand() / (double)RAND_MAX) * 2.0 - 1.0; + + int in_mesh = point_in_fixed_pobjectp(p, &cube_mesh); + int in_block = point_in_fixed_pobjectp(p, &cube_block); + + /* Allow mismatches very close to the surface (within 1e-6). */ + if (in_mesh != in_block) { + double dx = fmin(fabs(fabs(p.x) - 0.5), fmin(fabs(fabs(p.y) - 0.5), fabs(fabs(p.z) - 0.5))); + if (dx > 1e-6) mismatches++; + } + } + + ASSERT_TRUE("cube mesh matches block for random points", mismatches == 0); + + geometric_object_destroy(cube_mesh); + geometric_object_destroy(cube_block); + printf("done (%d mismatches)\n", mismatches); +} + +/************************************************************************/ +/* Test: line_segment intersection */ +/************************************************************************/ +static void test_cube_line_segment(void) { + printf("test_cube_line_segment... "); + geometric_object cube = make_cube_mesh(NULL); + + /* Line along x-axis through center: should intersect for length 1.0. */ + vector3 p = {0, 0, 0}; + vector3 d = {1, 0, 0}; + double len = intersect_line_segment_with_object(p, d, cube, -1.0, 1.0); + ASSERT_NEAR("x-axis segment length", len, 1.0, 0.01); + + /* Line along y-axis. */ + vector3 d2 = {0, 1, 0}; + double len2 = intersect_line_segment_with_object(p, d2, cube, -1.0, 1.0); + ASSERT_NEAR("y-axis segment length", len2, 1.0, 0.01); + + /* Line that misses the cube. */ + vector3 p3 = {2, 2, 2}; + double len3 = intersect_line_segment_with_object(p3, d, cube, -1.0, 1.0); + ASSERT_NEAR("miss segment length", len3, 0.0, 0.01); + + geometric_object_destroy(cube); + printf("done\n"); +} + +/************************************************************************/ +/* Test: copy and destroy */ +/************************************************************************/ +static void test_copy_destroy(void) { + printf("test_copy_destroy... "); + geometric_object cube = make_cube_mesh(NULL); + geometric_object cube_copy; + geometric_object_copy(&cube, &cube_copy); + + /* Verify copy works (point-in test). */ + vector3 p = {0, 0, 0}; + ASSERT_TRUE("copy: origin inside", point_in_fixed_pobjectp(p, &cube_copy)); + + vector3 p2 = {1, 0, 0}; + ASSERT_TRUE("copy: outside", !point_in_fixed_pobjectp(p2, &cube_copy)); + + double vol = geom_object_volume(cube_copy); + ASSERT_NEAR("copy: volume", vol, 1.0, TOLERANCE); + + geometric_object_destroy(cube); + geometric_object_destroy(cube_copy); + printf("done\n"); +} + +/************************************************************************/ +/* Test: display info */ +/************************************************************************/ +static void test_display_info(void) { + printf("test_display_info... "); + geometric_object cube = make_cube_mesh(NULL); + display_geometric_object_info(0, cube); + geometric_object_destroy(cube); + printf("done\n"); +} + +/************************************************************************/ +/* Test: point_in_mesh for tetrahedron */ +/************************************************************************/ +static void test_tetra_point_in(void) { + printf("test_tetra_point_in... "); + geometric_object tetra = make_tetra_mesh(NULL); + + /* Centroid (0,0,0) should be inside. */ + vector3 p_in = {0, 0, 0}; + ASSERT_TRUE("tetra: centroid inside", point_in_fixed_pobjectp(p_in, &tetra)); + + /* Point far outside. */ + vector3 p_out = {3, 3, 3}; + ASSERT_TRUE("tetra: far point outside", !point_in_fixed_pobjectp(p_out, &tetra)); + + geometric_object_destroy(tetra); + printf("done\n"); +} + +/************************************************************************/ +/* Test: make_mesh_with_center */ +/************************************************************************/ +static void test_mesh_with_center(void) { + printf("test_mesh_with_center... "); + vector3 verts[8] = { + {-0.5, -0.5, -0.5}, { 0.5, -0.5, -0.5}, { 0.5, 0.5, -0.5}, {-0.5, 0.5, -0.5}, + {-0.5, -0.5, 0.5}, { 0.5, -0.5, 0.5}, { 0.5, 0.5, 0.5}, {-0.5, 0.5, 0.5} + }; + int tris[12 * 3] = { + 0,2,1, 0,3,2, + 4,5,6, 4,6,7, + 0,1,5, 0,5,4, + 2,3,7, 2,7,6, + 0,4,7, 0,7,3, + 1,2,6, 1,6,5 + }; + vector3 center = {5, 5, 5}; + geometric_object cube = make_mesh_with_center(NULL, center, verts, 8, tris, 12); + + /* Point at center (5,5,5) should be inside. */ + vector3 p_in = {5, 5, 5}; + ASSERT_TRUE("centered mesh: center inside", point_in_fixed_pobjectp(p_in, &cube)); + + /* Point at origin should be outside. */ + vector3 p_out = {0, 0, 0}; + ASSERT_TRUE("centered mesh: origin outside", !point_in_fixed_pobjectp(p_out, &cube)); + + /* Volume should still be 1. */ + double vol = geom_object_volume(cube); + ASSERT_NEAR("centered mesh volume", vol, 1.0, TOLERANCE); + + geometric_object_destroy(cube); + printf("done\n"); +} + +/************************************************************************/ +int main(void) { + geom_initialize(); + + test_cube_point_in(); + test_cube_volume(); + test_tetra_volume(); + test_cube_bounding_box(); + test_cube_normals(); + test_cube_vs_block(); + test_cube_line_segment(); + test_copy_destroy(); + test_display_info(); + test_tetra_point_in(); + test_mesh_with_center(); + + printf("\n%d test failures\n", test_failures); + return test_failures > 0 ? 1 : 0; +} From bdae298e7d641c9b1051e96535beec3c496a329a Mon Sep 17 00:00:00 2001 From: Luochenghuang Date: Tue, 7 Apr 2026 13:18:07 -0700 Subject: [PATCH 2/8] Fix point_in_mesh robustness and reinit_mesh performance Two fixes: 1. point_in_mesh: use 3-ray majority vote instead of a single ray. A single ray can pass through a shared mesh edge, causing the deduplication to either miss or double-count a crossing and flip the inside/outside parity. Three rays in irrational directions make it extremely unlikely that more than one hits a degenerate case at the same query point. 2. reinit_mesh: skip BVH rebuild if already initialized. Unlike blocks/cylinders, mesh uses absolute coordinates and doesn't depend on the lattice basis. geom_fix_object_ptr is called hundreds of times during meep's init_sim via geometry copies; without the guard, each call rebuilt the BVH (~150ms for 9K triangles), causing ~40s of wasted time. --- utils/geom.c | 38 ++++++++++++++++++++++++++++++++------ 1 file changed, 32 insertions(+), 6 deletions(-) diff --git a/utils/geom.c b/utils/geom.c index 82d8eec..825e7ee 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -83,6 +83,7 @@ static void display_prism_info(int indentby, geometric_object *o); static void init_prism(geometric_object *o); static void reinit_prism(geometric_object *o); static void init_mesh(geometric_object *o); +static void reinit_mesh(geometric_object *o); static boolean point_in_mesh(const mesh *m, vector3 p); static vector3 normal_to_mesh(const mesh *m, vector3 p); static void get_mesh_bounding_box(const mesh *m, geom_box *box); @@ -183,7 +184,7 @@ void geom_fix_object_ptr(geometric_object *o) { break; } case GEOM MESH: { - init_mesh(o); + reinit_mesh(o); break; } case GEOM GEOMETRIC_OBJECT_SELF: @@ -2452,11 +2453,22 @@ static int count_ray_mesh_intersections(const mesh *m, vector3 origin, vector3 d static boolean point_in_mesh(const mesh *m, vector3 p) { if (!m->is_closed) return 0; - /* Cast a ray along a non-axis-aligned direction to reduce chance of - hitting mesh edges exactly (where deduplication is most fragile). */ - vector3 dir = {0.57735026918962576, 0.57735026918962576, 0.57735026918962576}; - int crossings = count_ray_mesh_intersections(m, p, dir); - return (crossings % 2) == 1; + /* Cast rays along multiple directions and take majority vote. + A single ray can give wrong parity when it passes through a mesh + edge or vertex, causing deduplication to miss or double-count + a crossing. Three irrational directions make it extremely unlikely + that more than one ray hits a degenerate case at the same point. */ + static const vector3 dirs[3] = { + {0.57735026918962576, 0.57735026918962576, 0.57735026918962576}, /* (1,1,1)/√3 */ + {0.80178372573727319, 0.53452248382484879, 0.26726124191242440}, /* (3,2,1)/√14 */ + {0.12309149097933272, 0.49236596391733088, 0.86164043685532904}, /* (1,4,7)/√66 */ + }; + int votes = 0; + for (int d = 0; d < 3; d++) { + int crossings = count_ray_mesh_intersections(m, p, dirs[d]); + if (crossings % 2 == 1) votes++; + } + return votes >= 2; } static vector3 normal_to_mesh(const mesh *m, vector3 p) { @@ -2622,6 +2634,20 @@ static void init_mesh(geometric_object *o) { /* make_mesh constructors */ /***************************************************************/ +static void reinit_mesh(geometric_object *o) { + mesh *m = o->subclass.mesh_data; + /* Skip if already initialized — mesh uses absolute coordinates + and doesn't depend on the lattice basis, unlike blocks/cylinders. */ + if (m->bvh != NULL) return; + free(m->face_normals); + free(m->face_areas); + free(m->bvh_face_ids); + m->face_normals = NULL; + m->face_areas = NULL; + m->bvh_face_ids = NULL; + init_mesh(o); +} + static int mesh_is_auto_center(vector3 c) { return isnan(c.x) && isnan(c.y) && isnan(c.z); } From 32451aac2cf9d9841f8b6b8e0a699bfea45405e7 Mon Sep 17 00:00:00 2001 From: Luochenghuang Date: Tue, 7 Apr 2026 13:40:11 -0700 Subject: [PATCH 3/8] Lazy majority vote and mesh closure validation - point_in_mesh: cast 1 ray first, only fall back to 3-ray majority vote if deduplication removed near-duplicate t-values (indicating an edge/vertex hit). 3x faster for the common case (no edge hits). - init_mesh: validate mesh closure by building an edge-to-face-count hash table. Every edge must be shared by exactly 2 faces. Sets is_closed=false and prints a warning for open/non-manifold meshes. point_in_mesh returns false for all points on open meshes. - Add 2 new tests: open mesh detection, closed mesh detection (13 tests total, all passing). --- utils/geom.c | 111 ++++++++++++++++++++++++++++++++++++++-------- utils/test-mesh.c | 47 ++++++++++++++++++++ 2 files changed, 139 insertions(+), 19 deletions(-) diff --git a/utils/geom.c b/utils/geom.c index 825e7ee..e70bbb2 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -2427,8 +2427,11 @@ static int remove_duplicate_intersections(double *slist, int count, double tol) #define MESH_MAX_INTERSECTIONS 256 /* Count ray-mesh intersections using BVH traversal. - Collects all intersection t-values, deduplicates, and returns the count. */ -static int count_ray_mesh_intersections(const mesh *m, vector3 origin, vector3 dir) { + Collects all intersection t-values, deduplicates, and returns the count. + If had_degenerate is non-NULL, sets it to 1 if any near-duplicate + t-values were found (indicating a potential edge/vertex hit). */ +static int count_ray_mesh_intersections_ex(const mesh *m, vector3 origin, vector3 dir, + int *had_degenerate) { double slist[MESH_MAX_INTERSECTIONS]; int count = mesh_ray_all_intersections(m, origin, dir, slist, MESH_MAX_INTERSECTIONS); @@ -2440,32 +2443,46 @@ static int count_ray_mesh_intersections(const mesh *m, vector3 origin, vector3 d if (nforward > 1) { qsort(slist, nforward, sizeof(double), mesh_dcmp); - nforward = remove_duplicate_intersections(slist, nforward, 1e-10); + int deduped = remove_duplicate_intersections(slist, nforward, 1e-10); + if (had_degenerate) *had_degenerate = (deduped != nforward); + nforward = deduped; + } else { + if (had_degenerate) *had_degenerate = 0; } return nforward; } +static int count_ray_mesh_intersections(const mesh *m, vector3 origin, vector3 dir) { + return count_ray_mesh_intersections_ex(m, origin, dir, NULL); +} + /***************************************************************/ /* Core mesh geometric operations */ /***************************************************************/ +static const vector3 mesh_ray_dirs[3] = { + {0.57735026918962576, 0.57735026918962576, 0.57735026918962576}, /* (1,1,1)/√3 */ + {0.80178372573727319, 0.53452248382484879, 0.26726124191242440}, /* (3,2,1)/√14 */ + {0.12309149097933272, 0.49236596391733088, 0.86164043685532904}, /* (1,4,7)/√66 */ +}; + static boolean point_in_mesh(const mesh *m, vector3 p) { if (!m->is_closed) return 0; - /* Cast rays along multiple directions and take majority vote. - A single ray can give wrong parity when it passes through a mesh - edge or vertex, causing deduplication to miss or double-count - a crossing. Three irrational directions make it extremely unlikely - that more than one ray hits a degenerate case at the same point. */ - static const vector3 dirs[3] = { - {0.57735026918962576, 0.57735026918962576, 0.57735026918962576}, /* (1,1,1)/√3 */ - {0.80178372573727319, 0.53452248382484879, 0.26726124191242440}, /* (3,2,1)/√14 */ - {0.12309149097933272, 0.49236596391733088, 0.86164043685532904}, /* (1,4,7)/√66 */ - }; - int votes = 0; - for (int d = 0; d < 3; d++) { - int crossings = count_ray_mesh_intersections(m, p, dirs[d]); + /* Fast path: cast one ray. If no degenerate edge/vertex hits were + detected during deduplication, trust the result immediately. + Only fall back to 3-ray majority vote when a degenerate case + is detected (near-duplicate t-values removed). */ + int had_degenerate; + int crossings = count_ray_mesh_intersections_ex(m, p, mesh_ray_dirs[0], &had_degenerate); + if (!had_degenerate) + return (crossings % 2) == 1; + + /* Degenerate hit detected — cast two more rays and take majority vote. */ + int votes = (crossings % 2 == 1) ? 1 : 0; + for (int d = 1; d < 3; d++) { + crossings = count_ray_mesh_intersections(m, p, mesh_ray_dirs[d]); if (crossings % 2 == 1) votes++; } return votes >= 2; @@ -2593,9 +2610,65 @@ static void init_mesh(geometric_object *o) { m->centroid = vector3_plus(m->centroid, m->vertices.items[i]); m->centroid = vector3_scale(1.0 / nv, m->centroid); - /* Check if mesh is closed (every edge shared by exactly 2 faces) - and ensure consistent outward-facing normals. */ - m->is_closed = 1; + /* Check if mesh is closed: every edge must be shared by exactly 2 faces. + An edge is identified by a sorted pair of vertex indices. We use a + hash table mapping edge → face count. */ + { + /* Hash table size: next power of 2 above 3*nf (each face has 3 edges). */ + int htsize = 1; + while (htsize < 4 * nf) htsize <<= 1; + int htmask = htsize - 1; + + /* Each bucket stores (v_lo, v_hi, count). Use open addressing. */ + int *ht_vlo = (int *)calloc(htsize, sizeof(int)); + int *ht_vhi = (int *)calloc(htsize, sizeof(int)); + int *ht_cnt = (int *)calloc(htsize, sizeof(int)); + /* Mark empty slots with -1. */ + for (int i = 0; i < htsize; i++) ht_vlo[i] = -1; + + m->is_closed = 1; + for (int f = 0; f < nf && m->is_closed; f++) { + for (int e = 0; e < 3; e++) { + int va = m->face_indices[3 * f + e]; + int vb = m->face_indices[3 * f + (e + 1) % 3]; + int vlo = (va < vb) ? va : vb; + int vhi = (va < vb) ? vb : va; + unsigned int h = (unsigned int)(vlo * 73856093u ^ vhi * 19349663u) & htmask; + /* Linear probing. */ + while (1) { + if (ht_vlo[h] == -1) { + /* Empty slot: insert new edge. */ + ht_vlo[h] = vlo; + ht_vhi[h] = vhi; + ht_cnt[h] = 1; + break; + } else if (ht_vlo[h] == vlo && ht_vhi[h] == vhi) { + /* Found existing edge. */ + ht_cnt[h]++; + if (ht_cnt[h] > 2) m->is_closed = 0; /* non-manifold edge */ + break; + } + h = (h + 1) & htmask; + } + } + } + /* Check that all edges have exactly 2 faces. */ + if (m->is_closed) { + for (int i = 0; i < htsize; i++) { + if (ht_vlo[i] != -1 && ht_cnt[i] != 2) { + m->is_closed = 0; + break; + } + } + } + free(ht_vlo); + free(ht_vhi); + free(ht_cnt); + + if (!m->is_closed) + ctl_printf("WARNING: mesh is not closed (not all edges shared by exactly 2 faces).\n" + " point_in_mesh results may be incorrect.\n"); + } /* Use signed volume to check/fix normal orientation. */ double signed_vol = 0; diff --git a/utils/test-mesh.c b/utils/test-mesh.c index d3ce69c..8f500d0 100644 --- a/utils/test-mesh.c +++ b/utils/test-mesh.c @@ -357,6 +357,51 @@ static void test_mesh_with_center(void) { printf("done\n"); } +/************************************************************************/ +/* Test: open mesh detection */ +/************************************************************************/ +static void test_open_mesh(void) { + printf("test_open_mesh... "); + /* A tetrahedron with one face duplicated and one face missing — not closed. + 4 faces total (passes the >= 4 check) but has boundary edges. */ + vector3 verts[4] = { + {0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1} + }; + int tris[4 * 3] = { + 0, 1, 2, + 0, 1, 3, + 0, 2, 3, + 0, 2, 3 /* duplicate of face 2 instead of the missing face (1,2,3) */ + }; + geometric_object open = make_mesh(NULL, verts, 4, tris, 4); + mesh *m = open.subclass.mesh_data; + ASSERT_TRUE("open mesh detected as not closed", !m->is_closed); + + /* point_in_mesh should return false for open meshes. */ + vector3 p = {0.1, 0.1, 0.1}; + ASSERT_TRUE("open mesh: point_in returns false", !point_in_fixed_pobjectp(p, &open)); + + geometric_object_destroy(open); + printf("done\n"); +} + +/************************************************************************/ +/* Test: closed mesh correctly detected */ +/************************************************************************/ +static void test_closed_detection(void) { + printf("test_closed_detection... "); + geometric_object cube = make_cube_mesh(NULL); + mesh *m = cube.subclass.mesh_data; + ASSERT_TRUE("cube mesh detected as closed", m->is_closed); + geometric_object_destroy(cube); + + geometric_object tetra = make_tetra_mesh(NULL); + m = tetra.subclass.mesh_data; + ASSERT_TRUE("tetra mesh detected as closed", m->is_closed); + geometric_object_destroy(tetra); + printf("done\n"); +} + /************************************************************************/ int main(void) { geom_initialize(); @@ -372,6 +417,8 @@ int main(void) { test_display_info(); test_tetra_point_in(); test_mesh_with_center(); + test_open_mesh(); + test_closed_detection(); printf("\n%d test failures\n", test_failures); return test_failures > 0 ? 1 : 0; From 573183da532e42d6a32ad49b1fd8b5c58f105aba Mon Sep 17 00:00:00 2001 From: Luochenghuang Date: Tue, 7 Apr 2026 13:50:45 -0700 Subject: [PATCH 4/8] Add closure validation edge case tests Test boundary edges (1 face per edge), non-manifold edges (4 faces sharing an edge), and isolated vertices (vertex unreferenced by any face). 16 tests total, all passing. --- utils/test-mesh.c | 84 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 84 insertions(+) diff --git a/utils/test-mesh.c b/utils/test-mesh.c index 8f500d0..164bbf5 100644 --- a/utils/test-mesh.c +++ b/utils/test-mesh.c @@ -402,6 +402,87 @@ static void test_closed_detection(void) { printf("done\n"); } +/************************************************************************/ +/* Test: boundary edge (1 face on an edge) */ +/************************************************************************/ +static void test_boundary_edge(void) { + printf("test_boundary_edge... "); + /* Take a valid tetrahedron but remove the last face, replace with + a face that shares no edge with it — leaving boundary edges + where the removed face was. */ + vector3 verts[5] = { + {0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {2, 2, 2} + }; + int tris[4 * 3] = { + 0, 1, 2, /* face 0 */ + 0, 1, 3, /* face 1 */ + 0, 2, 3, /* face 2 */ + /* missing: 1, 2, 3 — replaced with unconnected face */ + 0, 1, 4 /* face 3: edge (1,2) and (2,3) now have only 1 face */ + }; + geometric_object obj = make_mesh(NULL, verts, 5, tris, 4); + mesh *m = obj.subclass.mesh_data; + ASSERT_TRUE("boundary edge mesh detected as not closed", !m->is_closed); + geometric_object_destroy(obj); + printf("done\n"); +} + +/************************************************************************/ +/* Test: non-manifold edge (3+ faces sharing an edge) */ +/************************************************************************/ +static void test_nonmanifold_edge(void) { + printf("test_nonmanifold_edge... "); + /* A tetrahedron (4 faces) plus 2 extra faces sharing edge (0,1), + giving that edge 4 adjacent faces instead of 2. */ + vector3 verts[6] = { + {0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}, + {0.5, -1, 0.5}, {0.5, -1, -0.5} + }; + int tris[6 * 3] = { + 0, 1, 2, /* tetra face 0 — uses edge (0,1) */ + 0, 1, 3, /* tetra face 1 — uses edge (0,1) */ + 0, 2, 3, /* tetra face 2 */ + 1, 2, 3, /* tetra face 3 */ + 0, 1, 4, /* extra face — edge (0,1) now has 3 faces */ + 0, 1, 5, /* extra face — edge (0,1) now has 4 faces */ + }; + geometric_object obj = make_mesh(NULL, verts, 6, tris, 6); + mesh *m = obj.subclass.mesh_data; + ASSERT_TRUE("non-manifold edge mesh detected as not closed", !m->is_closed); + geometric_object_destroy(obj); + printf("done\n"); +} + +/************************************************************************/ +/* Test: isolated vertex (vertex not referenced by any face) */ +/************************************************************************/ +static void test_isolated_vertex(void) { + printf("test_isolated_vertex... "); + /* A valid tetrahedron plus an extra vertex that no face references. + The mesh should still be detected as closed (isolated vertices + don't affect edge manifold-ness). */ + vector3 verts[5] = { + {1, 1, 1}, {1, -1, -1}, {-1, 1, -1}, {-1, -1, 1}, + {99, 99, 99} /* isolated vertex */ + }; + int tris[4 * 3] = { + 0, 1, 2, + 0, 3, 1, + 0, 2, 3, + 1, 3, 2 + }; + geometric_object obj = make_mesh(NULL, verts, 5, tris, 4); + mesh *m = obj.subclass.mesh_data; + ASSERT_TRUE("mesh with isolated vertex still closed", m->is_closed); + + /* Point inside tetrahedron should still work. */ + vector3 p = {0, 0, 0}; + ASSERT_TRUE("isolated vertex: centroid inside", point_in_fixed_pobjectp(p, &obj)); + + geometric_object_destroy(obj); + printf("done\n"); +} + /************************************************************************/ int main(void) { geom_initialize(); @@ -419,6 +500,9 @@ int main(void) { test_mesh_with_center(); test_open_mesh(); test_closed_detection(); + test_boundary_edge(); + test_nonmanifold_edge(); + test_isolated_vertex(); printf("\n%d test failures\n", test_failures); return test_failures > 0 ? 1 : 0; From d9043b60988a68181a804a3b3a57f22eca8be5b3 Mon Sep 17 00:00:00 2001 From: Luochenghuang Date: Tue, 7 Apr 2026 15:24:58 -0700 Subject: [PATCH 5/8] Eliminate redundant BVH build and point_in_mesh fallback MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 1. make_mesh_with_center: shift vertices before calling init_mesh so the BVH is only built once (was building twice when center was specified). 2. intersect_line_segment_with_mesh: determine inside/outside at the segment start from the full intersection list parity (count crossings before parameter 'a') instead of calling point_in_mesh as a fallback. Eliminates millions of extra ray casts during subpixel smoothing. ~20% faster set_epsilon on the Utah teapot (28.6s → 22.9s). --- utils/geom.c | 78 +++++++++++++++++++++++----------------------------- 1 file changed, 34 insertions(+), 44 deletions(-) diff --git a/utils/geom.c b/utils/geom.c index e70bbb2..8dacc12 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -2542,27 +2542,27 @@ static double intersect_line_segment_with_mesh(const mesh *m, vector3 p, vector3 double slist[MESH_MAX_INTERSECTIONS]; int num = intersect_line_with_mesh(m, p, d, slist); - /* Find first intersection > a. */ - int na = -1; - for (int i = 0; na == -1 && i < num; i++) - if (slist[i] > a) na = i; - if (na == -1) { - /* No intersections after a. Check if p+a*d is inside. */ - vector3 pt = vector3_plus(p, vector3_scale(a, d)); - if (point_in_mesh(m, pt)) - return b - a; - return 0.0; - } - - /* Determine if we start inside (even-odd parity). */ - int inside = ((na % 2) == 0) ? 0 : 1; + /* The sorted intersection list gives all surface crossings along the + full ray. At t=-inf we are outside, so the parity after k crossings + tells us inside/outside: odd = inside, even = outside. + + Count crossings before 'a' to determine if we start inside [a,b]. + Then walk crossings within [a,b] toggling parity and accumulating + interior length. No point_in_mesh fallback needed. */ + int crossings_before_a = 0; + for (int i = 0; i < num && slist[i] <= a; i++) + crossings_before_a++; + + int inside = (crossings_before_a % 2 == 1); double last_s = a, ds = 0.0; - for (int i = na; i < num; i++) { - double this_s = fmin(b, slist[i]); - if (inside) ds += (this_s - last_s); - if (b <= slist[i]) break; - inside = 1 - inside; - last_s = this_s; + for (int i = crossings_before_a; i < num; i++) { + if (slist[i] >= b) { + if (inside) ds += (b - last_s); + break; + } + if (inside) ds += (slist[i] - last_s); + inside = !inside; + last_s = slist[i]; } if (inside && last_s < b) ds += (b - last_s); return ds > 0.0 ? ds : 0.0; @@ -2752,34 +2752,24 @@ geometric_object make_mesh_with_center(material_type material, vector3 center, CHECK(m->face_indices, "out of memory"); memcpy(m->face_indices, triangles, 3 * num_triangles * sizeof(int)); - /* Initialize derived data (normals, BVH, etc.). */ - init_mesh(&o); - - /* Set center. */ - if (mesh_is_auto_center(center)) - o.center = m->centroid; - else { - /* Shift vertices so centroid equals center. */ - vector3 shift = vector3_minus(center, m->centroid); + /* Shift vertices before init so BVH is only built once. */ + if (!mesh_is_auto_center(center)) { + /* Compute centroid to determine the shift. */ + vector3 centroid = {0, 0, 0}; + for (int i = 0; i < num_vertices; i++) + centroid = vector3_plus(centroid, m->vertices.items[i]); + centroid = vector3_scale(1.0 / num_vertices, centroid); + vector3 shift = vector3_minus(center, centroid); for (int i = 0; i < num_vertices; i++) m->vertices.items[i] = vector3_plus(m->vertices.items[i], shift); - m->centroid = center; - o.center = center; - - /* Rebuild BVH after shifting vertices. */ - free(m->bvh); - free(m->bvh_face_ids); - m->bvh_face_ids = (int *)malloc(m->num_faces * sizeof(int)); - CHECK(m->bvh_face_ids, "out of memory"); - for (int i = 0; i < m->num_faces; i++) - m->bvh_face_ids[i] = i; - int max_nodes = 2 * m->num_faces; - m->bvh = (mesh_bvh_node *)malloc(max_nodes * sizeof(mesh_bvh_node)); - CHECK(m->bvh, "out of memory"); - m->num_bvh_nodes = 0; - mesh_bvh_build(m, m->bvh_face_ids, 0, m->num_faces, m->bvh, &m->num_bvh_nodes); } + /* Initialize derived data (normals, closure check, BVH). */ + init_mesh(&o); + + /* Set center from the (possibly shifted) centroid. */ + o.center = m->centroid; + return o; } From 0c458ec706dfcca69a1e34464f4bc9594a60c4e4 Mon Sep 17 00:00:00 2001 From: Luochenghuang Date: Tue, 7 Apr 2026 15:27:08 -0700 Subject: [PATCH 6/8] Minor cleanups: remove dead code, optimize BVH queries - Remove unused bvh_node_box helper function - find_closest_face: push nearer BVH child last so it's popped first, improving pruning and reducing nodes visited for normal queries - Extract mesh_triangle_vertices helper to avoid redundant vertex fetches. Add mesh_triangle_bbox_centroid to compute both in one pass during BVH binning (was fetching 3 vertices twice per face) --- utils/geom.c | 75 ++++++++++++++++++++++++++++++++++------------------ 1 file changed, 50 insertions(+), 25 deletions(-) diff --git a/utils/geom.c b/utils/geom.c index 8dacc12..8af9477 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -1935,28 +1935,24 @@ geometric_object make_ellipsoid(material_type material, vector3 center, vector3 /* Triangle utilities */ /***************************************************************/ -/* Helper to get a geom_box from a BVH node. */ -static geom_box bvh_node_box(const mesh_bvh_node *node) { - geom_box b; - b.low = node->bbox_low; - b.high = node->bbox_high; - return b; -} - /* Helper to set a BVH node's bbox from a geom_box. */ static void bvh_node_set_box(mesh_bvh_node *node, const geom_box *box) { node->bbox_low = box->low; node->bbox_high = box->high; } +/* Fetch the three vertices of a triangle. */ +static void mesh_triangle_vertices(const mesh *m, int face_id, + vector3 *v0, vector3 *v1, vector3 *v2) { + *v0 = m->vertices.items[m->face_indices[3 * face_id]]; + *v1 = m->vertices.items[m->face_indices[3 * face_id + 1]]; + *v2 = m->vertices.items[m->face_indices[3 * face_id + 2]]; +} + /* Compute the AABB of a single triangle. */ static void mesh_triangle_bbox(const mesh *m, int face_id, geom_box *box) { - int i0 = m->face_indices[3 * face_id]; - int i1 = m->face_indices[3 * face_id + 1]; - int i2 = m->face_indices[3 * face_id + 2]; - vector3 v0 = m->vertices.items[i0]; - vector3 v1 = m->vertices.items[i1]; - vector3 v2 = m->vertices.items[i2]; + vector3 v0, v1, v2; + mesh_triangle_vertices(m, face_id, &v0, &v1, &v2); box->low.x = fmin(v0.x, fmin(v1.x, v2.x)); box->low.y = fmin(v0.y, fmin(v1.y, v2.y)); box->low.z = fmin(v0.z, fmin(v1.z, v2.z)); @@ -1967,12 +1963,8 @@ static void mesh_triangle_bbox(const mesh *m, int face_id, geom_box *box) { /* Compute the centroid of a triangle. */ static vector3 mesh_triangle_centroid(const mesh *m, int face_id) { - int i0 = m->face_indices[3 * face_id]; - int i1 = m->face_indices[3 * face_id + 1]; - int i2 = m->face_indices[3 * face_id + 2]; - vector3 v0 = m->vertices.items[i0]; - vector3 v1 = m->vertices.items[i1]; - vector3 v2 = m->vertices.items[i2]; + vector3 v0, v1, v2; + mesh_triangle_vertices(m, face_id, &v0, &v1, &v2); vector3 c; c.x = (v0.x + v1.x + v2.x) / 3.0; c.y = (v0.y + v1.y + v2.y) / 3.0; @@ -1980,6 +1972,22 @@ static vector3 mesh_triangle_centroid(const mesh *m, int face_id) { return c; } +/* Compute the AABB and centroid of a triangle in one pass. */ +static void mesh_triangle_bbox_centroid(const mesh *m, int face_id, + geom_box *box, vector3 *centroid) { + vector3 v0, v1, v2; + mesh_triangle_vertices(m, face_id, &v0, &v1, &v2); + box->low.x = fmin(v0.x, fmin(v1.x, v2.x)); + box->low.y = fmin(v0.y, fmin(v1.y, v2.y)); + box->low.z = fmin(v0.z, fmin(v1.z, v2.z)); + box->high.x = fmax(v0.x, fmax(v1.x, v2.x)); + box->high.y = fmax(v0.y, fmax(v1.y, v2.y)); + box->high.z = fmax(v0.z, fmax(v1.z, v2.z)); + centroid->x = (v0.x + v1.x + v2.x) / 3.0; + centroid->y = (v0.y + v1.y + v2.y) / 3.0; + centroid->z = (v0.z + v1.z + v2.z) / 3.0; +} + /* Compute the surface area of a geom_box. */ static double geom_box_surface_area(const geom_box *b) { double dx = b->high.x - b->low.x; @@ -2047,7 +2055,9 @@ static int mesh_bvh_build(mesh *m, int *face_ids, int start, int count, double inv_range = MESH_BVH_NUM_BINS / (hi - lo); for (int i = 0; i < count; i++) { - vector3 c = mesh_triangle_centroid(m, face_ids[start + i]); + geom_box fb; + vector3 c; + mesh_triangle_bbox_centroid(m, face_ids[start + i], &fb, &c); double val; switch (axis) { case 0: val = c.x; break; @@ -2058,8 +2068,6 @@ static int mesh_bvh_build(mesh *m, int *face_ids, int start, int count, if (bin < 0) bin = 0; if (bin >= MESH_BVH_NUM_BINS) bin = MESH_BVH_NUM_BINS - 1; bin_counts[bin]++; - geom_box fb; - mesh_triangle_bbox(m, face_ids[start + i], &fb); geom_box_union(&bin_boxes[bin], &bin_boxes[bin], &fb); } @@ -2354,8 +2362,25 @@ static int find_closest_face(const mesh *m, vector3 p, double *dist2) { } } else { if (stack_top + 2 > 64) continue; - stack[stack_top++] = node->left_child; - stack[stack_top++] = node->right_child; + /* Push the farther child first so the nearer child is popped first, + giving better pruning of the farther subtree. */ + const mesh_bvh_node *left = &m->bvh[node->left_child]; + const mesh_bvh_node *right = &m->bvh[node->right_child]; + double ldx = fmax(0, fmax(left->bbox_low.x - p.x, p.x - left->bbox_high.x)); + double ldy = fmax(0, fmax(left->bbox_low.y - p.y, p.y - left->bbox_high.y)); + double ldz = fmax(0, fmax(left->bbox_low.z - p.z, p.z - left->bbox_high.z)); + double rdx = fmax(0, fmax(right->bbox_low.x - p.x, p.x - right->bbox_high.x)); + double rdy = fmax(0, fmax(right->bbox_low.y - p.y, p.y - right->bbox_high.y)); + double rdz = fmax(0, fmax(right->bbox_low.z - p.z, p.z - right->bbox_high.z)); + double ld2 = ldx*ldx + ldy*ldy + ldz*ldz; + double rd2 = rdx*rdx + rdy*rdy + rdz*rdz; + if (ld2 < rd2) { + stack[stack_top++] = node->right_child; /* farther pushed first */ + stack[stack_top++] = node->left_child; + } else { + stack[stack_top++] = node->left_child; + stack[stack_top++] = node->right_child; + } } } From 27d23f0ce5eb46dd0cd8be67f00f19020602cf6c Mon Sep 17 00:00:00 2001 From: Luochenghuang Date: Tue, 7 Apr 2026 15:30:27 -0700 Subject: [PATCH 7/8] Add mesh API documentation to ctlgeom.h --- utils/ctlgeom.h | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/utils/ctlgeom.h b/utils/ctlgeom.h index 3ba8fbb..bf9ddb1 100644 --- a/utils/ctlgeom.h +++ b/utils/ctlgeom.h @@ -173,11 +173,14 @@ GEOMETRIC_OBJECT make_slanted_prism_with_center(MATERIAL_TYPE material, vector3 const vector3 *vertices, int num_vertices, double height, vector3 axis, double sidewall_angle); -// mesh with center computed as centroid of vertices +// Closed triangulated 3D mesh. Vertices are shared; triangles index into the +// vertex array (0-based, 3 ints per triangle). A BVH is built internally for +// O(log N) queries. The mesh must be watertight (every edge shared by exactly +// 2 faces); open meshes are detected at init time and produce a warning. GEOMETRIC_OBJECT make_mesh(MATERIAL_TYPE material, const vector3 *vertices, int num_vertices, const int *triangles, int num_triangles); -// mesh with explicit center (vertices shifted so centroid = center) +// As make_mesh, but translates vertices so the centroid equals center. GEOMETRIC_OBJECT make_mesh_with_center(MATERIAL_TYPE material, vector3 center, const vector3 *vertices, int num_vertices, const int *triangles, int num_triangles); From b14367f96b83c5581bc81fa45b90c0e63b571ce7 Mon Sep 17 00:00:00 2001 From: Luochenghuang Date: Tue, 7 Apr 2026 15:48:05 -0700 Subject: [PATCH 8/8] Include generated ctlgeom-types.h and geom-ctl-io.c These files are normally auto-generated from geom.scm via gen-ctl-io (requires Guile). Including them allows building libctlgeom without Guile, which is the common case for meep's Python-only builds. --- utils/ctlgeom-types.h | 209 ++++++++++++++ utils/geom-ctl-io.c | 656 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 865 insertions(+) create mode 100644 utils/ctlgeom-types.h create mode 100644 utils/geom-ctl-io.c diff --git a/utils/ctlgeom-types.h b/utils/ctlgeom-types.h new file mode 100644 index 0000000..ec911f6 --- /dev/null +++ b/utils/ctlgeom-types.h @@ -0,0 +1,209 @@ +/* THIS FILE WAS AUTOMATICALLY GENERATED. DO NOT MODIFY! */ +/* generated from the file: ./geom.scm */ + +#ifndef CTL_IO_H +#define CTL_IO_H + +#include + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + + /******* Type declarations *******/ + + typedef struct geometric_object_struct { + void* material; + vector3 center; + enum { + GEOMETRIC_OBJECT_SELF, PRISM, BLOCK, SPHERE, CYLINDER, COMPOUND_GEOMETRIC_OBJECT, MESH + } which_subclass; + union { + struct prism_struct *prism_data; + struct block_struct *block_data; + struct sphere_struct *sphere_data; + struct cylinder_struct *cylinder_data; + struct compound_geometric_object_struct *compound_geometric_object_data; + struct mesh_struct *mesh_data; + } subclass; + } geometric_object; + + typedef struct { + int num_items; + geometric_object *items; + } geometric_object_list; + + typedef struct compound_geometric_object_struct { + geometric_object_list component_objects; + } compound_geometric_object; + + typedef struct cylinder_struct { + vector3 axis; + number radius; + number height; + enum { + CYLINDER_SELF, WEDGE, CONE + } which_subclass; + union { + struct wedge_struct *wedge_data; + struct cone_struct *cone_data; + } subclass; + } cylinder; + + typedef struct cone_struct { + number radius2; + } cone; + + typedef struct wedge_struct { + number wedge_angle; + vector3 wedge_start; + vector3 e1; + vector3 e2; + } wedge; + + typedef struct sphere_struct { + number radius; + } sphere; + + typedef struct block_struct { + vector3 e1; + vector3 e2; + vector3 e3; + vector3 size; + matrix3x3 projection_matrix; + enum { + BLOCK_SELF, ELLIPSOID + } which_subclass; + union { + struct ellipsoid_struct *ellipsoid_data; + } subclass; + } block; + + typedef struct { + int num_items; + vector3 *items; + } vector3_list; + + typedef struct { + int num_items; + number *items; + } number_list; + + typedef struct prism_struct { + vector3_list vertices; + number height; + vector3 axis; + number sidewall_angle; + vector3_list vertices_p; + vector3_list top_polygon_diff_vectors_p; + vector3_list top_polygon_diff_vectors_scaled_p; + vector3_list vertices_top_p; + vector3_list vertices_top; + vector3 centroid; + number_list workspace; + matrix3x3 m_c2p; + matrix3x3 m_p2c; + } prism; + + typedef struct mesh_bvh_node { + vector3 bbox_low; + vector3 bbox_high; + int left_child; + int right_child; + int face_start; + int face_count; + } mesh_bvh_node; + + typedef struct mesh_struct { + vector3_list vertices; + int num_faces; + int *face_indices; + vector3 *face_normals; + number *face_areas; + int num_bvh_nodes; + mesh_bvh_node *bvh; + int *bvh_face_ids; + boolean is_closed; + vector3 centroid; + } mesh; + + typedef struct ellipsoid_struct { + vector3 inverse_semi_axes; + } ellipsoid; + + typedef struct lattice_struct { + vector3 basis1; + vector3 basis2; + vector3 basis3; + vector3 size; + vector3 basis_size; + vector3 b1; + vector3 b2; + vector3 b3; + matrix3x3 basis; + matrix3x3 metric; + } lattice; + + /******* Input variables *******/ + extern integer dimensions; + extern void* default_material; + extern vector3 geometry_center; + extern lattice geometry_lattice; + extern geometric_object_list geometry; + extern boolean ensure_periodicity; + + /******* Output variables *******/ + + /******* class copy function prototypes *******/ + + extern void lattice_copy(const lattice * o0, lattice * o); + extern void ellipsoid_copy(const ellipsoid * o0, ellipsoid * o); + extern void mesh_copy(const mesh * o0, mesh * o); + extern void prism_copy(const prism * o0, prism * o); + extern void block_copy(const block * o0, block * o); + extern void sphere_copy(const sphere * o0, sphere * o); + extern void wedge_copy(const wedge * o0, wedge * o); + extern void cone_copy(const cone * o0, cone * o); + extern void cylinder_copy(const cylinder * o0, cylinder * o); + extern void compound_geometric_object_copy(const compound_geometric_object * o0, compound_geometric_object * o); + extern void geometric_object_copy(const geometric_object * o0, geometric_object * o); + + /******* class equal function prototypes *******/ + + extern boolean lattice_equal(const lattice * o0, const lattice * o); + extern boolean ellipsoid_equal(const ellipsoid * o0, const ellipsoid * o); + extern boolean mesh_equal(const mesh * o0, const mesh * o); + extern boolean prism_equal(const prism * o0, const prism * o); + extern boolean block_equal(const block * o0, const block * o); + extern boolean sphere_equal(const sphere * o0, const sphere * o); + extern boolean wedge_equal(const wedge * o0, const wedge * o); + extern boolean cone_equal(const cone * o0, const cone * o); + extern boolean cylinder_equal(const cylinder * o0, const cylinder * o); + extern boolean compound_geometric_object_equal(const compound_geometric_object * o0, const compound_geometric_object * o); + extern boolean geometric_object_equal(const geometric_object * o0, const geometric_object * o); + + /******* class destruction function prototypes *******/ + + extern void lattice_destroy(lattice o); + extern void ellipsoid_destroy(ellipsoid o); + extern void mesh_destroy(mesh o); + extern void prism_destroy(prism o); + extern void block_destroy(block o); + extern void sphere_destroy(sphere o); + extern void wedge_destroy(wedge o); + extern void cone_destroy(cone o); + extern void cylinder_destroy(cylinder o); + extern void compound_geometric_object_destroy(compound_geometric_object o); + extern void geometric_object_destroy(geometric_object o); + + +#ifdef __cplusplus +} /* extern "C" */ +#endif /* __cplusplus */ + +#endif /* CTL_IO_H */ +#ifndef LIBCTL_MAJOR_VERSION +# define LIBCTL_MAJOR_VERSION 4 +# define LIBCTL_MINOR_VERSION 5 +# define LIBCTL_BUGFIX_VERSION 1 +#endif diff --git a/utils/geom-ctl-io.c b/utils/geom-ctl-io.c new file mode 100644 index 0000000..434f2fc --- /dev/null +++ b/utils/geom-ctl-io.c @@ -0,0 +1,656 @@ +/* THIS FILE WAS AUTOMATICALLY GENERATED. DO NOT MODIFY! */ +/* generated from the file: ./geom.scm */ + +#include +#include +#include +#include "ctlgeom-types.h" + +#ifdef CXX_CTL_IO +using namespace ctlio; +#endif + + +#include "geom-ctl-io-defaults.c" +#if 0 + +integer dimensions; +SCM default_material; +vector3 geometry_center; +lattice geometry_lattice; +geometric_object_list geometry; +boolean ensure_periodicity; + +#endif + + +/******* class copy functions *******/ + +void +lattice_copy(const lattice * o0, lattice * o) +{ + o->basis1 = o0->basis1; + o->basis2 = o0->basis2; + o->basis3 = o0->basis3; + o->size = o0->size; + o->basis_size = o0->basis_size; + o->b1 = o0->b1; + o->b2 = o0->b2; + o->b3 = o0->b3; + o->basis = o0->basis; + o->metric = o0->metric; +} + +void +ellipsoid_copy(const ellipsoid * o0, ellipsoid * o) +{ + o->inverse_semi_axes = o0->inverse_semi_axes; +} + +void +mesh_copy(const mesh * o0, mesh * o) +{ + { + int i_t; + o->vertices.num_items = o0->vertices.num_items; + o->vertices.items = ((vector3 *) malloc(sizeof(vector3) * (o->vertices.num_items))); + for (i_t = 0; i_t < o->vertices.num_items; i_t++) { + o->vertices.items[i_t] = o0->vertices.items[i_t]; + } + } + o->num_faces = o0->num_faces; + o->face_indices = (int *) malloc(sizeof(int) * 3 * o->num_faces); + memcpy(o->face_indices, o0->face_indices, sizeof(int) * 3 * o->num_faces); + o->face_normals = (vector3 *) malloc(sizeof(vector3) * o->num_faces); + memcpy(o->face_normals, o0->face_normals, sizeof(vector3) * o->num_faces); + o->face_areas = (number *) malloc(sizeof(number) * o->num_faces); + memcpy(o->face_areas, o0->face_areas, sizeof(number) * o->num_faces); + o->num_bvh_nodes = o0->num_bvh_nodes; + o->bvh = (mesh_bvh_node *) malloc(sizeof(mesh_bvh_node) * o->num_bvh_nodes); + memcpy(o->bvh, o0->bvh, sizeof(mesh_bvh_node) * o->num_bvh_nodes); + o->bvh_face_ids = (int *) malloc(sizeof(int) * o->num_faces); + memcpy(o->bvh_face_ids, o0->bvh_face_ids, sizeof(int) * o->num_faces); + o->is_closed = o0->is_closed; + o->centroid = o0->centroid; +} + +void +prism_copy(const prism * o0, prism * o) +{ + { + int i_t; + o->vertices.num_items = o0->vertices.num_items; + o->vertices.items = ((vector3 *) malloc(sizeof(vector3) * (o->vertices.num_items))); + for (i_t = 0; i_t < o->vertices.num_items; i_t++) { + o->vertices.items[i_t] = o0->vertices.items[i_t]; + } + } + o->height = o0->height; + o->axis = o0->axis; + o->sidewall_angle = o0->sidewall_angle; + { + int i_t; + o->vertices_p.num_items = o0->vertices_p.num_items; + o->vertices_p.items = ((vector3 *) malloc(sizeof(vector3) * (o->vertices_p.num_items))); + for (i_t = 0; i_t < o->vertices_p.num_items; i_t++) { + o->vertices_p.items[i_t] = o0->vertices_p.items[i_t]; + } + } + { + int i_t; + o->top_polygon_diff_vectors_p.num_items = o0->top_polygon_diff_vectors_p.num_items; + o->top_polygon_diff_vectors_p.items = ((vector3 *) malloc(sizeof(vector3) * (o->top_polygon_diff_vectors_p.num_items))); + for (i_t = 0; i_t < o->top_polygon_diff_vectors_p.num_items; i_t++) { + o->top_polygon_diff_vectors_p.items[i_t] = o0->top_polygon_diff_vectors_p.items[i_t]; + } + } + { + int i_t; + o->top_polygon_diff_vectors_scaled_p.num_items = o0->top_polygon_diff_vectors_scaled_p.num_items; + o->top_polygon_diff_vectors_scaled_p.items = ((vector3 *) malloc(sizeof(vector3) * (o->top_polygon_diff_vectors_scaled_p.num_items))); + for (i_t = 0; i_t < o->top_polygon_diff_vectors_scaled_p.num_items; i_t++) { + o->top_polygon_diff_vectors_scaled_p.items[i_t] = o0->top_polygon_diff_vectors_scaled_p.items[i_t]; + } + } + { + int i_t; + o->vertices_top_p.num_items = o0->vertices_top_p.num_items; + o->vertices_top_p.items = ((vector3 *) malloc(sizeof(vector3) * (o->vertices_top_p.num_items))); + for (i_t = 0; i_t < o->vertices_top_p.num_items; i_t++) { + o->vertices_top_p.items[i_t] = o0->vertices_top_p.items[i_t]; + } + } + { + int i_t; + o->vertices_top.num_items = o0->vertices_top.num_items; + o->vertices_top.items = ((vector3 *) malloc(sizeof(vector3) * (o->vertices_top.num_items))); + for (i_t = 0; i_t < o->vertices_top.num_items; i_t++) { + o->vertices_top.items[i_t] = o0->vertices_top.items[i_t]; + } + } + o->centroid = o0->centroid; + { + int i_t; + o->workspace.num_items = o0->workspace.num_items; + o->workspace.items = ((number *) malloc(sizeof(number) * (o->workspace.num_items))); + for (i_t = 0; i_t < o->workspace.num_items; i_t++) { + o->workspace.items[i_t] = o0->workspace.items[i_t]; + } + } + o->m_c2p = o0->m_c2p; + o->m_p2c = o0->m_p2c; +} + +void +block_copy(const block * o0, block * o) +{ + o->e1 = o0->e1; + o->e2 = o0->e2; + o->e3 = o0->e3; + o->size = o0->size; + o->projection_matrix = o0->projection_matrix; + if (o0->which_subclass == ELLIPSOID) { + o->which_subclass = ELLIPSOID; + o->subclass.ellipsoid_data = ((ellipsoid *) malloc(sizeof(ellipsoid))); + ellipsoid_copy(o0->subclass.ellipsoid_data, o->subclass.ellipsoid_data); + } else + o->which_subclass = BLOCK_SELF; +} + +void +sphere_copy(const sphere * o0, sphere * o) +{ + o->radius = o0->radius; +} + +void +wedge_copy(const wedge * o0, wedge * o) +{ + o->wedge_angle = o0->wedge_angle; + o->wedge_start = o0->wedge_start; + o->e1 = o0->e1; + o->e2 = o0->e2; +} + +void +cone_copy(const cone * o0, cone * o) +{ + o->radius2 = o0->radius2; +} + +void +cylinder_copy(const cylinder * o0, cylinder * o) +{ + o->axis = o0->axis; + o->radius = o0->radius; + o->height = o0->height; + if (o0->which_subclass == WEDGE) { + o->which_subclass = WEDGE; + o->subclass.wedge_data = ((wedge *) malloc(sizeof(wedge))); + wedge_copy(o0->subclass.wedge_data, o->subclass.wedge_data); + } else if (o0->which_subclass == CONE) { + o->which_subclass = CONE; + o->subclass.cone_data = ((cone *) malloc(sizeof(cone))); + cone_copy(o0->subclass.cone_data, o->subclass.cone_data); + } else + o->which_subclass = CYLINDER_SELF; +} + +void +compound_geometric_object_copy(const compound_geometric_object * o0, compound_geometric_object * o) +{ + { + int i_t; + o->component_objects.num_items = o0->component_objects.num_items; + o->component_objects.items = ((geometric_object *) malloc(sizeof(geometric_object) * (o->component_objects.num_items))); + for (i_t = 0; i_t < o->component_objects.num_items; i_t++) { + geometric_object_copy(&o0->component_objects.items[i_t], &o->component_objects.items[i_t]); + } + } +} + +void +geometric_object_copy(const geometric_object * o0, geometric_object * o) +{ + o->material = o0->material; + o->center = o0->center; + if (o0->which_subclass == PRISM) { + o->which_subclass = PRISM; + o->subclass.prism_data = ((prism *) malloc(sizeof(prism))); + prism_copy(o0->subclass.prism_data, o->subclass.prism_data); + } else if (o0->which_subclass == BLOCK) { + o->which_subclass = BLOCK; + o->subclass.block_data = ((block *) malloc(sizeof(block))); + block_copy(o0->subclass.block_data, o->subclass.block_data); + } else if (o0->which_subclass == SPHERE) { + o->which_subclass = SPHERE; + o->subclass.sphere_data = ((sphere *) malloc(sizeof(sphere))); + sphere_copy(o0->subclass.sphere_data, o->subclass.sphere_data); + } else if (o0->which_subclass == CYLINDER) { + o->which_subclass = CYLINDER; + o->subclass.cylinder_data = ((cylinder *) malloc(sizeof(cylinder))); + cylinder_copy(o0->subclass.cylinder_data, o->subclass.cylinder_data); + } else if (o0->which_subclass == COMPOUND_GEOMETRIC_OBJECT) { + o->which_subclass = COMPOUND_GEOMETRIC_OBJECT; + o->subclass.compound_geometric_object_data = ((compound_geometric_object *) malloc(sizeof(compound_geometric_object))); + compound_geometric_object_copy(o0->subclass.compound_geometric_object_data, o->subclass.compound_geometric_object_data); + } else if (o0->which_subclass == MESH) { + o->which_subclass = MESH; + o->subclass.mesh_data = ((mesh *) malloc(sizeof(mesh))); + mesh_copy(o0->subclass.mesh_data, o->subclass.mesh_data); + } else + o->which_subclass = GEOMETRIC_OBJECT_SELF; +} + +/******* class equal functions *******/ + +boolean +lattice_equal(const lattice * o0, const lattice * o) +{ + if (!vector3_equal(o->basis1, o0->basis1)) + return 0; + if (!vector3_equal(o->basis2, o0->basis2)) + return 0; + if (!vector3_equal(o->basis3, o0->basis3)) + return 0; + if (!vector3_equal(o->size, o0->size)) + return 0; + if (!vector3_equal(o->basis_size, o0->basis_size)) + return 0; + if (!vector3_equal(o->b1, o0->b1)) + return 0; + if (!vector3_equal(o->b2, o0->b2)) + return 0; + if (!vector3_equal(o->b3, o0->b3)) + return 0; + if (!matrix3x3_equal(o->basis, o0->basis)) + return 0; + if (!matrix3x3_equal(o->metric, o0->metric)) + return 0; + ; + return 1; +} + +boolean +ellipsoid_equal(const ellipsoid * o0, const ellipsoid * o) +{ + if (!vector3_equal(o->inverse_semi_axes, o0->inverse_semi_axes)) + return 0; + ; + return 1; +} + +boolean +mesh_equal(const mesh * o0, const mesh * o) +{ + { + int i_t; + if (o->vertices.num_items != o0->vertices.num_items) + return 0; + for (i_t = 0; i_t < o->vertices.num_items; i_t++) { + if (!vector3_equal(o->vertices.items[i_t], o0->vertices.items[i_t])) + return 0; + } + } + if (o->num_faces != o0->num_faces) + return 0; + { + int i_t; + for (i_t = 0; i_t < 3 * o->num_faces; i_t++) { + if (o->face_indices[i_t] != o0->face_indices[i_t]) + return 0; + } + } + ; + return 1; +} + +boolean +prism_equal(const prism * o0, const prism * o) +{ + { + int i_t; + if (o->vertices.num_items != o0->vertices.num_items) + return 0; + for (i_t = 0; i_t < o->vertices.num_items; i_t++) { + if (!vector3_equal(o->vertices.items[i_t], o0->vertices.items[i_t])) + return 0; + } + } + if (o->height != o0->height) + return 0; + if (!vector3_equal(o->axis, o0->axis)) + return 0; + if (o->sidewall_angle != o0->sidewall_angle) + return 0; + { + int i_t; + if (o->vertices_p.num_items != o0->vertices_p.num_items) + return 0; + for (i_t = 0; i_t < o->vertices_p.num_items; i_t++) { + if (!vector3_equal(o->vertices_p.items[i_t], o0->vertices_p.items[i_t])) + return 0; + } + } + { + int i_t; + if (o->top_polygon_diff_vectors_p.num_items != o0->top_polygon_diff_vectors_p.num_items) + return 0; + for (i_t = 0; i_t < o->top_polygon_diff_vectors_p.num_items; i_t++) { + if (!vector3_equal(o->top_polygon_diff_vectors_p.items[i_t], o0->top_polygon_diff_vectors_p.items[i_t])) + return 0; + } + } + { + int i_t; + if (o->top_polygon_diff_vectors_scaled_p.num_items != o0->top_polygon_diff_vectors_scaled_p.num_items) + return 0; + for (i_t = 0; i_t < o->top_polygon_diff_vectors_scaled_p.num_items; i_t++) { + if (!vector3_equal(o->top_polygon_diff_vectors_scaled_p.items[i_t], o0->top_polygon_diff_vectors_scaled_p.items[i_t])) + return 0; + } + } + { + int i_t; + if (o->vertices_top_p.num_items != o0->vertices_top_p.num_items) + return 0; + for (i_t = 0; i_t < o->vertices_top_p.num_items; i_t++) { + if (!vector3_equal(o->vertices_top_p.items[i_t], o0->vertices_top_p.items[i_t])) + return 0; + } + } + { + int i_t; + if (o->vertices_top.num_items != o0->vertices_top.num_items) + return 0; + for (i_t = 0; i_t < o->vertices_top.num_items; i_t++) { + if (!vector3_equal(o->vertices_top.items[i_t], o0->vertices_top.items[i_t])) + return 0; + } + } + if (!vector3_equal(o->centroid, o0->centroid)) + return 0; + { + int i_t; + if (o->workspace.num_items != o0->workspace.num_items) + return 0; + for (i_t = 0; i_t < o->workspace.num_items; i_t++) { + if (o->workspace.items[i_t] != o0->workspace.items[i_t]) + return 0; + } + } + if (!matrix3x3_equal(o->m_c2p, o0->m_c2p)) + return 0; + if (!matrix3x3_equal(o->m_p2c, o0->m_p2c)) + return 0; + ; + return 1; +} + +boolean +block_equal(const block * o0, const block * o) +{ + if (!vector3_equal(o->e1, o0->e1)) + return 0; + if (!vector3_equal(o->e2, o0->e2)) + return 0; + if (!vector3_equal(o->e3, o0->e3)) + return 0; + if (!vector3_equal(o->size, o0->size)) + return 0; + if (!matrix3x3_equal(o->projection_matrix, o0->projection_matrix)) + return 0; + if (o0->which_subclass != o->which_subclass) + return 0; + if (o0->which_subclass == ELLIPSOID) { + if (!ellipsoid_equal(o0->subclass.ellipsoid_data, o->subclass.ellipsoid_data)) + return 0; + } else; + return 1; +} + +boolean +sphere_equal(const sphere * o0, const sphere * o) +{ + if (o->radius != o0->radius) + return 0; + ; + return 1; +} + +boolean +wedge_equal(const wedge * o0, const wedge * o) +{ + if (o->wedge_angle != o0->wedge_angle) + return 0; + if (!vector3_equal(o->wedge_start, o0->wedge_start)) + return 0; + if (!vector3_equal(o->e1, o0->e1)) + return 0; + if (!vector3_equal(o->e2, o0->e2)) + return 0; + ; + return 1; +} + +boolean +cone_equal(const cone * o0, const cone * o) +{ + if (o->radius2 != o0->radius2) + return 0; + ; + return 1; +} + +boolean +cylinder_equal(const cylinder * o0, const cylinder * o) +{ + if (!vector3_equal(o->axis, o0->axis)) + return 0; + if (o->radius != o0->radius) + return 0; + if (o->height != o0->height) + return 0; + if (o0->which_subclass != o->which_subclass) + return 0; + if (o0->which_subclass == WEDGE) { + if (!wedge_equal(o0->subclass.wedge_data, o->subclass.wedge_data)) + return 0; + } else if (o0->which_subclass == CONE) { + if (!cone_equal(o0->subclass.cone_data, o->subclass.cone_data)) + return 0; + } else; + return 1; +} + +boolean +compound_geometric_object_equal(const compound_geometric_object * o0, const compound_geometric_object * o) +{ + { + int i_t; + if (o->component_objects.num_items != o0->component_objects.num_items) + return 0; + for (i_t = 0; i_t < o->component_objects.num_items; i_t++) { + if (!geometric_object_equal(&o0->component_objects.items[i_t], &o->component_objects.items[i_t])) + return 0; + } + } + ; + return 1; +} + +boolean +geometric_object_equal(const geometric_object * o0, const geometric_object * o) +{ + if (o->material != o0->material) + return 0; + if (!vector3_equal(o->center, o0->center)) + return 0; + if (o0->which_subclass != o->which_subclass) + return 0; + if (o0->which_subclass == PRISM) { + if (!prism_equal(o0->subclass.prism_data, o->subclass.prism_data)) + return 0; + } else if (o0->which_subclass == BLOCK) { + if (!block_equal(o0->subclass.block_data, o->subclass.block_data)) + return 0; + } else if (o0->which_subclass == SPHERE) { + if (!sphere_equal(o0->subclass.sphere_data, o->subclass.sphere_data)) + return 0; + } else if (o0->which_subclass == CYLINDER) { + if (!cylinder_equal(o0->subclass.cylinder_data, o->subclass.cylinder_data)) + return 0; + } else if (o0->which_subclass == COMPOUND_GEOMETRIC_OBJECT) { + if (!compound_geometric_object_equal(o0->subclass.compound_geometric_object_data, o->subclass.compound_geometric_object_data)) + return 0; + } else if (o0->which_subclass == MESH) { + if (!mesh_equal(o0->subclass.mesh_data, o->subclass.mesh_data)) + return 0; + } else; + return 1; +} + +/******* class destruction functions *******/ + +void +lattice_destroy(lattice o) +{ +} + +void +ellipsoid_destroy(ellipsoid o) +{ +} + +void +mesh_destroy(mesh o) +{ + free(o.vertices.items); + free(o.face_indices); + free(o.face_normals); + free(o.face_areas); + free(o.bvh); + free(o.bvh_face_ids); +} + +void +prism_destroy(prism o) +{ + { + int index_t; + for (index_t = 0; index_t < o.vertices.num_items; index_t++) { + } + } + free(o.vertices.items); + { + int index_t; + for (index_t = 0; index_t < o.vertices_p.num_items; index_t++) { + } + } + free(o.vertices_p.items); + { + int index_t; + for (index_t = 0; index_t < o.top_polygon_diff_vectors_p.num_items; index_t++) { + } + } + free(o.top_polygon_diff_vectors_p.items); + { + int index_t; + for (index_t = 0; index_t < o.top_polygon_diff_vectors_scaled_p.num_items; index_t++) { + } + } + free(o.top_polygon_diff_vectors_scaled_p.items); + { + int index_t; + for (index_t = 0; index_t < o.vertices_top_p.num_items; index_t++) { + } + } + free(o.vertices_top_p.items); + { + int index_t; + for (index_t = 0; index_t < o.vertices_top.num_items; index_t++) { + } + } + free(o.vertices_top.items); + { + int index_t; + for (index_t = 0; index_t < o.workspace.num_items; index_t++) { + } + } + free(o.workspace.items); +} + +void +block_destroy(block o) +{ + if (o.which_subclass == ELLIPSOID) { + ellipsoid_destroy(*o.subclass.ellipsoid_data); + free(o.subclass.ellipsoid_data); + } else { + } +} + +void +sphere_destroy(sphere o) +{ +} + +void +wedge_destroy(wedge o) +{ +} + +void +cone_destroy(cone o) +{ +} + +void +cylinder_destroy(cylinder o) +{ + if (o.which_subclass == WEDGE) { + wedge_destroy(*o.subclass.wedge_data); + free(o.subclass.wedge_data); + } else if (o.which_subclass == CONE) { + cone_destroy(*o.subclass.cone_data); + free(o.subclass.cone_data); + } else { + } +} + +void +compound_geometric_object_destroy(compound_geometric_object o) +{ + { + int index_t; + for (index_t = 0; index_t < o.component_objects.num_items; index_t++) { + geometric_object_destroy(o.component_objects.items[index_t]); + } + } + free(o.component_objects.items); +} + +void +geometric_object_destroy(geometric_object o) +{ + if (o.which_subclass == PRISM) { + prism_destroy(*o.subclass.prism_data); + free(o.subclass.prism_data); + } else if (o.which_subclass == BLOCK) { + block_destroy(*o.subclass.block_data); + free(o.subclass.block_data); + } else if (o.which_subclass == SPHERE) { + sphere_destroy(*o.subclass.sphere_data); + free(o.subclass.sphere_data); + } else if (o.which_subclass == CYLINDER) { + cylinder_destroy(*o.subclass.cylinder_data); + free(o.subclass.cylinder_data); + } else if (o.which_subclass == COMPOUND_GEOMETRIC_OBJECT) { + compound_geometric_object_destroy(*o.subclass.compound_geometric_object_data); + free(o.subclass.compound_geometric_object_data); + } else if (o.which_subclass == MESH) { + mesh_destroy(*o.subclass.mesh_data); + free(o.subclass.mesh_data); + } else { + } +}