Open
Conversation
Implement a new MESH primitive alongside SPHERE, CYLINDER, BLOCK, and PRISM. A mesh is defined by a vertex array and triangle index array, with a BVH (bounding volume hierarchy) for O(log N) queries. Core operations implemented: - point_in_mesh: ray casting with parity counting and deduplication - normal_to_mesh: BVH-accelerated closest-face search - get_mesh_bounding_box: from BVH root node - get_mesh_volume: divergence theorem (signed tetrahedron volumes) - intersect_line_segment_with_mesh: for subpixel smoothing integration - init_mesh: validates geometry, computes normals, fixes winding, builds BVH BVH uses SAH (surface area heuristic) binned construction with flat array layout. Ray-triangle intersection uses Moller-Trumbore algorithm. Closest-point queries use Eberly's algorithm with BVH pruning. Also adds make_mesh/make_mesh_with_center constructors, mesh copy/equal/destroy in geom-ctl-io.c, and 11 unit tests covering point-in, volume, bbox, normals, line-segment intersection, and comparison against the Block primitive.
Two fixes: 1. point_in_mesh: use 3-ray majority vote instead of a single ray. A single ray can pass through a shared mesh edge, causing the deduplication to either miss or double-count a crossing and flip the inside/outside parity. Three rays in irrational directions make it extremely unlikely that more than one hits a degenerate case at the same query point. 2. reinit_mesh: skip BVH rebuild if already initialized. Unlike blocks/cylinders, mesh uses absolute coordinates and doesn't depend on the lattice basis. geom_fix_object_ptr is called hundreds of times during meep's init_sim via geometry copies; without the guard, each call rebuilt the BVH (~150ms for 9K triangles), causing ~40s of wasted time.
- point_in_mesh: cast 1 ray first, only fall back to 3-ray majority vote if deduplication removed near-duplicate t-values (indicating an edge/vertex hit). 3x faster for the common case (no edge hits). - init_mesh: validate mesh closure by building an edge-to-face-count hash table. Every edge must be shared by exactly 2 faces. Sets is_closed=false and prints a warning for open/non-manifold meshes. point_in_mesh returns false for all points on open meshes. - Add 2 new tests: open mesh detection, closed mesh detection (13 tests total, all passing).
Test boundary edges (1 face per edge), non-manifold edges (4 faces sharing an edge), and isolated vertices (vertex unreferenced by any face). 16 tests total, all passing.
1. make_mesh_with_center: shift vertices before calling init_mesh so the BVH is only built once (was building twice when center was specified). 2. intersect_line_segment_with_mesh: determine inside/outside at the segment start from the full intersection list parity (count crossings before parameter 'a') instead of calling point_in_mesh as a fallback. Eliminates millions of extra ray casts during subpixel smoothing. ~20% faster set_epsilon on the Utah teapot (28.6s → 22.9s).
- Remove unused bvh_node_box helper function - find_closest_face: push nearer BVH child last so it's popped first, improving pruning and reducing nodes visited for normal queries - Extract mesh_triangle_vertices helper to avoid redundant vertex fetches. Add mesh_triangle_bbox_centroid to compute both in one pass during BVH binning (was fetching 3 vertices twice per face)
These files are normally auto-generated from geom.scm via gen-ctl-io (requires Guile). Including them allows building libctlgeom without Guile, which is the common case for meep's Python-only builds.
006ae18 to
b14367f
Compare
stevengj
reviewed
Apr 10, 2026
| case 1: lo = node_box.low.y; hi = node_box.high.y; break; | ||
| default: lo = node_box.low.z; hi = node_box.high.z; break; | ||
| } | ||
| if (hi - lo < 1e-15) continue; |
Collaborator
There was a problem hiding this comment.
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
reviewed
Apr 10, 2026
| vector3 n = vector3_cross(e1, e2); | ||
| double len = vector3_norm(n); | ||
| m->face_areas[f] = 0.5 * len; | ||
| m->face_normals[f] = (len > 1e-30) ? vector3_scale(1.0 / len, n) : n; |
Collaborator
There was a problem hiding this comment.
These "magic small numbers" should be relative to a lengthscale, e.g. the diameter of the object.
stevengj
reviewed
Apr 10, 2026
| vector3 e2 = vector3_minus(v2, v0); | ||
| vector3 n = vector3_cross(e1, e2); | ||
| double len = vector3_norm(n); | ||
| m->face_areas[f] = 0.5 * len; |
Collaborator
There was a problem hiding this comment.
Would be good to document that
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Closes #67. Adds a new
MESHgeometric 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, andPRISM. A flat-array BVH with SAH binning accelerates all queries to O(log N).All seven core operations are implemented:
point_in_meshnormal_to_meshintersect_line_segment_with_meshget_mesh_volumeget_mesh_bounding_boxdisplay_mesh_infoinit_meshMeep's subpixel smoothing works automatically. No meep C++ changes required.
API
Robustness
point_in_meshreturns false for open meshes.reinit_meshguard: skips redundant BVH rebuilds ingeom_fix_object_ptr(mesh uses absolute coordinates, not lattice-relative).Test plan
Blockprimitive: 10K random points, 0 mismatches