Slanted prism (sidewall_angle != 0) with a collinear vertex: wrong geometry, or an abort with subpixel averaging #100

I was working on some simulations with meep and encountered a bug with libctl not creating some geometries correctly. I had AI debug it and it left the following report. I am not familiar enough with the codebase to know whether this is completely accurate, but I wanted to leave an issue in case someone found it useful.

Summary

If a prism with a nonzero sidewall_angle has a vertex that lies exactly on the line between its two neighbours, point_in_prism gives wrong answers for every height z > 0. In the case below, the prism's interior is reported as outside, and a region outside the polygon is reported as inside. With subpixel averaging off, Meep then renders only a thin wedge of the prism. With subpixel averaging on (Meep's default), it doesn't get that far: init_sim aborts with degenerate plane in project_point_into_plane. With sidewall_angle = 0 the same vertex list works correctly.

Collinear vertices are common in practice. They come from layout tools, from merging polygons, and from deliberately subdividing long edges, which is the obvious workaround for the long-edge halo fixed in #99 until that fix is released.

Reproducer (Meep 1.34.0)

import math
import numpy as np
import meep as mp

mp.verbosity(0)

H, ANGLE = 0.5, math.radians(30)   # 30 degree sidewall, tapering inward
Z = 0.25                           # probe height; expected half-widths shrink by Z*tan(30) = 0.144

def inside_intervals(prism, along, lo=-3.0, hi=3.0, n=6001):
    """Intervals along the x or y axis (through the origin, at height Z) where is_point_in_object is True."""
    ts = np.linspace(lo, hi, n)
    pts = [mp.Vector3(t, 0, Z) if along == "x" else mp.Vector3(0, t, Z) for t in ts]
    inside = np.array([mp.is_point_in_object(p, prism) for p in pts])
    edges = np.flatnonzero(np.diff(inside.astype(int)))
    bounds = [lo] * int(inside[0]) + [round(float(ts[e + 1]), 3) for e in edges] + [hi] * int(inside[-1])
    return list(zip(bounds[::2], bounds[1::2]))

def rendered_extent(prism):
    """x-range where meep actually puts the prism material, along y=0, z=Z."""
    sim = mp.Simulation(cell_size=mp.Vector3(6, 3, 1), resolution=50, geometry=[prism], eps_averaging=False)
    sim.init_sim()
    eps = sim.get_array(center=mp.Vector3(0, 0, Z), size=mp.Vector3(6, 0, 0), component=mp.Dielectric)
    x = np.linspace(-3, 3, len(eps))[eps > 1.5]
    return (round(float(x.min()), 2), round(float(x.max()), 2)) if len(x) else None

# The same 4 x 2 rectangle, with and without one extra vertex at the midpoint of its x = +2 edge.
rectangle = [(-2, -1), (2, -1), (2, 1), (-2, 1)]
rectangle_with_collinear_vertex = [(-2, -1), (2, -1), (2, 0), (2, 1), (-2, 1)]

print(f"expected at z={Z}: x in [-2, 2], y in [-1, 1] for sidewall_angle=0; x in [-1.856, 1.856], y in [-0.856, 0.856] for 30 deg")
for name, xy in [("rectangle", rectangle), ("rectangle + collinear vertex", rectangle_with_collinear_vertex)]:
    for angle in (0.0, ANGLE):
        prism = mp.Prism([mp.Vector3(*v) for v in xy], height=H, sidewall_angle=angle, material=mp.Medium(index=2))
        print(f"{name}, sidewall_angle={math.degrees(angle):.0f} deg:\n"
              f"    is_point_in_object inside along x: {inside_intervals(prism, 'x')}\n"
              f"    is_point_in_object inside along y: {inside_intervals(prism, 'y')}\n"
              f"    rendered epsilon along x:          {rendered_extent(prism)}")

Output

expected at z=0.25: x in [-2, 2], y in [-1, 1] for sidewall_angle=0; x in [-1.856, 1.856], y in [-0.856, 0.856] for 30 deg
rectangle, sidewall_angle=0 deg:
    is_point_in_object inside along x: [(-2.002, 2.003)]
    is_point_in_object inside along y: [(-1.004, 1.005)]
    rendered epsilon along x:          (-2.0, 2.0)
rectangle, sidewall_angle=30 deg:
    is_point_in_object inside along x: [(-1.857, 1.858)]
    is_point_in_object inside along y: [(-0.859, 0.86)]
    rendered epsilon along x:          (-1.86, 1.86)
rectangle + collinear vertex, sidewall_angle=0 deg:
    is_point_in_object inside along x: [(-2.002, 2.001)]
    is_point_in_object inside along y: [(-1.004, 1.005)]
    rendered epsilon along x:          (-2.0, 2.0)
rectangle + collinear vertex, sidewall_angle=30 deg:
    is_point_in_object inside along x: [(-3.0, -1.853)]
    is_point_in_object inside along y: [(-0.859, -0.851), (0.852, 0.86)]
    rendered epsilon along x:          (-2.0, -1.86)

In the last case, everything to the left of the left sidewall counts as inside, out to the edge of the sampled range, and the interior counts as outside. The rendered ε is only the sliver between the bounding box (x = −2) and the slanted left wall. The small boundary slop of about 0.004 in the other cases is the edge-length halo fixed in #99.

With subpixel averaging, the same prism aborts instead:

import math
import meep as mp

vertices = [(-2, -1), (2, -1), (2, 0), (2, 1), (-2, 1)]   # rectangle + collinear vertex
prism = mp.Prism([mp.Vector3(*v) for v in vertices], height=0.5, sidewall_angle=math.radians(30),
                 material=mp.Medium(index=2))
sim = mp.Simulation(cell_size=mp.Vector3(6, 3, 1), resolution=20, geometry=[prism], eps_averaging=True)
sim.init_sim()

degenerate plane in project_point_into_plane

The process exits with status 1. Without the vertex at (2, 0), init_sim completes.

Cause

When sidewall_angle != 0, init_prism builds each top vertex by intersecting the two adjacent offset edges with the two-line intersection formula (utils/geom.c#L3985-L3986; geom.c is unchanged on master since that commit, the merge of #99):

number px = ((x1*y2-y1*x2)*(x3-x4)-(x1-x2)*(x3*y4-y3*x4)) / ((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4));
number py = ((x1*y2-y1*x2)*(y3-y4)-(y1-y2)*(x3*y4-y3*x4)) / ((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4));

There's no guard for parallel edges. When vertex nv is collinear with its neighbours, offset edges nv and nv+1 lie on the same line, both the numerator and the denominator are zero, and the top vertex becomes NaN. point_in_or_on_prism then interpolates every query polygon using that NaN vertex, which breaks the even-odd crossing count in node_in_or_on_polygon. Subpixel averaging goes through normal_to_prism, which builds the sidewall quadrilaterals from the same NaN vertex; their normal is NaN, and the CHECK(vector3_norm(v3) > 1.0e-6, ...) in normal_distance_to_plane fails (utils/geom.c#L3471). With sidewall_angle == 0 this block is skipped, which is why the untilted case is fine.

Suggested fix

When the adjacent edges are parallel, the offset lines coincide. The correct top vertex is then just the original vertex shifted by the shared offset, which is already stored as top_polygon_edges[nv].a2:

number denom = (x1-x2)*(y3-y4)-(y1-y2)*(x3-x4);
if (fabs(denom) <= 1e-12 * vector3_norm(top_polygon_edges[nv].v)
                          * vector3_norm(top_polygon_edges[(nv + 1 == num_vertices ? 0 : nv + 1)].v)) {
  // adjacent edges parallel: the offset lines coincide, so keep the offset vertex
  prsm->vertices_top_p.items[nv].x = top_polygon_edges[nv].a2.x;
  prsm->vertices_top_p.items[nv].y = top_polygon_edges[nv].a2.y;
}
else { /* existing px, py computation using denom */ }

denom is the cross product of the two edge vectors, so the test is on the sine of the angle between them. The relative tolerance also catches vertices that are collinear only up to rounding, as they typically are after a layout tool's integer database units are scaled to floats. An alternative is to drop collinear vertices in init_prism before building the top polygon.

A regression test would fit next to test_slanted_prism_cross_section from #99. It would build the rectangle above with and without the midpoint vertex and check that the two prisms give the same results.

Related

• Bug fix for prism cross section #99 / Prism side-wall angle does not give expected results. meep#3311: the long-edge halo. It's fixed on master but not in a release yet; v4.7.1 is the latest. Subdividing long edges avoids the halo on 4.7.1, but it creates exactly the collinear vertices that trigger this bug. Our workaround is to push the inserted vertices off their edge by 1e-6, alternately outward and inward, so that no three vertices are exactly collinear (shifting them all the same way would leave them on one line again).
• Artifacts around slanted prism objects #55 reports artifacts on slanted prisms too, but its hexagon has no collinear vertices, so it is a different problem.

Environment

• libctl 4.7.1 (conda-forge hb03c661_0); the relevant code is unchanged on master at 1e7b4c1
• Meep 1.34.0 (conda-forge, mpi_mpich build), Python 3.13.15, NumPy 2.5.3
• Linux x86_64 (RHEL 9)