fix: Bug in PolySlab's intersections_plane method in cases where the plane is coincident with PolySlab side faces

Author: dmarek-flex

CHANGELOG.md

@@ -17,6 +17,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Fixed `AutoImpedanceSpec` validation to check path intersections against all conductors, not just filtered ones, as well as the mode plane bounds.
 - Fixed adjoint gradients being treated as zero due to scale-dependent `np.allclose(..., atol=1e-8)` checks, which could skip adjoint simulations and return zero gradients.
 - Fixed interpolation handling for permittivity and conductivity gradients in CustomMedium.
+- Fixed `intersections_plane` method in `PolySlab`, which sometimes missed vertices for planes coincident with `PolySlab` side faces.
 
 ## [2.10.0] - 2025-12-18
 

tests/test_components/test_geometry.py

@@ -941,6 +941,75 @@ def test_polyslab_intersection_inf_bounds():
     assert poly.intersections_plane(x=0)[0] == shapely.box(-1, -LARGE_NUMBER, 1, 0)
 
 
+def test_polyslab_intersection_with_coincident_plane():
+    """Test if intersection returns the correct shape when the plane is coincident with the side face."""
+    poly = td.PolySlab(
+        vertices=[[500.0, -7500.0], [500.0, 7500.0], [-500.0, 7500.0], [-500.0, -7500.0]],
+        slab_bounds=[0, 50],
+        axis=2,
+    )
+    # Each case should give one side face of the polyslab
+    expected_x_face = shapely.box(-7500, 0, 7500, 50)  # y-extent × z-extent
+    expected_y_face = shapely.box(-500, 0, 500, 50)  # x-extent × z-extent
+
+    assert poly.intersections_plane(x=-500) == [expected_x_face]
+    assert poly.intersections_plane(x=500) == [expected_x_face]
+    assert poly.intersections_plane(y=-7500) == [expected_y_face]
+    assert poly.intersections_plane(y=7500) == [expected_y_face]
+
+
+def test_polyslab_intersection_rotated_square():
+    """Test PolySlab plane intersection with a rotated square (diamond shape)."""
+    # Create a diamond by rotating a square 45 degrees
+    size = 2.0
+    angle = np.pi / 4
+    base_vertices = np.array(
+        [[-size / 2, -size / 2], [size / 2, -size / 2], [size / 2, size / 2], [-size / 2, size / 2]]
+    )
+    cos_a, sin_a = np.cos(angle), np.sin(angle)
+    rotation = np.array([[cos_a, -sin_a], [sin_a, cos_a]])
+    rotated = base_vertices @ rotation.T
+    rotated = rotated - rotated.min(axis=0) + 0.5  # shift to positive quadrant
+    vertices = [tuple(v) for v in rotated]
+
+    polyslab = td.PolySlab(vertices=vertices, slab_bounds=(0, 3), axis=2)
+
+    all_verts = np.array(vertices)
+    left_tip_x = all_verts[:, 0].min()
+    bottom_tip_y = all_verts[:, 1].min()
+    x_center = (all_verts[:, 0].min() + all_verts[:, 0].max()) / 2
+
+    # Test 1: Cut at z=1.5 (middle of slab) - should give full diamond
+    cross_section = polyslab.intersections_plane(z=1.5)
+    assert len(cross_section) == 1
+    assert np.isclose(cross_section[0].area, 4.0)
+
+    # Test 2: Cut through center at x=x_center - should give rectangle
+    cross_section = polyslab.intersections_plane(x=x_center)
+    assert len(cross_section) == 1
+    assert cross_section[0].area > 0
+
+    # Test 3: Cut at left corner tip (tangent touch) - should give degenerate shape
+    cross_section = polyslab.intersections_plane(x=left_tip_x)
+    assert len(cross_section) == 1
+    assert np.isclose(cross_section[0].area, 0.0)
+
+    # Test 4: Cut near left corner (slightly inside) - should give small shape
+    cross_section = polyslab.intersections_plane(x=left_tip_x + 0.3)
+    assert len(cross_section) == 1
+    assert cross_section[0].area > 0
+
+    # Test 5: Cut at bottom corner (tangent touch) - should give degenerate shape
+    cross_section = polyslab.intersections_plane(y=bottom_tip_y)
+    assert len(cross_section) == 1
+    assert np.isclose(cross_section[0].area, 0.0)
+
+    # Test 6: Cut at z=0 (bottom boundary) - should give full diamond
+    cross_section = polyslab.intersections_plane(z=0)
+    assert len(cross_section) == 1
+    assert np.isclose(cross_section[0].area, 4.0)
+
+
 def test_from_shapely():
     ring = shapely.LinearRing([(-16, 9), (-8, 9), (-12, 2)])
     poly = shapely.Polygon([(-2, 0), (-10, 0), (-6, 7)])

tidy3d/components/geometry/polyslab.py

@@ -858,23 +858,22 @@ def _find_intersecting_ys_angle_vertical(
         x_vertices_f, _ = vertices_f.T
         x_vertices_axis, _ = vertices_axis.T
 
-        # find which segments intersect
-        f_left_to_intersect = x_vertices_f <= position
-        orig_right_to_intersect = x_vertices_axis > position
-        intersects_b = np.logical_and(f_left_to_intersect, orig_right_to_intersect)
+        # Find which segments intersect:
+        # 1. Strictly crossing: one endpoint strictly left, one strictly right
+        # 2. Touching: exactly one endpoint on the plane (xor), which excludes
+        #    edges lying entirely on the plane (both endpoints at position).
+        orig_on_plane = np.isclose(x_vertices_axis, position, rtol=_IS_CLOSE_RTOL)
+        f_on_plane = np.roll(orig_on_plane, shift=-1)
+        crosses_b = (x_vertices_axis > position) & (x_vertices_f < position)
+        crosses_f = (x_vertices_axis < position) & (x_vertices_f > position)
 
-        f_right_to_intersect = x_vertices_f > position
-        orig_left_to_intersect = x_vertices_axis <= position
-        intersects_f = np.logical_and(f_right_to_intersect, orig_left_to_intersect)
-
-        # exclude vertices at the position if exclude_on_vertices is True
         if exclude_on_vertices:
-            intersects_on = np.isclose(x_vertices_axis, position, rtol=_IS_CLOSE_RTOL)
-            intersects_f_on = np.isclose(x_vertices_f, position, rtol=_IS_CLOSE_RTOL)
-            intersects_both_off = np.logical_not(np.logical_or(intersects_on, intersects_f_on))
-            intersects_f &= intersects_both_off
-            intersects_b &= intersects_both_off
-        intersects_segment = np.logical_or(intersects_b, intersects_f)
+            # exclude vertices at the position
+            not_touching = np.logical_not(orig_on_plane | f_on_plane)
+            intersects_segment = (crosses_b | crosses_f) & not_touching
+        else:
+            single_touch = np.logical_xor(orig_on_plane, f_on_plane)
+            intersects_segment = crosses_b | crosses_f | single_touch
 
         iverts_b = vertices_axis[intersects_segment]
         iverts_f = vertices_f[intersects_segment]
@@ -893,10 +892,27 @@ def _find_intersecting_ys_angle_vertical(
         ints_y = np.array(ints_y)
         ints_angle = np.array(ints_angle)
 
-        sort_index = np.argsort(ints_y)
-        ints_y_sort = ints_y[sort_index]
+        # Get rid of duplicate intersection points (vertices counted twice if directly on position)
+        ints_y_sort, sort_index = np.unique(ints_y, return_index=True)
         ints_angle_sort = ints_angle[sort_index]
 
+        # For tangent touches (vertex on plane, both neighbors on same side),
+        # add y-value back to form a degenerate pair
+        if not exclude_on_vertices:
+            n = len(vertices_axis)
+            for idx in np.where(orig_on_plane)[0]:
+                prev_on = orig_on_plane[(idx - 1) % n]
+                next_on = orig_on_plane[(idx + 1) % n]
+                if not prev_on and not next_on:
+                    prev_side = x_vertices_axis[(idx - 1) % n] > position
+                    next_side = x_vertices_axis[(idx + 1) % n] > position
+                    if prev_side == next_side:
+                        ints_y_sort = np.append(ints_y_sort, vertices_axis[idx, 1])
+                        ints_angle_sort = np.append(ints_angle_sort, 0)
+
+        sort_index = np.argsort(ints_y_sort)
+        ints_y_sort = ints_y_sort[sort_index]
+        ints_angle_sort = ints_angle[sort_index]
         return ints_y_sort, ints_angle_sort
 
     def _find_intersecting_ys_angle_slant(