From ba1b826f08813e2eaba50d24479cb09535499b7a Mon Sep 17 00:00:00 2001
From: Luochenghuang <Luocheng@uw.edu>
Date: Fri, 29 May 2026 23:22:24 -0700
Subject: [PATCH 1/5] Refactor mesh storage with PIMPL: opaque internal cache
 behind scheme-driven public struct
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

PR #74 (mesh geometry import) hand-extended utils/ctlgeom-types.h with C-only internal fields (mesh_bvh_node, face_normals, face_areas, bvh, bvh_face_ids, is_closed, centroid, lengthscale, num_bvh_nodes, and a flat int* face_indices) because none of those had natural Scheme representations in geom.scm. That worked when Guile wasn't detected at configure time, but any from-git build with Guile installed (including the new GitHub Actions CI) regenerated ctlgeom-types.h from geom.scm, clobbered the hand-extensions, and failed to compile with 'unknown type name mesh_bvh_node'.

This commit applies the PIMPL (pointer-to-implementation) idiom to clean up that situation properly. The public mesh struct in ctlgeom-types.h now contains only Scheme-visible fields — vertices, face_indices (as a vector3_list, each vector3 packing one triangle's three int indices), and a void* internal pointer. All C-only state moves into mesh_internal in a new private header utils/mesh-internal.h, accessed via a static-inline mesh_priv(m) cast helper. ctlgeom-types.h goes back to being a regular auto-generated artifact derived from geom.scm; no Makefile.am workaround needed.

Lifecycle:
  - make_mesh / make_mesh_with_center: populate public fields, then call init_mesh which calloc()s mesh_internal, unpacks face_indices into a flat int[], computes face normals/areas, checks watertight closure, builds the BVH, and stores the pointer in m->internal.
  - reinit_mesh: fast-path returns when m->internal != NULL; otherwise calls init_mesh. Same single-threaded-init contract as before.
  - mesh_destroy (auto-generated): frees vertices.items + face_indices.items, then calls mesh_internal_free (defined in geom.c) to free the cache.
  - mesh_copy (auto-generated): deep-copies the public fields, sets internal=NULL, then calls mesh_init_internal (defined in geom.c) to eagerly rebuild the BVH on the copy. Eliminates the previous shallow-pointer-share footgun where two mesh copies pointed at the same BVH allocation.

geom.c gets ~105 mechanical m->X → mesh_priv(m)->X rewrites for the 10 internal fields (the cast inlines to nothing at runtime). test-mesh.c gets the same rewrite for 6 sites where tests poke directly at m->is_closed.

Both auto-generated files (ctlgeom-types.h and geom-ctl-io.c) now match what gen-ctl-io would emit from the simplified geom.scm, with one hand-tweak in geom-ctl-io.c: mesh_destroy and mesh_copy call into the externally-defined mesh_internal_free / mesh_init_internal so they can manage the opaque cache. This is the minimum manual surface required for any PIMPL design.

Verified: configure + make + make check all clean. test-prism passes (0/10000 point failures, 0/1000 segment failures, 0/65 prism helper failures). test-mesh passes 0 failures across all 19 cases including test_copy_destroy (which exercises mesh_copy + point queries on the copy + destroy of both original and copy).
---
 utils/ctlgeom-types.h |  21 +---
 utils/geom-ctl-io.c   |  46 ++++---
 utils/geom.c          | 285 ++++++++++++++++++++++++------------------
 utils/geom.scm        |   6 +-
 utils/mesh-internal.h |  45 +++++++
 utils/test-mesh.c     |  13 +-
 6 files changed, 245 insertions(+), 171 deletions(-)
 create mode 100644 utils/mesh-internal.h

diff --git a/utils/ctlgeom-types.h b/utils/ctlgeom-types.h
index 02744d8..f92be6d 100644
--- a/utils/ctlgeom-types.h
+++ b/utils/ctlgeom-types.h
@@ -105,27 +105,10 @@ extern		"C" {
 		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;
-		number		lengthscale;
+		vector3_list	face_indices;
+		void*		internal;
 	}		mesh;
 
 	typedef struct ellipsoid_struct {
diff --git a/utils/geom-ctl-io.c b/utils/geom-ctl-io.c
index 88d694c..9e61e6d 100644
--- a/utils/geom-ctl-io.c
+++ b/utils/geom-ctl-io.c
@@ -47,6 +47,9 @@ ellipsoid_copy(const ellipsoid * o0, ellipsoid * o)
 	o->inverse_semi_axes = o0->inverse_semi_axes;
 }
 
+/* Defined in geom.c; eagerly rebuilds the opaque mesh_internal cache. */
+extern void mesh_init_internal(mesh *m);
+
 void
 mesh_copy(const mesh * o0, mesh * o)
 {
@@ -58,21 +61,16 @@ mesh_copy(const mesh * o0, mesh * o)
 			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;
-	o->lengthscale = o0->lengthscale;
+	{
+		int		i_t;
+		o->face_indices.num_items = o0->face_indices.num_items;
+		o->face_indices.items = ((vector3 *) malloc(sizeof(vector3) * (o->face_indices.num_items)));
+		for (i_t = 0; i_t < o->face_indices.num_items; i_t++) {
+			o->face_indices.items[i_t] = o0->face_indices.items[i_t];
+		}
+	}
+	o->internal = NULL;
+	mesh_init_internal(o);  /* rebuild BVH eagerly; safe under later concurrent queries */
 }
 
 void
@@ -293,12 +291,12 @@ mesh_equal(const mesh * o0, const mesh * o)
 				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])
+		if (o->face_indices.num_items != o0->face_indices.num_items)
+			return 0;
+		for (i_t = 0; i_t < o->face_indices.num_items; i_t++) {
+			if (!vector3_equal(o->face_indices.items[i_t], o0->face_indices.items[i_t]))
 				return 0;
 		}
 	}
@@ -523,15 +521,15 @@ ellipsoid_destroy(ellipsoid o)
 {
 }
 
+/* Defined in geom.c; frees the opaque mesh_internal cache. Safe on NULL. */
+extern void mesh_internal_free(void *p);
+
 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);
+	free(o.face_indices.items);
+	mesh_internal_free(o.internal);
 }
 
 void
diff --git a/utils/geom.c b/utils/geom.c
index 1020b1e..0f08949 100644
--- a/utils/geom.c
+++ b/utils/geom.c
@@ -33,6 +33,7 @@
 static void material_type_copy(void **src, void **dest) { *dest = *src; }
 #endif
 #include "ctlgeom.h"
+#include "mesh-internal.h"
 
 #ifdef CXX_CTL_IO
 using namespace ctlio;
@@ -1969,9 +1970,9 @@ static void bvh_node_set_box(mesh_bvh_node *node, const geom_box *box) {
 /* 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]];
+  *v0 = m->vertices.items[mesh_priv(m)->face_indices[3 * face_id]];
+  *v1 = m->vertices.items[mesh_priv(m)->face_indices[3 * face_id + 1]];
+  *v2 = m->vertices.items[mesh_priv(m)->face_indices[3 * face_id + 2]];
 }
 
 /* Compute the AABB of a single triangle. */
@@ -2058,7 +2059,7 @@ static int mesh_bvh_build(mesh *m, int *face_ids, int start, int count,
   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 * m->lengthscale * m->lengthscale;
+  if (parent_area == 0) parent_area = 1e-30 * mesh_priv(m)->lengthscale * mesh_priv(m)->lengthscale;
 
   for (int axis = 0; axis < 3; axis++) {
     double lo, hi;
@@ -2067,7 +2068,7 @@ static int mesh_bvh_build(mesh *m, int *face_ids, int start, int count,
       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 * m->lengthscale) continue;
+    if (hi - lo < 1e-15 * mesh_priv(m)->lengthscale) continue;
 
     /* Bin face centroids. */
     int bin_counts[MESH_BVH_NUM_BINS];
@@ -2364,7 +2365,7 @@ static int find_closest_face(const mesh *m, vector3 p, double *dist2) {
 
   while (stack_top > 0) {
     int node_idx = stack[--stack_top];
-    const mesh_bvh_node *node = &m->bvh[node_idx];
+    const mesh_bvh_node *node = &mesh_priv(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));
@@ -2376,10 +2377,10 @@ static int find_closest_face(const mesh *m, vector3 p, double *dist2) {
     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]];
+        int fid = mesh_priv(m)->bvh_face_ids[node->face_start + i];
+        vector3 v0 = m->vertices.items[mesh_priv(m)->face_indices[3 * fid]];
+        vector3 v1 = m->vertices.items[mesh_priv(m)->face_indices[3 * fid + 1]];
+        vector3 v2 = m->vertices.items[mesh_priv(m)->face_indices[3 * fid + 2]];
         vector3 closest;
         double d2 = closest_point_on_triangle(p, v0, v1, v2, &closest);
         if (d2 < best_dist2) {
@@ -2391,8 +2392,8 @@ static int find_closest_face(const mesh *m, vector3 p, double *dist2) {
       if (stack_top + 2 > 64) continue;
       /* 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];
+      const mesh_bvh_node *left = &mesh_priv(m)->bvh[node->left_child];
+      const mesh_bvh_node *right = &mesh_priv(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));
@@ -2462,24 +2463,24 @@ static void mesh_ray_all_intersections(const mesh *m, vector3 origin, vector3 di
   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;
 
-  double det_eps = 1e-12 * m->lengthscale * m->lengthscale;
+  double det_eps = 1e-12 * mesh_priv(m)->lengthscale * mesh_priv(m)->lengthscale;
   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];
+    const mesh_bvh_node *node = &mesh_priv(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]];
+        int fid = mesh_priv(m)->bvh_face_ids[node->face_start + i];
+        vector3 v0 = m->vertices.items[mesh_priv(m)->face_indices[3 * fid]];
+        vector3 v1 = m->vertices.items[mesh_priv(m)->face_indices[3 * fid + 1]];
+        vector3 v2 = m->vertices.items[mesh_priv(m)->face_indices[3 * fid + 2]];
 
         double t;
         if (ray_triangle_intersect(origin, dir, v0, v1, v2, det_eps, &t, NULL, NULL))
@@ -2523,7 +2524,7 @@ static int count_ray_mesh_intersections_ex(const mesh *m, vector3 origin, vector
   mesh_ray_all_intersections(m, origin, dir, &hits);
 
   /* Keep only forward intersections (t > eps). */
-  double fwd_eps = 1e-12 * m->lengthscale;
+  double fwd_eps = 1e-12 * mesh_priv(m)->lengthscale;
   int nforward = 0;
   for (int i = 0; i < hits.count; i++)
     if (hits.data[i] > fwd_eps)
@@ -2531,7 +2532,7 @@ static int count_ray_mesh_intersections_ex(const mesh *m, vector3 origin, vector
 
   if (nforward > 1) {
     qsort(hits.data, nforward, sizeof(double), mesh_dcmp);
-    double dedup_tol = 1e-10 * m->lengthscale;
+    double dedup_tol = 1e-10 * mesh_priv(m)->lengthscale;
     int deduped = remove_duplicate_intersections(hits.data, nforward, dedup_tol);
     if (had_degenerate) *had_degenerate = (deduped != nforward);
     nforward = deduped;
@@ -2558,7 +2559,7 @@ static const vector3 mesh_ray_dirs[3] = {
 };
 
 static boolean point_in_mesh(const mesh *m, vector3 p) {
-  if (!m->is_closed) return 0;
+  if (!mesh_priv(m)->is_closed) return 0;
 
   /* Fast path: cast one ray. If no degenerate edge/vertex hits were
      detected during deduplication, trust the result immediately.
@@ -2585,13 +2586,13 @@ static vector3 normal_to_mesh(const mesh *m, vector3 p) {
     vector3 zero = {0, 0, 0};
     return zero;
   }
-  return m->face_normals[face];
+  return mesh_priv(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;
+  if (mesh_priv(m)->num_bvh_nodes > 0) {
+    box->low = mesh_priv(m)->bvh[0].bbox_low;
+    box->high = mesh_priv(m)->bvh[0].bbox_high;
   } else {
     box->low = box->high = m->vertices.items[0];
     for (int i = 1; i < m->vertices.num_items; i++)
@@ -2602,10 +2603,10 @@ static void get_mesh_bounding_box(const mesh *m, geom_box *box) {
 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]];
+  for (int f = 0; f < mesh_priv(m)->num_faces; f++) {
+    vector3 v0 = m->vertices.items[mesh_priv(m)->face_indices[3 * f]];
+    vector3 v1 = m->vertices.items[mesh_priv(m)->face_indices[3 * f + 1]];
+    vector3 v2 = m->vertices.items[mesh_priv(m)->face_indices[3 * f + 2]];
     vol += vector3_dot(v0, vector3_cross(v1, v2));
   }
   return fabs(vol) / 6.0;
@@ -2614,8 +2615,8 @@ static double get_mesh_volume(const mesh *m) {
 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)");
+             m->vertices.num_items, mesh_priv(m)->num_faces,
+             mesh_priv(m)->is_closed ? "closed" : "OPEN (WARNING)");
 }
 
 static double intersect_line_segment_with_mesh(const mesh *m, vector3 p, vector3 d,
@@ -2626,7 +2627,7 @@ static double intersect_line_segment_with_mesh(const mesh *m, vector3 p, vector3
 
   if (hits.count > 1) {
     qsort(hits.data, hits.count, sizeof(double), mesh_dcmp);
-    hits.count = remove_duplicate_intersections(hits.data, hits.count, 1e-10 * m->lengthscale);
+    hits.count = remove_duplicate_intersections(hits.data, hits.count, 1e-10 * mesh_priv(m)->lengthscale);
   }
 
   /* The sorted intersection list gives all surface crossings along the
@@ -2653,30 +2654,75 @@ static double intersect_line_segment_with_mesh(const mesh *m, vector3 p, vector3
   return ds > 0.0 ? ds : 0.0;
 }
 
+/* Forward declaration; init_mesh body lives below. */
+static void init_mesh(geometric_object *o);
+
+/* Build the opaque mesh_internal cache for a mesh whose public fields
+   (vertices, face_indices) have just been populated and whose internal
+   pointer is NULL. Used by the auto-generated mesh_copy in geom-ctl-io.c
+   to eagerly rebuild a copied mesh's BVH so subsequent _fixed_ queries
+   (which skip geom_fix_object_ptr / reinit_mesh) work correctly. */
+void mesh_init_internal(mesh *m) {
+  geometric_object o;
+  o.subclass.mesh_data = m;
+  init_mesh(&o);
+}
+
+/* Free the opaque mesh_internal cache and all its nested allocations.
+   Safe on NULL. Called from the auto-generated mesh_destroy in
+   geom-ctl-io.c via the extern declaration there. */
+void mesh_internal_free(void *p) {
+  if (!p) return;
+  mesh_internal *mi = (mesh_internal *)p;
+  free(mi->face_indices);
+  free(mi->face_normals);
+  free(mi->face_areas);
+  free(mi->bvh);
+  free(mi->bvh_face_ids);
+  free(mi);
+}
+
 /***************************************************************/
-/* init_mesh: validate, compute normals, build BVH             */
+/* init_mesh: allocate the opaque mesh_internal cache, unpack   */
+/* face_indices into a flat int array, compute face normals,    */
+/* areas, centroid, lengthscale, check closure, and build BVH.  */
 /*                                                              */
-/* NOT THREAD-SAFE: mutates m->face_indices (winding flip),     */
-/* allocates m->face_normals / m->face_areas / m->bvh /         */
-/* m->bvh_face_ids, and writes m->is_closed, m->centroid,       */
-/* m->lengthscale, m->num_bvh_nodes. Only invoked from the      */
-/* mesh constructors (make_mesh / make_mesh_with_center) and    */
-/* from reinit_mesh; both rely on it running single-threaded    */
-/* per mesh.                                                    */
+/* NOT THREAD-SAFE: writes the per-mesh internal cache.         */
+/* Invoked only from the mesh constructors and from reinit_mesh;*/
+/* both rely on it running single-threaded per mesh.            */
 /***************************************************************/
 
 static void init_mesh(geometric_object *o) {
   mesh *m = o->subclass.mesh_data;
   int nv = m->vertices.num_items;
-  int nf = m->num_faces;
+
+  /* Allocate the opaque internal cache. m->internal must be NULL on entry
+     (either the mesh was just constructed, or reinit_mesh cleared it). */
+  mesh_internal *p = (mesh_internal *)calloc(1, sizeof(mesh_internal));
+  CHECK(p, "out of memory");
+  m->internal = p;
+
+  /* Unpack face_indices: the public vector3_list stores 3 ints per
+     triangle packed into a vector3 (x, y, z are the indices as doubles,
+     exactly representable up to 2^53). The internal cache uses a flat
+     int[3*num_faces] for fast indexed access. */
+  int nf = m->face_indices.num_items;
+  p->num_faces = nf;
+  p->face_indices = (int *)malloc(3 * nf * sizeof(int));
+  CHECK(p->face_indices, "out of memory");
+  for (int f = 0; f < nf; f++) {
+    p->face_indices[3 * f]     = (int)m->face_indices.items[f].x;
+    p->face_indices[3 * f + 1] = (int)m->face_indices.items[f].y;
+    p->face_indices[3 * f + 2] = (int)m->face_indices.items[f].z;
+  }
 
   /* 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");
+    CHECK(mesh_priv(m)->face_indices[3 * f] >= 0 && mesh_priv(m)->face_indices[3 * f] < nv, "mesh face index out of range");
+    CHECK(mesh_priv(m)->face_indices[3 * f + 1] >= 0 && mesh_priv(m)->face_indices[3 * f + 1] < nv, "mesh face index out of range");
+    CHECK(mesh_priv(m)->face_indices[3 * f + 2] >= 0 && mesh_priv(m)->face_indices[3 * f + 2] < nv, "mesh face index out of range");
   }
 
   /* Compute characteristic lengthscale from vertex bounding box diagonal. */
@@ -2688,34 +2734,34 @@ static void init_mesh(geometric_object *o) {
       hi.x = fmax(hi.x, v.x); hi.y = fmax(hi.y, v.y); hi.z = fmax(hi.z, v.z);
     }
     double dx = hi.x - lo.x, dy = hi.y - lo.y, dz = hi.z - lo.z;
-    m->lengthscale = sqrt(dx * dx + dy * dy + dz * dz);
-    if (m->lengthscale == 0) m->lengthscale = 1.0;
+    mesh_priv(m)->lengthscale = sqrt(dx * dx + dy * dy + dz * dz);
+    if (mesh_priv(m)->lengthscale == 0) mesh_priv(m)->lengthscale = 1.0;
   }
 
   /* 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");
+  mesh_priv(m)->face_normals = (vector3 *)malloc(nf * sizeof(vector3));
+  CHECK(mesh_priv(m)->face_normals, "out of memory");
+  mesh_priv(m)->face_areas = (double *)malloc(nf * sizeof(double));
+  CHECK(mesh_priv(m)->face_areas, "out of memory");
 
-  double area_eps = 1e-20 * m->lengthscale * m->lengthscale;
+  double area_eps = 1e-20 * mesh_priv(m)->lengthscale * mesh_priv(m)->lengthscale;
   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 v0 = m->vertices.items[mesh_priv(m)->face_indices[3 * f]];
+    vector3 v1 = m->vertices.items[mesh_priv(m)->face_indices[3 * f + 1]];
+    vector3 v2 = m->vertices.items[mesh_priv(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 > area_eps) ? vector3_scale(1.0 / len, n) : n;
+    mesh_priv(m)->face_areas[f] = 0.5 * len;
+    mesh_priv(m)->face_normals[f] = (len > area_eps) ? vector3_scale(1.0 / len, n) : n;
   }
 
   /* Compute centroid. */
-  m->centroid.x = m->centroid.y = m->centroid.z = 0;
+  mesh_priv(m)->centroid.x = mesh_priv(m)->centroid.y = mesh_priv(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);
+    mesh_priv(m)->centroid = vector3_plus(mesh_priv(m)->centroid, m->vertices.items[i]);
+  mesh_priv(m)->centroid = vector3_scale(1.0 / nv, mesh_priv(m)->centroid);
 
   /* 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
@@ -2733,11 +2779,11 @@ static void init_mesh(geometric_object *o) {
     /* 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++) {
+    mesh_priv(m)->is_closed = 1;
+    for (int f = 0; f < nf && mesh_priv(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 va = mesh_priv(m)->face_indices[3 * f + e];
+        int vb = mesh_priv(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;
@@ -2752,7 +2798,7 @@ static void init_mesh(geometric_object *o) {
           } 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 */
+            if (ht_cnt[h] > 2) mesh_priv(m)->is_closed = 0; /* non-manifold edge */
             break;
           }
           h = (h + 1) & htmask;
@@ -2760,10 +2806,10 @@ static void init_mesh(geometric_object *o) {
       }
     }
     /* Check that all edges have exactly 2 faces. */
-    if (m->is_closed) {
+    if (mesh_priv(m)->is_closed) {
       for (int i = 0; i < htsize; i++) {
         if (ht_vlo[i] != -1 && ht_cnt[i] != 2) {
-          m->is_closed = 0;
+          mesh_priv(m)->is_closed = 0;
           break;
         }
       }
@@ -2772,7 +2818,7 @@ static void init_mesh(geometric_object *o) {
     free(ht_vhi);
     free(ht_cnt);
 
-    if (!m->is_closed)
+    if (!mesh_priv(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");
   }
@@ -2799,9 +2845,9 @@ static void init_mesh(geometric_object *o) {
     } while (0)
 
     for (int f = 0; f < nf; f++) {
-      int a = m->face_indices[3 * f];
-      int b = m->face_indices[3 * f + 1];
-      int c = m->face_indices[3 * f + 2];
+      int a = mesh_priv(m)->face_indices[3 * f];
+      int b = mesh_priv(m)->face_indices[3 * f + 1];
+      int c = mesh_priv(m)->face_indices[3 * f + 2];
       /* Union a-b. */
       UF_FIND(a); UF_FIND(b);
       if (a != b) {
@@ -2810,9 +2856,9 @@ static void init_mesh(geometric_object *o) {
         if (uf_rank[a] == uf_rank[b]) uf_rank[a]++;
       }
       /* Union a-c (a is already a root after the above). */
-      a = m->face_indices[3 * f];
+      a = mesh_priv(m)->face_indices[3 * f];
       UF_FIND(a);
-      c = m->face_indices[3 * f + 2];
+      c = mesh_priv(m)->face_indices[3 * f + 2];
       UF_FIND(c);
       if (a != c) {
         if (uf_rank[a] < uf_rank[c]) { int t = a; a = c; c = t; }
@@ -2838,21 +2884,21 @@ static void init_mesh(geometric_object *o) {
     double *comp_vol = (double *)calloc(num_components, sizeof(double));
     CHECK(comp_vol, "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]];
-      int ci = comp_id[m->face_indices[3 * f]];
+      vector3 v0 = m->vertices.items[mesh_priv(m)->face_indices[3 * f]];
+      vector3 v1 = m->vertices.items[mesh_priv(m)->face_indices[3 * f + 1]];
+      vector3 v2 = m->vertices.items[mesh_priv(m)->face_indices[3 * f + 2]];
+      int ci = comp_id[mesh_priv(m)->face_indices[3 * f]];
       comp_vol[ci] += vector3_dot(v0, vector3_cross(v1, v2));
     }
 
     /* Flip faces in components with negative signed volume. */
     for (int f = 0; f < nf; f++) {
-      int ci = comp_id[m->face_indices[3 * f]];
+      int ci = comp_id[mesh_priv(m)->face_indices[3 * f]];
       if (comp_vol[ci] < 0) {
-        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;
-        m->face_normals[f] = vector3_scale(-1.0, m->face_normals[f]);
+        int tmp = mesh_priv(m)->face_indices[3 * f + 1];
+        mesh_priv(m)->face_indices[3 * f + 1] = mesh_priv(m)->face_indices[3 * f + 2];
+        mesh_priv(m)->face_indices[3 * f + 2] = tmp;
+        mesh_priv(m)->face_normals[f] = vector3_scale(-1.0, mesh_priv(m)->face_normals[f]);
       }
     }
 
@@ -2864,15 +2910,15 @@ static void init_mesh(geometric_object *o) {
 
   /* 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");
+  mesh_priv(m)->bvh = (mesh_bvh_node *)malloc(max_nodes * sizeof(mesh_bvh_node));
+  CHECK(mesh_priv(m)->bvh, "out of memory");
+  mesh_priv(m)->bvh_face_ids = (int *)malloc(nf * sizeof(int));
+  CHECK(mesh_priv(m)->bvh_face_ids, "out of memory");
   for (int i = 0; i < nf; i++)
-    m->bvh_face_ids[i] = i;
+    mesh_priv(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);
+  mesh_priv(m)->num_bvh_nodes = 0;
+  mesh_bvh_build(m, mesh_priv(m)->bvh_face_ids, 0, nf, mesh_priv(m)->bvh, &mesh_priv(m)->num_bvh_nodes);
 }
 
 /***************************************************************/
@@ -2881,29 +2927,18 @@ static void init_mesh(geometric_object *o) {
 
 static void reinit_mesh(geometric_object *o) {
   mesh *m = o->subclass.mesh_data;
-  /* Skip rebuild if BVH is already built. Mesh vertices are not part of the
-     documented post-construction API, so a non-NULL bvh implies the cached
-     state is still valid. This preserves the fast path for geom_fix_object_ptr
-     on copied meshes (called hundreds of times during meep's init_sim).
-     NOT THREAD-SAFE for first-time init: the m->bvh != NULL check below is
-     unsynchronized, and init_mesh allocates m->bvh BEFORE populating its
-     contents. Concurrent first-time callers would race both on the check
-     (double-init, leak) and on the post-allocation window (second thread
-     reads bvh != NULL and traverses an empty BVH). Callers must ensure the
-     first reinit_mesh / init_mesh on each mesh runs single-threaded — in
-     practice this is done by calling geom_fix_objects() once at geometry
-     setup, before any parallel queries. After init, all subsequent
-     reinit_mesh calls are pure no-op fast-path returns and safe under
-     concurrency. */
-  if (m->bvh != NULL) return;
-  free(m->face_normals);
-  free(m->face_areas);
-  free(m->bvh);
-  free(m->bvh_face_ids);
-  m->face_normals = NULL;
-  m->face_areas = NULL;
-  m->bvh = NULL;
-  m->bvh_face_ids = NULL;
+  /* Fast path: if the internal cache exists, it is fully built and valid
+     (init_mesh only stores m->internal at the end, atomically from a
+     reader's perspective). mesh_copy clears m->internal to NULL on the
+     copy, which causes the first reinit_mesh on that copy to lazily
+     rebuild here. NOT THREAD-SAFE for first-time init on a given mesh:
+     the m->internal == NULL check below is unsynchronized. Callers must
+     ensure the first reinit_mesh / init_mesh on each mesh runs single-
+     threaded — in practice this is done by calling geom_fix_objects()
+     once at geometry setup, before any parallel queries. After init, all
+     subsequent reinit_mesh calls are no-op fast-path returns and safe
+     under concurrency. */
+  if (m->internal != NULL) return;
   init_mesh(o);
 }
 
@@ -2932,15 +2967,22 @@ geometric_object make_mesh_with_center(material_type material, vector3 center,
   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));
-
-  /* Shift vertices before init so BVH is only built once. */
+  /* Pack triangle indices into the public vector3_list. Each triangle's
+     three ints go into x, y, z of a vector3 (representable as doubles
+     up to 2^53). init_mesh will unpack into a flat int[] in the internal
+     cache for fast indexed access. */
+  m->face_indices.num_items = num_triangles;
+  m->face_indices.items = (vector3 *)malloc(num_triangles * sizeof(vector3));
+  CHECK(m->face_indices.items, "out of memory");
+  for (int t = 0; t < num_triangles; t++) {
+    m->face_indices.items[t].x = (double)triangles[3 * t];
+    m->face_indices.items[t].y = (double)triangles[3 * t + 1];
+    m->face_indices.items[t].z = (double)triangles[3 * t + 2];
+  }
+
+  /* Shift vertices before init so the BVH is built around the final
+     positions. */
   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]);
@@ -2950,11 +2992,12 @@ geometric_object make_mesh_with_center(material_type material, vector3 center,
       m->vertices.items[i] = vector3_plus(m->vertices.items[i], shift);
   }
 
-  /* Initialize derived data (normals, closure check, BVH). */
+  /* Initialize derived data (allocates m->internal, computes normals,
+     closure check, BVH). */
   init_mesh(&o);
 
   /* Set center from the (possibly shifted) centroid. */
-  o.center = m->centroid;
+  o.center = mesh_priv(m)->centroid;
 
   return o;
 }
diff --git a/utils/geom.scm b/utils/geom.scm
index cfc17a9..d952371 100644
--- a/utils/geom.scm
+++ b/utils/geom.scm
@@ -155,7 +155,11 @@
 (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-property face_indices '() (make-list-type 'vector3))
+; opaque pointer to mesh_internal (face normals, BVH, etc.) — allocated by
+; init_mesh in C and never touched from Scheme. 'SCM becomes void* in
+; ctlgeom-types.h via the existing sed rule in utils/Makefile.am.
+  (define-property internal '() 'SCM))
 
 (define-class ellipsoid block
   (define-derived-property inverse-semi-axes 'vector3
diff --git a/utils/mesh-internal.h b/utils/mesh-internal.h
new file mode 100644
index 0000000..da40349
--- /dev/null
+++ b/utils/mesh-internal.h
@@ -0,0 +1,45 @@
+/* mesh-internal.h: private C-side cache for the mesh geometric object.
+ *
+ * The public mesh struct (in ctlgeom-types.h, generated from geom.scm) holds
+ * only what's representable in Scheme: the user-supplied vertices + triangle
+ * indices and an opaque void *internal pointer. Everything else — face
+ * normals, areas, BVH acceleration structure, closure flag, centroid,
+ * characteristic lengthscale — lives in mesh_internal, allocated by
+ * init_mesh and reached via mesh_priv(m).
+ *
+ * This header is NOT installed and NOT in geom.scm; it is included only by
+ * geom.c and the mesh tests.
+ */
+
+#ifndef MESH_INTERNAL_H
+#define MESH_INTERNAL_H
+
+#include "ctlgeom-types.h"
+
+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_internal {
+  int            num_faces;
+  int           *face_indices;    /* unpacked flat: 3 ints per triangle */
+  vector3       *face_normals;
+  number        *face_areas;
+  int            num_bvh_nodes;
+  mesh_bvh_node *bvh;
+  int           *bvh_face_ids;
+  boolean        is_closed;
+  vector3        centroid;
+  number         lengthscale;
+} mesh_internal;
+
+static inline mesh_internal *mesh_priv(const mesh *m) {
+  return (mesh_internal *)m->internal;
+}
+
+#endif /* MESH_INTERNAL_H */
diff --git a/utils/test-mesh.c b/utils/test-mesh.c
index 0677538..cce3eb1 100644
--- a/utils/test-mesh.c
+++ b/utils/test-mesh.c
@@ -28,6 +28,7 @@
 #include <string.h>
 
 #include "ctlgeom.h"
+#include "mesh-internal.h"
 
 #define K_PI 3.141592653589793238462643383279502884197
 #define TOLERANCE 1e-6
@@ -441,7 +442,7 @@ static void test_open_mesh(void) {
   };
   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);
+  ASSERT_TRUE("open mesh detected as not closed", !mesh_priv(m)->is_closed);
 
   /* point_in_mesh should return false for open meshes. */
   vector3 p = {0.1, 0.1, 0.1};
@@ -458,12 +459,12 @@ 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);
+  ASSERT_TRUE("cube mesh detected as closed", mesh_priv(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);
+  ASSERT_TRUE("tetra mesh detected as closed", mesh_priv(m)->is_closed);
   geometric_object_destroy(tetra);
   printf("done\n");
 }
@@ -488,7 +489,7 @@ static void test_boundary_edge(void) {
   };
   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);
+  ASSERT_TRUE("boundary edge mesh detected as not closed", !mesh_priv(m)->is_closed);
   geometric_object_destroy(obj);
   printf("done\n");
 }
@@ -514,7 +515,7 @@ static void test_nonmanifold_edge(void) {
   };
   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);
+  ASSERT_TRUE("non-manifold edge mesh detected as not closed", !mesh_priv(m)->is_closed);
   geometric_object_destroy(obj);
   printf("done\n");
 }
@@ -539,7 +540,7 @@ static void test_isolated_vertex(void) {
   };
   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);
+  ASSERT_TRUE("mesh with isolated vertex still closed", mesh_priv(m)->is_closed);
 
   /* Point inside tetrahedron should still work. */
   vector3 p = {0, 0, 0};

From 3bad3395cb7fa0e8458108ce4e5c66c4b9b70e9e Mon Sep 17 00:00:00 2001
From: Luochenghuang <Luocheng@uw.edu>
Date: Fri, 29 May 2026 10:38:06 -0700
Subject: [PATCH 2/5] Eliminate slist_unique VLA in intersect_line_with_prism
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The duplicate-removal step previously wrote deduped values into a stack VLA slist_unique[num_vertices+2] and then did 'slist = slist_unique;' before returning the new count. That assignment is a dead store: slist is a pass-by-value pointer parameter, so reassigning it only rebinds the local variable. The caller's buffer is never touched and slist_unique is destroyed when the function returns, so the caller reads the first num_unique_elements entries of the sorted-but-undeduped buffer instead of the actual deduped values.

In intersect_line_segment_with_prism this corrupts the inside/outside parity walk: a near-duplicate that the dedup loop intended to drop stays in the buffer (producing a phantom thin 'inside' sliver of width ~1e-3*|s|), and a real far-end crossing is silently dropped because the caller stops at the smaller returned count. The bug triggers when two s-values fall within 1e-3*|s| of each other — edge-grazing rays, rays through a shared vertex, near-tangent rays. test-prism's random-ray distribution didn't statistically hit the window often enough at OMP=1 to flag it, but parallel runs (e.g. OMP=8) shuffle the libc rand() sequence and hit the window often enough that the strict 1e-6 segment-length check fails.

The fix eliminates the VLA entirely. The dedup is done in place with a write head 'iv' and read head 'nv': because writes go to slist[iv++] with iv <= nv, position nv-1 is either untouched or was previously written with its own value as a no-op, so reading slist[nv-1] still yields the original sorted value and the comparison semantics are bit-identical to the original intent. Net effect: removes the dead store, removes a potentially huge stack allocation, and the dedup now actually applies to the caller's buffer.

Verified with --enable-openmp: test-prism passes at both OMP_NUM_THREADS=1 (0/10000 points failed, 0/1000 segments failed, 0/65 prism helper failures) and OMP_NUM_THREADS=8 across multiple trials. test-mesh remains clean at OMP=8 as well.
---
 utils/geom.c | 32 ++++++++++++++------------------
 1 file changed, 14 insertions(+), 18 deletions(-)

diff --git a/utils/geom.c b/utils/geom.c
index 0f08949..5379cb6 100644
--- a/utils/geom.c
+++ b/utils/geom.c
@@ -3309,24 +3309,20 @@ int intersect_line_with_prism(prism *prsm, vector3 pc, vector3 dc, double *slist
   if (num_intersections > slist_len) return num_intersections;
 
   qsort((void *)slist, num_intersections, sizeof(double), dcmp);
-  // if num_intersections is zero then just return that
-  if (num_intersections == 0) return num_intersections;
-  else {
-    // remove duplicates from slist
-    double duplicate_tolerance = 1e-3;
-    int num_unique_elements = 1;
-    double slist_unique[num_vertices+2];
-    slist_unique[0] = slist[0];
-    for (nv = 1; nv < num_intersections; nv++) {
-      if (fabs(slist[nv] - slist[nv-1]) > duplicate_tolerance*fabs(slist[nv])) {
-        slist_unique[num_unique_elements] = slist[nv];
-        num_unique_elements++;
-      }
-    }
-    slist = slist_unique;
-    num_intersections = num_unique_elements;
-    return num_intersections;
-  }
+  // Remove near-duplicates in place: walk the sorted slist with a write
+  // head `iv` and read head `nv`; the write at slist[iv++] is always to a
+  // position <= nv, so the comparison slist[nv-1] reads the original
+  // sorted value (positions [iv..nv-1] are either untouched or written
+  // with their own value as a no-op).
+  double duplicate_tolerance = 1e-3;
+  int iv = num_intersections < 1 ? 0 : 1; /* min(1, num_intersections) */
+  for (nv = 1; nv < num_intersections; nv++) {
+    if (fabs(slist[nv] - slist[nv-1]) > duplicate_tolerance*fabs(slist[nv])) {
+      slist[iv++] = slist[nv];
+    }
+  }
+  num_intersections = iv;
+  return num_intersections;
 }
 
 /***************************************************************/

From 3499e017c0f84a6be512882bd23bedccd9dc6316 Mon Sep 17 00:00:00 2001
From: Luochenghuang <Luocheng@uw.edu>
Date: Mon, 1 Jun 2026 11:28:00 -0700
Subject: [PATCH 3/5] Inline mesh-internal.h into geom.c so MPB-style cp builds
 still work

The previous PIMPL commit placed the private mesh_bvh_node / mesh_internal definitions and the mesh_priv() accessor in a separate header utils/mesh-internal.h, included by geom.c. That breaks MPB, meep, and the libctl examples/ directory itself: all of those build by 'cp -f utils/geom.c their_tree/geom.c' and compile the copy in their own build context, without the companion header.

Until those downstreams are taught to link against libctlgeom directly (instead of copying geom.c), geom.c must remain self-contained. This commit inlines the private definitions directly at the top of geom.c (right after the existing includes) and removes the #include 'mesh-internal.h' line. The mesh-internal.h file is deleted.

To preserve test-mesh's ability to verify the closure check without exposing mesh_priv outside geom.c, the is_closed field is promoted to the public mesh struct via geom.scm. It is conceptually a derived public property anyway (users may want to query 'did init_mesh detect this mesh as watertight?'), and the auto-generated mesh_copy / mesh_equal handle it as a plain boolean value with no special code. ctlgeom-types.h and geom-ctl-io.c are updated to match. test-mesh.c drops the #include 'mesh-internal.h' line and reverts mesh_priv(m)->is_closed back to m->is_closed.

Build clean. test-prism and test-mesh both pass at OMP_NUM_THREADS=1 and 8: test-mesh 0 failures (including 0 mismatches in test_cube_vs_block and test_cube_segments_vs_block); test-prism 0/10000 point failures, 0/1000 segment failures, 0/65 prism helper failures.
---
 utils/ctlgeom-types.h |  1 +
 utils/geom-ctl-io.c   |  3 +++
 utils/geom.c          | 53 +++++++++++++++++++++++++++++++++++--------
 utils/geom.scm        |  5 +++-
 utils/mesh-internal.h | 45 ------------------------------------
 utils/test-mesh.c     | 13 +++++------
 6 files changed, 58 insertions(+), 62 deletions(-)
 delete mode 100644 utils/mesh-internal.h

diff --git a/utils/ctlgeom-types.h b/utils/ctlgeom-types.h
index f92be6d..9a2693f 100644
--- a/utils/ctlgeom-types.h
+++ b/utils/ctlgeom-types.h
@@ -108,6 +108,7 @@ extern		"C" {
 	typedef struct mesh_struct {
 		vector3_list	vertices;
 		vector3_list	face_indices;
+		boolean		is_closed;
 		void*		internal;
 	}		mesh;
 
diff --git a/utils/geom-ctl-io.c b/utils/geom-ctl-io.c
index 9e61e6d..da27d8e 100644
--- a/utils/geom-ctl-io.c
+++ b/utils/geom-ctl-io.c
@@ -69,6 +69,7 @@ mesh_copy(const mesh * o0, mesh * o)
 			o->face_indices.items[i_t] = o0->face_indices.items[i_t];
 		}
 	}
+	o->is_closed = o0->is_closed;
 	o->internal = NULL;
 	mesh_init_internal(o);  /* rebuild BVH eagerly; safe under later concurrent queries */
 }
@@ -300,6 +301,8 @@ mesh_equal(const mesh * o0, const mesh * o)
 				return 0;
 		}
 	}
+	if (o->is_closed != o0->is_closed)
+		return 0;
 	;
 	return 1;
 }
diff --git a/utils/geom.c b/utils/geom.c
index 5379cb6..fcfe8a6 100644
--- a/utils/geom.c
+++ b/utils/geom.c
@@ -33,7 +33,42 @@
 static void material_type_copy(void **src, void **dest) { *dest = *src; }
 #endif
 #include "ctlgeom.h"
-#include "mesh-internal.h"
+
+/***************************************************************/
+/* Private mesh internals.                                      */
+/*                                                              */
+/* Inlined here rather than placed in a separate header because */
+/* MPB/meep and the libctl examples/ directory build by copying */
+/* this geom.c file out of the libctl tree, so any companion    */
+/* header would not travel with it. Keep this section in sync   */
+/* with any future split if/when those downstreams are taught   */
+/* to link against libctlgeom directly.                         */
+/***************************************************************/
+
+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_internal {
+  int            num_faces;
+  int           *face_indices;    /* unpacked flat: 3 ints per triangle */
+  vector3       *face_normals;
+  number        *face_areas;
+  int            num_bvh_nodes;
+  mesh_bvh_node *bvh;
+  int           *bvh_face_ids;
+  vector3        centroid;
+  number         lengthscale;
+} mesh_internal;
+
+static inline mesh_internal *mesh_priv(const mesh *m) {
+  return (mesh_internal *)m->internal;
+}
 
 #ifdef CXX_CTL_IO
 using namespace ctlio;
@@ -2559,7 +2594,7 @@ static const vector3 mesh_ray_dirs[3] = {
 };
 
 static boolean point_in_mesh(const mesh *m, vector3 p) {
-  if (!mesh_priv(m)->is_closed) return 0;
+  if (!m->is_closed) return 0;
 
   /* Fast path: cast one ray. If no degenerate edge/vertex hits were
      detected during deduplication, trust the result immediately.
@@ -2616,7 +2651,7 @@ 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, mesh_priv(m)->num_faces,
-             mesh_priv(m)->is_closed ? "closed" : "OPEN (WARNING)");
+             m->is_closed ? "closed" : "OPEN (WARNING)");
 }
 
 static double intersect_line_segment_with_mesh(const mesh *m, vector3 p, vector3 d,
@@ -2779,8 +2814,8 @@ static void init_mesh(geometric_object *o) {
     /* Mark empty slots with -1. */
     for (int i = 0; i < htsize; i++) ht_vlo[i] = -1;
 
-    mesh_priv(m)->is_closed = 1;
-    for (int f = 0; f < nf && mesh_priv(m)->is_closed; f++) {
+    m->is_closed = 1;
+    for (int f = 0; f < nf && m->is_closed; f++) {
       for (int e = 0; e < 3; e++) {
         int va = mesh_priv(m)->face_indices[3 * f + e];
         int vb = mesh_priv(m)->face_indices[3 * f + (e + 1) % 3];
@@ -2798,7 +2833,7 @@ static void init_mesh(geometric_object *o) {
           } else if (ht_vlo[h] == vlo && ht_vhi[h] == vhi) {
             /* Found existing edge. */
             ht_cnt[h]++;
-            if (ht_cnt[h] > 2) mesh_priv(m)->is_closed = 0; /* non-manifold edge */
+            if (ht_cnt[h] > 2) m->is_closed = 0; /* non-manifold edge */
             break;
           }
           h = (h + 1) & htmask;
@@ -2806,10 +2841,10 @@ static void init_mesh(geometric_object *o) {
       }
     }
     /* Check that all edges have exactly 2 faces. */
-    if (mesh_priv(m)->is_closed) {
+    if (m->is_closed) {
       for (int i = 0; i < htsize; i++) {
         if (ht_vlo[i] != -1 && ht_cnt[i] != 2) {
-          mesh_priv(m)->is_closed = 0;
+          m->is_closed = 0;
           break;
         }
       }
@@ -2818,7 +2853,7 @@ static void init_mesh(geometric_object *o) {
     free(ht_vhi);
     free(ht_cnt);
 
-    if (!mesh_priv(m)->is_closed)
+    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");
   }
diff --git a/utils/geom.scm b/utils/geom.scm
index d952371..9c71662 100644
--- a/utils/geom.scm
+++ b/utils/geom.scm
@@ -156,7 +156,10 @@
 ; fields to be filled in by users
   (define-property vertices '() (make-list-type 'vector3))
   (define-property face_indices '() (make-list-type 'vector3))
-; opaque pointer to mesh_internal (face normals, BVH, etc.) — allocated by
+; computed by init_mesh: true if the mesh is watertight (every edge shared
+; by exactly two triangles). Read-only from the user's perspective.
+  (define-property is_closed 0 'boolean)
+; opaque pointer to mesh_internal (face normals, BVH, etc.) allocated by
 ; init_mesh in C and never touched from Scheme. 'SCM becomes void* in
 ; ctlgeom-types.h via the existing sed rule in utils/Makefile.am.
   (define-property internal '() 'SCM))
diff --git a/utils/mesh-internal.h b/utils/mesh-internal.h
deleted file mode 100644
index da40349..0000000
--- a/utils/mesh-internal.h
+++ /dev/null
@@ -1,45 +0,0 @@
-/* mesh-internal.h: private C-side cache for the mesh geometric object.
- *
- * The public mesh struct (in ctlgeom-types.h, generated from geom.scm) holds
- * only what's representable in Scheme: the user-supplied vertices + triangle
- * indices and an opaque void *internal pointer. Everything else — face
- * normals, areas, BVH acceleration structure, closure flag, centroid,
- * characteristic lengthscale — lives in mesh_internal, allocated by
- * init_mesh and reached via mesh_priv(m).
- *
- * This header is NOT installed and NOT in geom.scm; it is included only by
- * geom.c and the mesh tests.
- */
-
-#ifndef MESH_INTERNAL_H
-#define MESH_INTERNAL_H
-
-#include "ctlgeom-types.h"
-
-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_internal {
-  int            num_faces;
-  int           *face_indices;    /* unpacked flat: 3 ints per triangle */
-  vector3       *face_normals;
-  number        *face_areas;
-  int            num_bvh_nodes;
-  mesh_bvh_node *bvh;
-  int           *bvh_face_ids;
-  boolean        is_closed;
-  vector3        centroid;
-  number         lengthscale;
-} mesh_internal;
-
-static inline mesh_internal *mesh_priv(const mesh *m) {
-  return (mesh_internal *)m->internal;
-}
-
-#endif /* MESH_INTERNAL_H */
diff --git a/utils/test-mesh.c b/utils/test-mesh.c
index cce3eb1..0677538 100644
--- a/utils/test-mesh.c
+++ b/utils/test-mesh.c
@@ -28,7 +28,6 @@
 #include <string.h>
 
 #include "ctlgeom.h"
-#include "mesh-internal.h"
 
 #define K_PI 3.141592653589793238462643383279502884197
 #define TOLERANCE 1e-6
@@ -442,7 +441,7 @@ static void test_open_mesh(void) {
   };
   geometric_object open = make_mesh(NULL, verts, 4, tris, 4);
   mesh *m = open.subclass.mesh_data;
-  ASSERT_TRUE("open mesh detected as not closed", !mesh_priv(m)->is_closed);
+  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};
@@ -459,12 +458,12 @@ 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", mesh_priv(m)->is_closed);
+  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", mesh_priv(m)->is_closed);
+  ASSERT_TRUE("tetra mesh detected as closed", m->is_closed);
   geometric_object_destroy(tetra);
   printf("done\n");
 }
@@ -489,7 +488,7 @@ static void test_boundary_edge(void) {
   };
   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", !mesh_priv(m)->is_closed);
+  ASSERT_TRUE("boundary edge mesh detected as not closed", !m->is_closed);
   geometric_object_destroy(obj);
   printf("done\n");
 }
@@ -515,7 +514,7 @@ static void test_nonmanifold_edge(void) {
   };
   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", !mesh_priv(m)->is_closed);
+  ASSERT_TRUE("non-manifold edge mesh detected as not closed", !m->is_closed);
   geometric_object_destroy(obj);
   printf("done\n");
 }
@@ -540,7 +539,7 @@ static void test_isolated_vertex(void) {
   };
   geometric_object obj = make_mesh(NULL, verts, 5, tris, 4);
   mesh *m = obj.subclass.mesh_data;
-  ASSERT_TRUE("mesh with isolated vertex still closed", mesh_priv(m)->is_closed);
+  ASSERT_TRUE("mesh with isolated vertex still closed", m->is_closed);
 
   /* Point inside tetrahedron should still work. */
   vector3 p = {0, 0, 0};

From 62857bd73941df03fa37ba624103c60f39698ed4 Mon Sep 17 00:00:00 2001
From: Luochenghuang <Luocheng@uw.edu>
Date: Mon, 1 Jun 2026 13:47:02 -0700
Subject: [PATCH 4/5] Treat geom-ctl-io.c as checked-in, not BUILT_SOURCES

The PIMPL refactor added two hand-tweaks to utils/geom-ctl-io.c that gen-ctl-io would not emit: extern declarations for mesh_init_internal and mesh_internal_free, plus calls to them inside the auto-generated mesh_copy and mesh_destroy bodies. These hooks complete the opaque mesh_internal lifecycle (eagerly rebuild the cache on mesh_copy so the copy has its own state; free the cache on mesh_destroy).

But geom-ctl-io.c was still classified as a BUILT_SOURCE with an active regen rule. On any build that detects Guile (WITH_GUILE=true, true on the new GitHub Actions CI), the regen rule sed-rewrites geom-ctl-io.c from ctl-io.c and silently overwrites the hand-tweaks. Compilation still succeeds, but the runtime contract breaks: mesh_destroy stops freeing the internal cache (silent leak per mesh) and mesh_copy reverts to a shallow pointer copy (two meshes sharing one mesh_internal allocation, with a dangling-pointer footgun on first destroy). PR #74 didn't expose this because legacy Guile detection in configure.ac fails on modern distros so WITH_GUILE was effectively false; the new CI is the first build where Guile is actually detected.

This commit removes geom-ctl-io.c from BUILT_SOURCES, deletes its regen rule, adds it to EXTRA_DIST so it ships in tarballs, and drops it from clean-local so 'make clean' doesn't wipe it. ctlgeom-types.h stays as-is because the PIMPL refactor doesn't hand-edit it (the regenerated content matches the committed content byte-for-byte).

Verified: 'make -n geom-ctl-io.c' reports 'Nothing to be done' (regen rule successfully gone). test-mesh and test-prism both pass cleanly at OMP_NUM_THREADS=1 and 8.
---
 utils/Makefile.am | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/utils/Makefile.am b/utils/Makefile.am
index 5f523bc..e5e7fe0 100644
--- a/utils/Makefile.am
+++ b/utils/Makefile.am
@@ -10,7 +10,12 @@ nodist_include_HEADERS = ctlgeom-types.h
 lib_LTLIBRARIES = libctlgeom.la
 noinst_PROGRAMS = geomtst
 
-EXTRA_DIST = gen-ctl-io.in README geom.scm geom-ctl-io-defaults.c nlopt.c
+# geom-ctl-io.c is checked in and hand-maintained — it carries PIMPL
+# lifecycle hooks (extern declarations + calls to mesh_init_internal and
+# mesh_internal_free in mesh_copy / mesh_destroy) that gen-ctl-io does
+# not emit. Treating it as a regular distributed source (not BUILT) keeps
+# Guile-detected builds from regenerating and clobbering those hooks.
+EXTRA_DIST = gen-ctl-io.in README geom.scm geom-ctl-io-defaults.c nlopt.c geom-ctl-io.c
 
 libctlgeom_la_SOURCES = geom.c $(top_srcdir)/src/ctl-math.c $(top_srcdir)/src/integrator.c geom-ctl-io.c ctlgeom-types.h
 libctlgeom_la_LDFLAGS = -version-info @SHARED_VERSION_INFO@
@@ -34,7 +39,7 @@ TESTS = test-prism test-mesh
 
 dist_man_MANS = gen-ctl-io.1
 
-BUILT_SOURCES = gen-ctl-io geom-ctl-io.c ctlgeom-types.h nlopt-constants.scm
+BUILT_SOURCES = gen-ctl-io ctlgeom-types.h nlopt-constants.scm
 
 nlopt-constants.scm:
 	echo "#include <nlopt.h>" > nlopt-constants.h
@@ -58,10 +63,7 @@ ctlgeom-types.h: ctl-io.h
 	echo '#  define LIBCTL_BUGFIX_VERSION '$(LIBCTL_BUGFIX_VERSION) >> $@
 	echo '#endif' >> $@
 
-geom-ctl-io.c: ctl-io.c
-	sed 's,ctl-io\.h,ctlgeom-types.h,;s,/.* Input variables .*/,@#include "geom-ctl-io-defaults.c"@#if 0@,;s,/.* Output variables .*/,#endif@,' ctl-io.c | tr '@' '\n' > $@
-
 endif
 
 clean-local:
-	rm -f ctl-io.[ch] nlopt-constants.scm ctlgeom-types.h geom-ctl-io.c
+	rm -f ctl-io.[ch] nlopt-constants.scm ctlgeom-types.h

From 50e549803d4f926f1d09e946e383f34f3ce551eb Mon Sep 17 00:00:00 2001
From: Luochenghuang <Luocheng@uw.edu>
Date: Thu, 4 Jun 2026 10:57:34 -0700
Subject: [PATCH 5/5] Refresh ctlgeom-types.h and geom-ctl-io.c against current
 geom.scm

PR #74 committed both ctlgeom-types.h and geom-ctl-io.c, but the versions it shipped were generated from a pre-PR-#75 geom.scm that still had the prism workspace field. After #75 removed workspace from geom.scm and from the C code in geom.c, the committed copies became stale. Without-Guile builds silently kept compiling because the two stale files were consistent with each other (prism_copy/destroy/equal operated on a workspace field that nothing else touched). With-Guile builds papered over it by regenerating ctlgeom-types.h to match the post-#75 geom.scm, which silently differed from the committed copy.

The previous commit (treat geom-ctl-io.c as checked-in) exposed the inconsistency by removing the regeneration safety net: with WITH_GUILE=true, ctlgeom-types.h gets regenerated fresh from geom.scm (no workspace, no number_list typedef), while the committed geom-ctl-io.c still tried to access prism->workspace, producing 'has no member named workspace' compile errors in prism_copy / prism_equal / prism_destroy.

This commit aligns both committed files with current geom.scm:
  - ctlgeom-types.h: replaced with the gen-ctl-io output. Drops the unused number_list typedef and the workspace field from prism struct. Same mesh struct as before (vertices, face_indices, is_closed, internal).
  - geom-ctl-io.c: removed the stale workspace handling from prism_copy / prism_equal / prism_destroy. Kept the PIMPL hand-tweaks (extern mesh_init_internal / mesh_internal_free + their calls in mesh_copy / mesh_destroy) since the file is no longer auto-regenerated.

Verified locally in both build modes:
  - WITHOUT Guile (WITH_GUILE=false, no regen): make + make check + OMP=1/8 clean.
  - WITH Guile (WITH_GUILE=true, ctlgeom-types.h regenerates and matches committed): make + make install + examples/ build (which exercises the cp utils/geom.c -> examples/geom.c MPB-style copy path) + make check + OMP=1/8 clean.
---
 utils/ctlgeom-types.h | 345 ++++++++++++++++++++----------------------
 utils/geom-ctl-io.c   |  23 ---
 2 files changed, 167 insertions(+), 201 deletions(-)

diff --git a/utils/ctlgeom-types.h b/utils/ctlgeom-types.h
index 9a2693f..6322150 100644
--- a/utils/ctlgeom-types.h
+++ b/utils/ctlgeom-types.h
@@ -7,188 +7,177 @@
 #include <ctl-math.h>
 
 #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_struct {
-		vector3_list	vertices;
-		vector3_list	face_indices;
-		boolean		is_closed;
-		void*		internal;
-	}		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);
+extern "C" {
+#endif                          /* __cplusplus */
+
+/******* Type declarations *******/
+
+typedef struct geometric_object_struct {
+void* material;
+vector3 center;
+enum { GEOMETRIC_OBJECT_SELF, MESH, PRISM, BLOCK, SPHERE, CYLINDER, COMPOUND_GEOMETRIC_OBJECT } which_subclass;
+union {
+struct mesh_struct *mesh_data;
+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;
+} 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 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;
+matrix3x3 m_c2p;
+matrix3x3 m_p2c;
+} prism;
+
+typedef struct mesh_struct {
+vector3_list vertices;
+vector3_list face_indices;
+boolean is_closed;
+void* internal;
+} 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 */
+}                               /* extern "C" */
+#endif                          /* __cplusplus */
+
+#endif                          /* CTL_IO_H */
 
-#endif				/* CTL_IO_H */
 #ifndef LIBCTL_MAJOR_VERSION
 #  define LIBCTL_MAJOR_VERSION 4
-#  define LIBCTL_MINOR_VERSION 5
-#  define LIBCTL_BUGFIX_VERSION 1
+#  define LIBCTL_MINOR_VERSION 6
+#  define LIBCTL_BUGFIX_VERSION 0
 #endif
diff --git a/utils/geom-ctl-io.c b/utils/geom-ctl-io.c
index da27d8e..8ac80e9 100644
--- a/utils/geom-ctl-io.c
+++ b/utils/geom-ctl-io.c
@@ -129,14 +129,6 @@ prism_copy(const prism * o0, prism * o)
 		}
 	}
 	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;
 }
@@ -372,15 +364,6 @@ prism_equal(const prism * o0, const prism * o)
 	}
 	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))
@@ -574,12 +557,6 @@ prism_destroy(prism o)
 		}
 	}
 	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