Add sphere convergence test: mesh icosphere vs native Sphere

Verifies that a mesh icosphere converges to meep's native Sphere
primitive as subdivision increases. L2 error decreases monotonically
at ~2x per subdivision level, with mismatched pixels dropping from
2276 (20 faces) to 12 (5120 faces).

Author: Luochenghuang

utils/test_sphere_convergence.py

@@ -0,0 +1,127 @@
+"""Convergence test: mesh icosphere vs native Sphere.
+
+Compares the epsilon grid of a mesh icosphere against meep's native
+Sphere at increasing resolution. Verifies that the mesh converges
+to the same geometry as the analytic primitive.
+"""
+import meep as mp
+import numpy as np
+import math
+
+def make_icosphere(radius, subdivisions, center=mp.Vector3(0, 0, 0)):
+    """Create an icosphere mesh by subdividing an icosahedron."""
+    # Golden ratio
+    phi = (1 + math.sqrt(5)) / 2
+
+    # Icosahedron vertices (normalized to unit sphere, then scaled)
+    raw_verts = [
+        (-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_verts:
+        norm = math.sqrt(v[0]**2 + v[1]**2 + v[2]**2)
+        verts.append((v[0]/norm * radius, v[1]/norm * radius, v[2]/norm * 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),
+    ]
+
+    # Subdivide
+    for _ in range(subdivisions):
+        edge_midpoints = {}
+        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_midpoints:
+                    va, vb = verts[ea], verts[eb]
+                    mx = (va[0]+vb[0])/2
+                    my = (va[1]+vb[1])/2
+                    mz = (va[2]+vb[2])/2
+                    norm = math.sqrt(mx*mx + my*my + mz*mz)
+                    edge_midpoints[edge] = len(verts)
+                    verts.append((mx/norm*radius, my/norm*radius, mz/norm*radius))
+                mids.append(edge_midpoints[edge])
+            m01, m12, m20 = mids
+            new_tris.extend([
+                (i0, m01, m20), (i1, m12, m01),
+                (i2, m20, m12), (m01, m12, m20),
+            ])
+        tris = new_tris
+
+    return mp.Mesh(
+        vertices=[mp.Vector3(*v) for v in verts],
+        triangles=tris,
+        center=center,
+        material=mp.Medium(epsilon=12),
+    )
+
+
+def get_eps_slice(geometry, cell_size, resolution, n_sample=100):
+    """Get a 2D epsilon slice at z=0."""
+    sim = mp.Simulation(
+        cell_size=cell_size,
+        geometry=geometry,
+        resolution=resolution,
+        dimensions=3,
+        eps_averaging=False,
+    )
+    sim.init_sim()
+    xtics = np.linspace(-cell_size.x/2, cell_size.x/2, n_sample)
+    ytics = np.linspace(-cell_size.y/2, cell_size.y/2, n_sample)
+    ztics = np.array([0.0])
+    eps = np.real(sim.get_epsilon_grid(xtics, ytics, ztics)).reshape(n_sample, n_sample)
+    return eps
+
+
+# Parameters
+radius = 0.4
+cell_size = mp.Vector3(1.5, 1.5, 1.5)
+n_sample = 200
+
+# Reference: native sphere
+print("=== Native Sphere (reference) ===", flush=True)
+sphere_geom = [mp.Sphere(radius=radius, center=mp.Vector3(0,0,0),
+                          material=mp.Medium(epsilon=12))]
+eps_sphere = get_eps_slice(sphere_geom, cell_size, resolution=4, n_sample=n_sample)
+
+# Mesh spheres at increasing subdivision
+print("\n=== Mesh Icosphere convergence ===", flush=True)
+print(f"{'subdivs':>8} {'verts':>8} {'faces':>8} {'L2 error':>12} {'max error':>12}", flush=True)
+
+errors = []
+for subdivs in range(5):
+    mesh_sphere = make_icosphere(radius, subdivs)
+    nv = len(mesh_sphere.vertices)
+    nf = len(mesh_sphere.triangles)
+
+    eps_mesh = get_eps_slice([mesh_sphere], cell_size, resolution=4, n_sample=n_sample)
+
+    # Compare against native sphere
+    diff = np.abs(eps_mesh - eps_sphere)
+    l2_err = np.sqrt(np.mean(diff**2))
+    max_err = np.max(diff)
+    mismatch = np.sum(np.abs(eps_mesh - eps_sphere) > 0.5)
+
+    print(f"{subdivs:>8} {nv:>8} {nf:>8} {l2_err:>12.6f} {max_err:>12.6f}  ({mismatch} mismatched pixels)", flush=True)
+    errors.append((nv, nf, l2_err, mismatch))
+
+# Check convergence: error should decrease with subdivision
+print("\n=== Convergence check ===", flush=True)
+for i in range(1, len(errors)):
+    ratio = errors[i-1][2] / errors[i][2] if errors[i][2] > 0 else float('inf')
+    print(f"  subdivs {i-1}->{i}: L2 error ratio = {ratio:.2f}x", flush=True)
+
+# Check: error should decrease monotonically and mismatched pixels should shrink
+monotonic = all(errors[i][2] < errors[i-1][2] for i in range(1, len(errors)))
+final_mismatches = errors[-1][3]
+if monotonic and final_mismatches < 50:
+    print(f"\nPASS: L2 error converging monotonically, {final_mismatches} mismatched pixels at highest subdivision", flush=True)
+else:
+    print(f"\nFAIL: monotonic={monotonic}, final mismatches={final_mismatches}", flush=True)