Fix PolySlab intersections when bounds contain inf

Author: weiliangjin2021

CHANGELOG.md

@@ -18,6 +18,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Fixed
 - `tidy3d.plugins.design.Results` store the `BatchData` for batch runs in the `.batch_data` field.
 - Prompting user to check solver log when loading solver data if warnings were found in the log, or if the simulation diverged or errored.
+- Bug in PolySlab intersection if slab bounds are `inf` on one side.
 
 ### Changed
 - Slightly reorganized `web.run` logging when `verbose=True` to make it clearer.

tests/test_components/test_geometry.py

@@ -11,6 +11,7 @@
 import warnings
 
 import tidy3d as td
+from tidy3d.constants import LARGE_NUMBER
 from tidy3d.exceptions import SetupError, Tidy3dKeyError, ValidationError
 from tidy3d.components.geometry.base import Planar
 from tidy3d.components.geometry.utils import flatten_groups, traverse_geometries
@@ -710,6 +711,22 @@ def test_polyslab_axis(axis):
     assert not ps.intersects_plane(x=plane_coord[0], y=plane_coord[1], z=plane_coord[2])
 
 
+def test_polyslab_intersection_inf_bounds():
+    """Test if intersection returns correct shapes when one of the slab_bounds is Inf."""
+    # 1) [0, inf]
+    poly = td.PolySlab(
+        vertices=[[2, -1], [-2, -1], [-2, 1], [2, 1]],
+        slab_bounds=[0, td.inf],
+    )
+    assert len(poly.intersections_plane(x=0)) == 1
+    assert poly.intersections_plane(x=0)[0] == shapely.box(-1, 0.0, 1, LARGE_NUMBER)
+
+    # 2) [-inf, 0]
+    poly = poly.updated_copy(slab_bounds=[-td.inf, 0])
+    assert len(poly.intersections_plane(x=0)) == 1
+    assert poly.intersections_plane(x=0)[0] == shapely.box(-1, -LARGE_NUMBER, 1, 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

@@ -16,7 +16,7 @@
 from ..types import MatrixReal4x4, Shapely
 from ...log import log
 from ...exceptions import SetupError, ValidationError
-from ...constants import MICROMETER, fp_eps
+from ...constants import MICROMETER, fp_eps, LARGE_NUMBER
 from ...packaging import verify_packages_import
 
 from . import base
@@ -353,6 +353,8 @@ def center_axis(self) -> float:
         zmin, zmax = self.slab_bounds
         if np.isneginf(zmin) and np.isposinf(zmax):
             return 0.0
+        zmin = max(zmin, -LARGE_NUMBER)
+        zmax = min(zmax, LARGE_NUMBER)
         return (zmax + zmin) / 2.0
 
     @property
@@ -386,7 +388,7 @@ def middle_polygon(self) -> np.ndarray:
             The vertices of the polygon at the middle.
         """
 
-        dist = self._extrusion_length_to_offset_distance(self.length_axis / 2)
+        dist = self._extrusion_length_to_offset_distance(self.finite_length_axis / 2)
         if self.reference_plane == "bottom":
             return self._shift_vertices(self.reference_polygon, dist)[0]
         if self.reference_plane == "top":
@@ -405,7 +407,7 @@ def base_polygon(self) -> np.ndarray:
         """
         if self.reference_plane == "bottom":
             return self.reference_polygon
-        dist = self._extrusion_length_to_offset_distance(-self.length_axis / 2)
+        dist = self._extrusion_length_to_offset_distance(-self.finite_length_axis / 2)
         return self._shift_vertices(self.middle_polygon, dist)[0]
 
     @cached_property
@@ -419,7 +421,7 @@ def top_polygon(self) -> np.ndarray:
         """
         if self.reference_plane == "top":
             return self.reference_polygon
-        dist = self._extrusion_length_to_offset_distance(self.length_axis / 2)
+        dist = self._extrusion_length_to_offset_distance(self.finite_length_axis / 2)
         return self._shift_vertices(self.middle_polygon, dist)[0]
 
     @cached_property
@@ -461,7 +463,7 @@ def inside(
 
         z0 = self.center_axis
         dist_z = np.abs(z - z0)
-        inside_height = dist_z <= (self.length_axis / 2)
+        inside_height = dist_z <= (self.finite_length_axis / 2)
 
         # avoid going into face checking if no points are inside slab bounds
         if not np.any(inside_height):
@@ -646,14 +648,14 @@ def _intersections_side(self, position, axis) -> list:
 
         # find out all z_i where the plane will intersect the vertex
         z0 = self.center_axis
-        z_base = z0 - self.length_axis / 2
+        z_base = z0 - self.finite_length_axis / 2
 
         axis_ordered = self._order_axis(axis)
         height_list = self._find_intersecting_height(position, axis_ordered)
         polys = []
 
         # looping through z_i to assemble the polygons
-        height_list = np.append(height_list, self.length_axis)
+        height_list = np.append(height_list, self.finite_length_axis)
         h_base = 0.0
         for h_top in height_list:
             # length within between top and bottom
@@ -670,7 +672,7 @@ def _intersections_side(self, position, axis) -> list:
                 )
             else:
                 # for slanted sidewall, move up by `fp_eps` in case vertices are degenerate at the base.
-                dist = -(h_base - self.length_axis / 2 + fp_eps) * self._tanq
+                dist = -(h_base - self.finite_length_axis / 2 + fp_eps) * self._tanq
                 vertices = self._shift_vertices(self.middle_polygon, dist)[0]
                 ints_y, ints_angle = self._find_intersecting_ys_angle_slant(
                     vertices, position, axis_ordered
@@ -744,14 +746,14 @@ def _find_intersecting_height(self, position: float, axis: int) -> np.ndarray:
 
         # distance to the plane in the direction of vertex shifting
         distance = self.middle_polygon[:, axis] - position
-        height = distance / self._tanq / shift_val + self.length_axis / 2
+        height = distance / self._tanq / shift_val + self.finite_length_axis / 2
         height = np.unique(height)
         # further filter very close ones
         is_not_too_close = np.insert((np.diff(height) > fp_eps), 0, True)
         height = height[is_not_too_close]
 
         height = height[height > fp_eps]
-        height = height[height < self.length_axis - fp_eps]
+        height = height[height < self.finite_length_axis - fp_eps]
         return height
 
     def _find_intersecting_ys_angle_vertical(
@@ -955,11 +957,11 @@ def bounds(self) -> Bound:
         max_offset = self.dilation
         if not isclose(self.sidewall_angle, 0):
             if self.reference_plane == "bottom":
-                max_offset += max(0, -self._tanq * self.length_axis)
+                max_offset += max(0, -self._tanq * self.finite_length_axis)
             elif self.reference_plane == "top":
-                max_offset += max(0, self._tanq * self.length_axis)
+                max_offset += max(0, self._tanq * self.finite_length_axis)
             elif self.reference_plane == "middle":
-                max_offset += max(0, abs(self._tanq) * self.length_axis / 2)
+                max_offset += max(0, abs(self._tanq) * self.finite_length_axis / 2)
 
         # special care when dilated
         if max_offset > 0:
@@ -1513,7 +1515,7 @@ def _dilation_length(self) -> List[float]:
         """dilation length from reference plane to the top/bottom of the polyslab."""
 
         # for "bottom", only needs to compute the offset length to the top
-        dist = [self._extrusion_length_to_offset_distance(self.length_axis)]
+        dist = [self._extrusion_length_to_offset_distance(self.finite_length_axis)]
         # reverse the dilation value if the reference plane is on the top
         if self.reference_plane == "top":
             dist = [-dist[0]]
@@ -1527,9 +1529,9 @@ def _dilation_value_at_reference_to_coord(self, dilation: float) -> float:
 
         z_coord = -dilation / self._tanq + self.slab_bounds[0]
         if self.reference_plane == "middle":
-            return z_coord + self.length_axis / 2
+            return z_coord + self.finite_length_axis / 2
         if self.reference_plane == "top":
-            return z_coord + self.length_axis
+            return z_coord + self.finite_length_axis
         # bottom case
         return z_coord