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..5379cb6 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;
 }
@@ -3266,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;
 }
 
 /***************************************************************/
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};