Mesh geometry import

[Luochenghuang]

Summary

Closes #67. Adds a new MESH geometric object type for importing arbitrary triangulated 3D surfaces (e.g., from STL files) into libctl's geometry system.

Design

The mesh follows the same tagged-union dispatch as SPHERE, CYLINDER, BLOCK, and PRISM. A flat-array BVH with SAH binning accelerates all queries to O(log N).

All seven core operations are implemented:

Operation	Implementation
point_in_mesh	Ray casting with lazy 3-ray majority vote (Moller-Trumbore + BVH)
normal_to_mesh	BVH nearest-face search (Eberly closest-point-on-triangle)
intersect_line_segment_with_mesh	Sorted ray-mesh intersections with even-odd parity accumulation
get_mesh_volume	Divergence theorem (signed tetrahedron volumes)
get_mesh_bounding_box	Cached from BVH root node
display_mesh_info	Vertex/face counts and closure status
init_mesh	Validates, computes normals, checks closure, fixes winding, builds BVH

Meep's subpixel smoothing works automatically. No meep C++ changes required.

API

// Closed triangulated 3D mesh. The mesh must be watertight (every edge
// shared by exactly 2 faces); open meshes are detected and produce a warning.
geometric_object make_mesh(material_type material,
                           const vector3 *vertices, int num_vertices,
                           const int *triangles, int num_triangles);

geometric_object make_mesh_with_center(material_type material, vector3 center,
                                       const vector3 *vertices, int num_vertices,
                                       const int *triangles, int num_triangles);

Robustness

• Lazy majority vote: casts 1 ray; only falls back to 3-ray vote when deduplication detects a shared-edge hit. Correct and fast for the common case.
• Mesh closure validation: edge-adjacency hash table detects boundary and non-manifold edges. point_in_mesh returns false for open meshes.
• reinit_mesh guard: skips redundant BVH rebuilds in geom_fix_object_ptr (mesh uses absolute coordinates, not lattice-relative).

Test plan

• 16 unit tests: point-in, volume, bbox, normals, line-segment intersection, copy/destroy, closure detection (open, boundary edge, non-manifold, isolated vertex)
• Cube mesh vs Block primitive: 10K random points, 0 mismatches

[stevengj]

This should be 1e-15 * lengthscale where lengthscale is some characteristic scale, e.g. the diameter of the whole objects.

Similarly, parent_area = 1e-30 * lengthscale * lengthscale.

[stevengj]

These "magic small numbers" should be relative to a lengthscale, e.g. the diameter of the object.

[Luochenghuang]

Hi Prof. @stevengj,

I wanted to get your preference on how to handle STL/mesh file loading on the Python side.

---

Option A: Lightweight custom parser (zero dependencies)

• I've written an 85-line STL parser using only Python's struct module, no pip dependencies at all
• Supports both ASCII and binary STL, deduplicates vertices, returns arrays ready for make_mesh()
• Pros: No new dependencies, simple, easy to audit, ships with libctl
• Cons: Only handles STL

---

Option B: Use an existing library (e.g. trimesh, meshio, or numpy-stl)

• trimesh: full mesh processing (40+ formats, boolean ops, repair, validation). Heavy dependency tree
• Pros: Broad format support, battle-tested, handles edge cases we might miss, trimesh can also validate/repair meshes before passing to libctl
• Cons: Adds a dependency to meep's install, potential version conflicts, overkill if users only need STL

---

Which approach do you prefer? Also happy to hear if there are specific mesh formats beyond STL that meep users commonly work with.

Thanks!

[stevengj]

Let's just have the C api in libctl, and in meep a corresponding MeshObject constructor that takes a list of triangles. Then we can just document an example about how to use an existing library to load the list of triangles from a file and pass it to the MeshObject constructor. This way we don't need to add an explicit dependency to Meep itself — the user can call whatever library they want.

[stevengj]

I think if (m->bvh != NULL) return; is sufficient here since we provide no documented way for the caller to mutate the vertices, and none of our other functions do so. So you can get rid of all the vertex_hash stuff.

(Otherwise, a hash collision could create a bug here.)