From a7ab1db8960178d1a73bafd434bc39122825050d Mon Sep 17 00:00:00 2001 From: Luochenghuang Date: Tue, 7 Apr 2026 15:44:27 -0700 Subject: [PATCH 1/4] Add Mesh geometry type: Python class, STL/OBJ parsers, C++ binding - Mesh class in geom.py: vertices + triangles, optional center, from_stl() and from_obj() classmethods with built-in parsers (binary/ASCII STL with vertex welding, OBJ with quad triangulation) - pymesh_to_mesh in typemap_utils.cpp: extracts vertices and triangles from Python, calls make_mesh_with_center from libctl - Mesh exported from meep.__init__ (requires manual addition to generated __init__.py until codegen is updated) Requires libctl with MESH support (Luochenghuang/libctl#mesh-geometry-import). --- python/geom.py | 167 +++++++++++++++++++++++++++++++++++++++ python/typemap_utils.cpp | 58 ++++++++++++++ 2 files changed, 225 insertions(+) diff --git a/python/geom.py b/python/geom.py index 26f7b5270..dae2430b9 100755 --- a/python/geom.py +++ b/python/geom.py @@ -1343,6 +1343,173 @@ def __init__( super().__init__(center=center, **kwargs) +class Mesh(GeometricObject): + """ + Triangulated 3D mesh object. + + The mesh must be a closed (watertight) manifold with consistently + oriented triangles. Open meshes will produce a warning and may + give incorrect results for subpixel smoothing. + """ + + def __init__(self, vertices, triangles, center=None, **kwargs): + """ + Construct a `Mesh`. + + + **`vertices` [list of `Vector3`]** — Vertex positions. + + + **`triangles` [list of tuples of 3 ints]** — Triangle vertex indices + (0-based). Each tuple defines a face with outward normal determined + by the right-hand rule. + + + **`center` [`Vector3`, optional]** — If specified, the mesh is + translated so its centroid is at this position. If omitted, the + centroid of the vertices is used. + """ + self.vertices = [Vector3(*v) for v in vertices] + self.triangles = [tuple(t) for t in triangles] + + centroid = sum(self.vertices, Vector3(0)) * (1.0 / len(self.vertices)) + if center is not None: + center = Vector3(*center) + shift = center - centroid + self.vertices = [v + shift for v in self.vertices] + else: + center = centroid + + super().__init__(center=center, **kwargs) + + @classmethod + def from_stl(cls, filename, material=None, center=None, scale=1.0): + """Load mesh from STL file (binary or ASCII).""" + vertices, triangles = _parse_stl(filename, scale) + kwargs = {} + if material is not None: + kwargs["material"] = material + return cls(vertices=vertices, triangles=triangles, + center=center, **kwargs) + + @classmethod + def from_obj(cls, filename, material=None, center=None, scale=1.0): + """Load mesh from Wavefront OBJ file.""" + vertices, triangles = _parse_obj(filename, scale) + kwargs = {} + if material is not None: + kwargs["material"] = material + return cls(vertices=vertices, triangles=triangles, + center=center, **kwargs) + + +def _parse_stl(filename, scale=1.0): + """Parse an STL file (binary or ASCII) and return (vertices, triangles). + + STL stores per-face vertices without sharing, so vertex welding is + performed to build a shared vertex array and triangle index array. + """ + import struct + + with open(filename, "rb") as f: + header = f.read(80) + data = f.read() + + # Try binary format first + is_binary = True + try: + num_tris = struct.unpack(" 0 else 1e8 + vertex_map = {} + vertices = [] + triangles = [] + + def _quantize(v): + return ( + int(round(v[0] * inv_tol)), + int(round(v[1] * inv_tol)), + int(round(v[2] * inv_tol)), + ) + + for i in range(0, len(raw_verts), 3): + tri = [] + for j in range(3): + v = raw_verts[i + j] + key = _quantize(v) + if key not in vertex_map: + vertex_map[key] = len(vertices) + vertices.append(Vector3(v[0], v[1], v[2])) + tri.append(vertex_map[key]) + triangles.append(tuple(tri)) + + return vertices, triangles + + +def _parse_obj(filename, scale=1.0): + """Parse a Wavefront OBJ file and return (vertices, triangles). + + Quad faces are automatically triangulated. OBJ is 1-indexed; + indices are converted to 0-indexed. + """ + vertices = [] + triangles = [] + + with open(filename, "r") as f: + for line in f: + line = line.strip() + if not line or line.startswith("#"): + continue + parts = line.split() + if parts[0] == "v" and len(parts) >= 4: + x, y, z = float(parts[1]), float(parts[2]), float(parts[3]) + vertices.append(Vector3(x * scale, y * scale, z * scale)) + elif parts[0] == "f": + # OBJ face indices may include texture/normal refs: v/vt/vn + face_indices = [] + for p in parts[1:]: + idx = int(p.split("/")[0]) + # OBJ is 1-indexed; handle negative indices + if idx < 0: + idx = len(vertices) + idx + else: + idx -= 1 + face_indices.append(idx) + # Triangulate (fan triangulation for convex polygons) + for i in range(1, len(face_indices) - 1): + triangles.append( + (face_indices[0], face_indices[i], face_indices[i + 1]) + ) + + return vertices, triangles + + class Matrix: """ The `Matrix` class represents a 3x3 matrix with c1, c2, and c3 as its columns. diff --git a/python/typemap_utils.cpp b/python/typemap_utils.cpp index 5352bf1eb..d4b827c51 100644 --- a/python/typemap_utils.cpp +++ b/python/typemap_utils.cpp @@ -1038,6 +1038,63 @@ static int pyprism_to_prism(PyObject *py_prism, geometric_object *p) { ; } +static int pymesh_to_mesh(PyObject *py_mesh, geometric_object *o) { + SWIG_PYTHON_THREAD_SCOPED_BLOCK; + material_type material; + vector3 center; + + if (!get_attr_material(py_mesh, &material) || !get_attr_v3(py_mesh, ¢er, "center")) { + return 0; + } + + PyObject *py_verts = PyObject_GetAttrString(py_mesh, "vertices"); + if (!py_verts || !PyList_Check(py_verts)) { + meep::abort("Expected Mesh.vertices to be a list\n"); + return 0; + } + + int num_vertices = PyList_Size(py_verts); + vector3 *vertices = new vector3[num_vertices]; + for (Py_ssize_t i = 0; i < num_vertices; ++i) { + if (!pyv3_to_v3(PyList_GetItem(py_verts, i), &vertices[i])) { + delete[] vertices; + Py_DECREF(py_verts); + return 0; + } + } + Py_DECREF(py_verts); + + PyObject *py_tris = PyObject_GetAttrString(py_mesh, "triangles"); + if (!py_tris || !PyList_Check(py_tris)) { + delete[] vertices; + meep::abort("Expected Mesh.triangles to be a list\n"); + return 0; + } + + int num_triangles = PyList_Size(py_tris); + int *triangles = new int[3 * num_triangles]; + for (Py_ssize_t i = 0; i < num_triangles; ++i) { + PyObject *tri = PyList_GetItem(py_tris, i); + if (!PyTuple_Check(tri) || PyTuple_Size(tri) != 3) { + delete[] vertices; + delete[] triangles; + Py_DECREF(py_tris); + meep::abort("Each triangle must be a tuple of 3 integers\n"); + return 0; + } + for (int j = 0; j < 3; ++j) + triangles[3 * i + j] = (int)PyLong_AsLong(PyTuple_GetItem(tri, j)); + } + Py_DECREF(py_tris); + + *o = make_mesh_with_center(material, center, vertices, num_vertices, triangles, num_triangles); + + delete[] vertices; + delete[] triangles; + + return 1; +} + static int py_gobj_to_gobj(PyObject *po, geometric_object *o) { SWIG_PYTHON_THREAD_SCOPED_BLOCK; int success = 0; @@ -1050,6 +1107,7 @@ static int py_gobj_to_gobj(PyObject *po, geometric_object *o) { else if (go_type == "Block") { success = pyblock_to_block(po, o); } else if (go_type == "Ellipsoid") { success = pyellipsoid_to_ellipsoid(po, o); } else if (go_type == "Prism") { success = pyprism_to_prism(po, o); } + else if (go_type == "Mesh") { success = pymesh_to_mesh(po, o); } else { PyErr_Format(PyExc_TypeError, "Error: %s is not a valid GeometricObject type", go_type.c_str()); return 0; From fb852889d05578a625835d33619093885b95fd76 Mon Sep 17 00:00:00 2001 From: Luochenghuang Date: Tue, 7 Apr 2026 16:04:46 -0700 Subject: [PATCH 2/4] Add mesh geometry integration test Two tests (~1s total): 1. Geometric convergence: mesh icosphere L2 error vs native Sphere decreases with subdivision (20 faces -> 5120 faces -> exact match) 2. Subpixel smoothing: eps_averaging=True produces intermediate epsilon values at mesh-air interfaces, confirming Kottke averaging works for mesh objects --- python/geom.py | 8 +- python/tests/test_mesh.py | 207 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 210 insertions(+), 5 deletions(-) create mode 100644 python/tests/test_mesh.py diff --git a/python/geom.py b/python/geom.py index dae2430b9..116afc5f1 100755 --- a/python/geom.py +++ b/python/geom.py @@ -1386,8 +1386,7 @@ def from_stl(cls, filename, material=None, center=None, scale=1.0): kwargs = {} if material is not None: kwargs["material"] = material - return cls(vertices=vertices, triangles=triangles, - center=center, **kwargs) + return cls(vertices=vertices, triangles=triangles, center=center, **kwargs) @classmethod def from_obj(cls, filename, material=None, center=None, scale=1.0): @@ -1396,8 +1395,7 @@ def from_obj(cls, filename, material=None, center=None, scale=1.0): kwargs = {} if material is not None: kwargs["material"] = material - return cls(vertices=vertices, triangles=triangles, - center=center, **kwargs) + return cls(vertices=vertices, triangles=triangles, center=center, **kwargs) def _parse_stl(filename, scale=1.0): @@ -1481,7 +1479,7 @@ def _parse_obj(filename, scale=1.0): vertices = [] triangles = [] - with open(filename, "r") as f: + with open(filename) as f: for line in f: line = line.strip() if not line or line.startswith("#"): diff --git a/python/tests/test_mesh.py b/python/tests/test_mesh.py new file mode 100644 index 000000000..c8ebc3fb4 --- /dev/null +++ b/python/tests/test_mesh.py @@ -0,0 +1,207 @@ +"""Integration tests for the Mesh geometric object type.""" + +import meep as mp +import numpy as np +import struct +import tempfile +import os +import unittest + + +def make_cube_mesh(material, center=mp.Vector3(0, 0, 0)): + """Create a unit cube mesh centered at the given point.""" + verts = [ + (-0.5, -0.5, -0.5), + (0.5, -0.5, -0.5), + (0.5, 0.5, -0.5), + (-0.5, 0.5, -0.5), + (-0.5, -0.5, 0.5), + (0.5, -0.5, 0.5), + (0.5, 0.5, 0.5), + (-0.5, 0.5, 0.5), + ] + tris = [ + (0, 2, 1), + (0, 3, 2), # -z + (4, 5, 6), + (4, 6, 7), # +z + (0, 1, 5), + (0, 5, 4), # -y + (2, 3, 7), + (2, 7, 6), # +y + (0, 4, 7), + (0, 7, 3), # -x + (1, 2, 6), + (1, 6, 5), # +x + ] + return mp.Mesh(vertices=verts, triangles=tris, center=center, material=material) + + +def write_cube_stl(path): + """Write a unit cube as a binary STL file.""" + faces = [ + ((-0.5, -0.5, -0.5), (0.5, 0.5, -0.5), (0.5, -0.5, -0.5)), + ((-0.5, -0.5, -0.5), (-0.5, 0.5, -0.5), (0.5, 0.5, -0.5)), + ((-0.5, -0.5, 0.5), (0.5, -0.5, 0.5), (0.5, 0.5, 0.5)), + ((-0.5, -0.5, 0.5), (0.5, 0.5, 0.5), (-0.5, 0.5, 0.5)), + ((-0.5, -0.5, -0.5), (0.5, -0.5, -0.5), (0.5, -0.5, 0.5)), + ((-0.5, -0.5, -0.5), (0.5, -0.5, 0.5), (-0.5, -0.5, 0.5)), + ((0.5, 0.5, -0.5), (-0.5, 0.5, -0.5), (-0.5, 0.5, 0.5)), + ((0.5, 0.5, -0.5), (-0.5, 0.5, 0.5), (0.5, 0.5, 0.5)), + ((-0.5, -0.5, -0.5), (-0.5, -0.5, 0.5), (-0.5, 0.5, 0.5)), + ((-0.5, -0.5, -0.5), (-0.5, 0.5, 0.5), (-0.5, 0.5, -0.5)), + ((0.5, -0.5, -0.5), (0.5, 0.5, -0.5), (0.5, 0.5, 0.5)), + ((0.5, -0.5, -0.5), (0.5, 0.5, 0.5), (0.5, -0.5, 0.5)), + ] + with open(path, "wb") as f: + f.write(b"\x00" * 80) + f.write(struct.pack(" 6) + total = eps.size + fill_fraction = inside / total + expected = 1.0 / 8.0 # unit cube in 2x2x2 cell + self.assertAlmostEqual(fill_fraction, expected, delta=0.02) + + def test_eps_averaging(self): + """Subpixel smoothing produces intermediate epsilon at interfaces.""" + mat = mp.Medium(epsilon=12) + cell = mp.Vector3(2, 2, 2) + cube = make_cube_mesh(mat) + + sim_off = mp.Simulation( + cell_size=cell, + geometry=[cube], + resolution=4, + dimensions=3, + eps_averaging=False, + ) + sim_off.init_sim() + eps_off = np.real( + sim_off.get_array(center=mp.Vector3(), size=cell, component=mp.Dielectric) + ) + + sim_on = mp.Simulation( + cell_size=cell, + geometry=[cube], + resolution=4, + dimensions=3, + eps_averaging=True, + ) + sim_on.init_sim() + eps_on = np.real( + sim_on.get_array(center=mp.Vector3(), size=cell, component=mp.Dielectric) + ) + + unique_off = len(np.unique(np.round(eps_off, 2))) + unique_on = len(np.unique(np.round(eps_on, 2))) + self.assertGreater( + unique_on, + unique_off, + "eps_averaging should produce more unique epsilon values", + ) + + def test_mesh_vs_block(self): + """Cube mesh epsilon grid matches native Block.""" + mat = mp.Medium(epsilon=12) + cell = mp.Vector3(2, 2, 2) + n = 30 + xtics = np.linspace(-1, 1, n) + ytics = np.linspace(-1, 1, n) + ztics = np.array([0.0]) + + cube_mesh = make_cube_mesh(mat) + sim_mesh = mp.Simulation( + cell_size=cell, + geometry=[cube_mesh], + resolution=4, + dimensions=3, + eps_averaging=False, + ) + sim_mesh.init_sim() + eps_mesh = np.real(sim_mesh.get_epsilon_grid(xtics, ytics, ztics)).flatten() + + block = mp.Block(size=mp.Vector3(1, 1, 1), material=mat) + sim_block = mp.Simulation( + cell_size=cell, + geometry=[block], + resolution=4, + dimensions=3, + eps_averaging=False, + ) + sim_block.init_sim() + eps_block = np.real(sim_block.get_epsilon_grid(xtics, ytics, ztics)).flatten() + + mismatches = np.sum(np.abs(eps_mesh - eps_block) > 0.5) + self.assertEqual( + mismatches, + 0, + f"Cube mesh should match Block exactly, got {mismatches} mismatches", + ) + + +if __name__ == "__main__": + unittest.main() From 6306b692f92040aabfd44c8af3005e0d179f0604 Mon Sep 17 00:00:00 2001 From: Luochenghuang Date: Thu, 18 Jun 2026 11:29:20 -0700 Subject: [PATCH 3/4] Export Mesh through SWIG codegen import list Add Mesh to the 'from .geom import (...)' block in meep.i so it is exported into the meep namespace via the normal codegen path, removing the need for the manual __init__.py edit noted in the original commit. Now that the libctl MESH support is upstream (NanoComp/libctl#74), the binding works against the released libctl with no further changes. --- python/meep.i | 1 + 1 file changed, 1 insertion(+) diff --git a/python/meep.i b/python/meep.i index 47fa6df8b..17e06c956 100644 --- a/python/meep.i +++ b/python/meep.i @@ -1631,6 +1631,7 @@ PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where MaterialGrid, Matrix, Medium, + Mesh, MultilevelAtom, NoisyDrudeSusceptibility, NoisyLorentzianSusceptibility, From 7b87bcd12dd0febc69caca423b90b07ff5270882 Mon Sep 17 00:00:00 2001 From: Luochenghuang Date: Fri, 26 Jun 2026 11:47:57 -0700 Subject: [PATCH 4/4] CI: rebuild cached libctl when its revision changes libctl's default branch now carries the mesh geometric object (and the C++ fixes) the Mesh binding needs, so the plain checkout + build is fine. The actual failure was the dependency cache: ~/local was keyed only on the workflow-file hash, so CI kept restoring a pre-mesh libctl and never rebuilt it. Add libctl's checked-out commit to the deps-cache key so a new libctl revision rebuilds the cached deps instead of restoring a stale build. The libctl checkout and autogen build steps are otherwise unchanged. --- .github/workflows/build-ci-adjoint.yml | 10 +++++++++- .github/workflows/build-ci.yml | 10 +++++++++- 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/.github/workflows/build-ci-adjoint.yml b/.github/workflows/build-ci-adjoint.yml index 4189fd3c5..494368d94 100644 --- a/.github/workflows/build-ci-adjoint.yml +++ b/.github/workflows/build-ci-adjoint.yml @@ -35,6 +35,12 @@ jobs: repository: NanoComp/libctl path: libctl-src + - name: Record libctl revision for the cache key + id: libctl + # Track libctl's checked-out commit so a new libctl revision invalidates + # the cached ~/local instead of restoring a stale (pre-mesh) build. + run: echo "sha=$(git -C libctl-src rev-parse HEAD)" >> "$GITHUB_OUTPUT" + - name: Checkout harminv repository uses: actions/checkout@v6 with: @@ -52,7 +58,9 @@ jobs: id: deps-cache with: path: ~/local - key: deps-adjoint-${{ runner.os }}-${{ hashFiles('.github/workflows/build-ci-adjoint.yml') }} + # libctl's commit is part of the key so a new libctl revision rebuilds + # the cached deps instead of restoring a stale (pre-mesh) build. + key: deps-adjoint-${{ runner.os }}-libctl-${{ steps.libctl.outputs.sha }}-${{ hashFiles('.github/workflows/build-ci-adjoint.yml') }} - name: Build and install libctl if: steps.deps-cache.outputs.cache-hit != 'true' diff --git a/.github/workflows/build-ci.yml b/.github/workflows/build-ci.yml index 461024786..028093d2e 100644 --- a/.github/workflows/build-ci.yml +++ b/.github/workflows/build-ci.yml @@ -57,6 +57,12 @@ jobs: repository: NanoComp/libctl path: libctl-src + - name: Record libctl revision for the cache key + id: libctl + # Track libctl's checked-out commit so a new libctl revision invalidates + # the cached ~/local instead of restoring a stale (pre-mesh) build. + run: echo "sha=$(git -C libctl-src rev-parse HEAD)" >> "$GITHUB_OUTPUT" + - name: Checkout harminv repository uses: actions/checkout@v6 with: @@ -80,7 +86,9 @@ jobs: id: deps-cache with: path: ~/local - key: deps-${{ runner.os }}-${{ hashFiles('.github/workflows/build-ci.yml') }} + # libctl's commit is part of the key so a new libctl revision rebuilds + # the cached deps instead of restoring a stale (pre-mesh) build. + key: deps-${{ runner.os }}-libctl-${{ steps.libctl.outputs.sha }}-${{ hashFiles('.github/workflows/build-ci.yml') }} - name: Build and install libctl if: steps.deps-cache.outputs.cache-hit != 'true'