Replace sphere convergence test with fast test_mesh.py

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 (23 unique values vs 4
   without smoothing), confirming Kottke averaging works for meshes

Author: Luochenghuang

utils/test_mesh.py

@@ -0,0 +1,113 @@
+"""Fast convergence test: mesh icosphere vs native Sphere.
+
+Compares epsilon grids at z=0 to verify the mesh geometry converges
+to the analytic sphere. Also tests eps_averaging by comparing the
+internal Yee grid epsilon (via sim.get_epsilon()) for a single resolution.
+Runs in ~1 second.
+"""
+import meep as mp
+import numpy as np
+import math
+import time
+
+
+def make_icosphere(radius, subdivisions):
+    phi = (1 + math.sqrt(5)) / 2
+    raw = [
+        (-1, phi, 0), (1, phi, 0), (-1, -phi, 0), (1, -phi, 0),
+        (0, -1, phi), (0, 1, phi), (0, -1, -phi), (0, 1, -phi),
+        (phi, 0, -1), (phi, 0, 1), (-phi, 0, -1), (-phi, 0, 1),
+    ]
+    verts = []
+    for v in raw:
+        n = math.sqrt(v[0]**2 + v[1]**2 + v[2]**2)
+        verts.append((v[0]/n*radius, v[1]/n*radius, v[2]/n*radius))
+    tris = [
+        (0,11,5),(0,5,1),(0,1,7),(0,7,10),(0,10,11),
+        (1,5,9),(5,11,4),(11,10,2),(10,7,6),(7,1,8),
+        (3,9,4),(3,4,2),(3,2,6),(3,6,8),(3,8,9),
+        (4,9,5),(2,4,11),(6,2,10),(8,6,7),(9,8,1),
+    ]
+    for _ in range(subdivisions):
+        edge_mid = {}
+        new_tris = []
+        for i0, i1, i2 in tris:
+            mids = []
+            for ea, eb in [(i0,i1),(i1,i2),(i2,i0)]:
+                edge = (min(ea,eb), max(ea,eb))
+                if edge not in edge_mid:
+                    va, vb = verts[ea], verts[eb]
+                    mx, my, mz = (va[0]+vb[0])/2, (va[1]+vb[1])/2, (va[2]+vb[2])/2
+                    n = math.sqrt(mx*mx+my*my+mz*mz)
+                    edge_mid[edge] = len(verts)
+                    verts.append((mx/n*radius, my/n*radius, mz/n*radius))
+                mids.append(edge_mid[edge])
+            m01, m12, m20 = mids
+            new_tris.extend([(i0,m01,m20),(i1,m12,m01),(i2,m20,m12),(m01,m12,m20)])
+        tris = new_tris
+    return verts, tris
+
+
+t0 = time.time()
+radius = 0.4
+cell = mp.Vector3(1.5, 1.5, 1.5)
+res = 4
+n_sample = 50
+mat = mp.Medium(epsilon=12)
+
+# --- Test 1: Geometric convergence (eps_averaging=False) ---
+print("=== Test 1: Geometric convergence ===", flush=True)
+
+# Reference: native sphere
+sim = mp.Simulation(cell_size=cell, geometry=[mp.Sphere(radius=radius, material=mat)],
+                    resolution=res, dimensions=3, eps_averaging=False)
+sim.init_sim()
+xtics = np.linspace(-cell.x/2, cell.x/2, n_sample)
+ytics = np.linspace(-cell.y/2, cell.y/2, n_sample)
+ztics = np.array([0.0])
+eps_sphere = np.real(sim.get_epsilon_grid(xtics, ytics, ztics)).flatten()
+
+errors = []
+for subdivs in [0, 2, 4]:
+    verts, tris = make_icosphere(radius, subdivs)
+    mesh_obj = mp.Mesh(vertices=verts, triangles=tris, material=mat)
+    sim = mp.Simulation(cell_size=cell, geometry=[mesh_obj],
+                        resolution=res, dimensions=3, eps_averaging=False)
+    sim.init_sim()
+    eps_mesh = np.real(sim.get_epsilon_grid(xtics, ytics, ztics)).flatten()
+    l2 = np.sqrt(np.mean((eps_mesh - eps_sphere)**2))
+    mismatch = np.sum(np.abs(eps_mesh - eps_sphere) > 0.5)
+    errors.append(l2)
+    print(f"  subdivs={subdivs}: {len(verts)} verts, {len(tris)} faces, "
+          f"L2={l2:.4f}, mismatches={mismatch}", flush=True)
+
+assert errors[-1] < errors[0], "FAIL: error should decrease with subdivision"
+print("  PASS: error decreases with subdivision", flush=True)
+
+# --- Test 2: eps_averaging produces non-binary epsilon ---
+print("\n=== Test 2: eps_averaging produces smoothed epsilon ===", flush=True)
+
+verts, tris = make_icosphere(radius, 2)
+mesh_obj = mp.Mesh(vertices=verts, triangles=tris, material=mat)
+
+# Without averaging
+sim_off = mp.Simulation(cell_size=cell, geometry=[mesh_obj],
+                        resolution=res, 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))
+unique_off = len(np.unique(np.round(eps_off, 2)))
+
+# With averaging
+sim_on = mp.Simulation(cell_size=cell, geometry=[mesh_obj],
+                       resolution=res, 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_on = len(np.unique(np.round(eps_on, 2)))
+
+print(f"  eps_averaging=False: {unique_off} unique epsilon values", flush=True)
+print(f"  eps_averaging=True:  {unique_on} unique epsilon values", flush=True)
+assert unique_on > unique_off, "FAIL: smoothing should produce more unique epsilon values"
+print("  PASS: smoothing produces intermediate epsilon values at interfaces", flush=True)
+
+dt = time.time() - t0
+print(f"\nAll tests passed in {dt:.1f}s", flush=True)