diff --git a/CHANGELOG.md b/CHANGELOG.md index 10d390105..daa533145 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/tests/test_components/test_geometry.py b/tests/test_components/test_geometry.py index c08a2b022..8ede268e9 100644 --- a/tests/test_components/test_geometry.py +++ b/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)]) diff --git a/tidy3d/components/geometry/polyslab.py b/tidy3d/components/geometry/polyslab.py index bf3f141cf..ccd9e9ec5 100644 --- a/tidy3d/components/geometry/polyslab.py +++ b/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