From 0a1e13b6f200f53e1faab2c5beae612fc28b9b03 Mon Sep 17 00:00:00 2001 From: Weiliang Jin Date: Mon, 24 Nov 2025 08:49:27 -0800 Subject: [PATCH] feat(meshing): disable polygon cleanup in corner finder Corner finder: reduce the number of vertices for smoothly curved cross section Mesher and gap refinement: performance improvement in parse_structures --- tests/test_components/test_geometry.py | 29 ++++ tidy3d/components/geometry/base.py | 200 +++++++++++++++++++---- tidy3d/components/geometry/mesh.py | 24 ++- tidy3d/components/geometry/polyslab.py | 25 ++- tidy3d/components/geometry/primitives.py | 63 +++++-- tidy3d/components/geometry/utils.py | 9 +- tidy3d/components/grid/corner_finder.py | 11 +- tidy3d/components/grid/grid_spec.py | 4 +- tidy3d/components/grid/mesher.py | 102 +++++++++--- 9 files changed, 397 insertions(+), 70 deletions(-) diff --git a/tests/test_components/test_geometry.py b/tests/test_components/test_geometry.py index d00a11515b..ca444807c8 100644 --- a/tests/test_components/test_geometry.py +++ b/tests/test_components/test_geometry.py @@ -236,6 +236,35 @@ def test_intersections_plane_inf(): assert len(c.intersections_plane(y=0)) == 1 +@pytest.mark.parametrize("component", [BOX, CYLINDER, SPHERE, POLYSLAB, UNION, GROUP]) +@pytest.mark.parametrize("cleanup", [True, False]) +def test_intersections_plane_cleanup_param(component, cleanup): + """Test that cleanup parameter is accepted by all geometry types.""" + shapes = component.intersections_plane(z=0, cleanup=cleanup) + assert isinstance(shapes, list) + shapes_default = component.intersections_plane(z=0) + # Both should return valid shapely objects + for shape in shapes: + assert hasattr(shape, "is_valid") + for shape in shapes_default: + assert hasattr(shape, "is_valid") + + +@pytest.mark.parametrize("component", [SPHERE, CYLINDER]) +@pytest.mark.parametrize("quad_segs", [8, 50, 200]) +def test_intersections_plane_quad_segs(component, quad_segs): + """Test that quad_segs parameter controls discretization of circular shapes.""" + shapes = component.intersections_plane(z=0, quad_segs=quad_segs) + assert len(shapes) > 0 + # For circular shapes, quad_segs should affect the number of vertices + # Higher quad_segs should generally give more vertices + shape = shapes[0] + num_coords = len(shape.exterior.coords) + assert num_coords > 4 * quad_segs, ( + f"Expected more than {4 * quad_segs} coords, got {num_coords}" + ) + + def test_center_not_inf_validate(): with pytest.raises(pydantic.ValidationError): _ = td.Box(center=(td.inf, 0, 0)) diff --git a/tidy3d/components/geometry/base.py b/tidy3d/components/geometry/base.py index 41682ae736..a6eae24230 100644 --- a/tidy3d/components/geometry/base.py +++ b/tidy3d/components/geometry/base.py @@ -229,7 +229,12 @@ def inside_meshgrid( @abstractmethod def intersections_tilted_plane( - self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4 + self, + normal: Coordinate, + origin: Coordinate, + to_2D: MatrixReal4x4, + cleanup: bool = True, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Return a list of shapely geometries at the plane specified by normal and origin. @@ -241,6 +246,11 @@ def intersections_tilted_plane( Vector defining the plane origin. to_2D : MatrixReal4x4 Transformation matrix to apply to resulting shapes. + cleanup : bool = True + If True, removes extremely small features from each polygon's boundary. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. If ``None``, uses + high-quality visualization settings. Returns ------- @@ -251,7 +261,12 @@ def intersections_tilted_plane( """ def intersections_plane( - self, x: Optional[float] = None, y: Optional[float] = None, z: Optional[float] = None + self, + x: Optional[float] = None, + y: Optional[float] = None, + z: Optional[float] = None, + cleanup: bool = True, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Returns list of shapely geometries at plane specified by one non-None value of x,y,z. @@ -263,6 +278,11 @@ def intersections_plane( Position of plane in y direction, only one of x,y,z can be specified to define plane. z : float = None Position of plane in z direction, only one of x,y,z can be specified to define plane. + cleanup : bool = True + If True, removes extremely small features from each polygon's boundary. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. If ``None``, uses + high-quality visualization settings. Returns ------- @@ -278,7 +298,9 @@ def intersections_plane( if axis != 2: last, indices = self.pop_axis((0, 1, 2), axis) to_2D = to_2D[[*list(indices), last, 3]] - return self.intersections_tilted_plane(normal, origin, to_2D) + return self.intersections_tilted_plane( + normal, origin, to_2D, cleanup=cleanup, quad_segs=quad_segs + ) def intersections_2dbox(self, plane: Box) -> list[Shapely]: """Returns list of shapely geometries representing the intersections of the geometry with @@ -1564,7 +1586,12 @@ class SimplePlaneIntersection(Geometry, ABC): """A geometry where intersections with an axis aligned plane may be computed efficiently.""" def intersections_tilted_plane( - self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4 + self, + normal: Coordinate, + origin: Coordinate, + to_2D: MatrixReal4x4, + cleanup: bool = True, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Return a list of shapely geometries at the plane specified by normal and origin. Checks special cases before relying on the complete computation. @@ -1577,6 +1604,11 @@ def intersections_tilted_plane( Vector defining the plane origin. to_2D : MatrixReal4x4 Transformation matrix to apply to resulting shapes. + cleanup : bool = True + If True, removes extremely small features from each polygon's boundary. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. If ``None``, uses + high-quality visualization settings. Returns ------- @@ -1591,7 +1623,7 @@ def intersections_tilted_plane( axis = np.argmax(np.abs(normal)).item() coord = "xyz"[axis] kwargs = {coord: origin[axis]} - section = self.intersections_plane(**kwargs) + section = self.intersections_plane(cleanup=cleanup, quad_segs=quad_segs, **kwargs) # Apply transformation in the plane by removing row and column to_2D_in_plane = np.delete(np.delete(to_2D, 2, 0), axis, 1) @@ -1603,11 +1635,17 @@ def transform(p_array: NDArray) -> NDArray: transformed_section = shapely.transform(section, transformation=transform) return transformed_section # Otherwise compute the arbitrary intersection - return self._do_intersections_tilted_plane(normal=normal, origin=origin, to_2D=to_2D) + return self._do_intersections_tilted_plane( + normal=normal, origin=origin, to_2D=to_2D, quad_segs=quad_segs + ) @abstractmethod def _do_intersections_tilted_plane( - self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4 + self, + normal: Coordinate, + origin: Coordinate, + to_2D: MatrixReal4x4, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Return a list of shapely geometries at the plane specified by normal and origin. @@ -1619,6 +1657,8 @@ def _do_intersections_tilted_plane( Vector defining the plane origin. to_2D : MatrixReal4x4 Transformation matrix to apply to resulting shapes. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. Returns ------- @@ -1704,7 +1744,12 @@ def reference_axis_pos(self) -> float: return self.center_axis def intersections_plane( - self, x: Optional[float] = None, y: Optional[float] = None, z: Optional[float] = None + self, + x: Optional[float] = None, + y: Optional[float] = None, + z: Optional[float] = None, + cleanup: bool = True, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Returns shapely geometry at plane specified by one non None value of x,y,z. @@ -1716,29 +1761,37 @@ def intersections_plane( Position of plane in y direction, only one of x,y,z can be specified to define plane. z : float Position of plane in z direction, only one of x,y,z can be specified to define plane. + cleanup : bool = True + If True, removes extremely small features from each polygon's boundary. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. If ``None``, uses + high-quality visualization settings. Returns ------- List[shapely.geometry.base.BaseGeometry] List of 2D shapes that intersect plane. For more details refer to - `Shapely's Documentation `_. + `Shapely's Documentation ``. """ axis, position = self.parse_xyz_kwargs(x=x, y=y, z=z) if not self.intersects_axis_position(axis, position): return [] if axis == self.axis: - return self._intersections_normal(position) + return self._intersections_normal(position, quad_segs=quad_segs) return self._intersections_side(position, axis) @abstractmethod - def _intersections_normal(self, z: float) -> list: + def _intersections_normal(self, z: float, quad_segs: Optional[int] = None) -> list: """Find shapely geometries intersecting planar geometry with axis normal to slab. Parameters ---------- z : float Position along the axis normal to slab + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. If ``None``, uses + high-quality visualization settings. Returns ------- @@ -2021,7 +2074,11 @@ def surfaces_with_exclusion(cls, size: Size, center: Coordinate, **kwargs: Any) @verify_packages_import(["trimesh"]) def _do_intersections_tilted_plane( - self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4 + self, + normal: Coordinate, + origin: Coordinate, + to_2D: MatrixReal4x4, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Return a list of shapely geometries at the plane specified by normal and origin. @@ -2033,6 +2090,8 @@ def _do_intersections_tilted_plane( Vector defining the plane origin. to_2D : MatrixReal4x4 Transformation matrix to apply to resulting shapes. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. Not used for Box geometry. Returns ------- @@ -2071,7 +2130,12 @@ def _do_intersections_tilted_plane( return path.polygons_full def intersections_plane( - self, x: Optional[float] = None, y: Optional[float] = None, z: Optional[float] = None + self, + x: Optional[float] = None, + y: Optional[float] = None, + z: Optional[float] = None, + cleanup: bool = True, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Returns shapely geometry at plane specified by one non None value of x,y,z. @@ -2083,6 +2147,10 @@ def intersections_plane( Position of plane in y direction, only one of x,y,z can be specified to define plane. z : float = None Position of plane in z direction, only one of x,y,z can be specified to define plane. + cleanup : bool = True + If True, removes extremely small features from each polygon's boundary. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. Not used for Box geometry. Returns ------- @@ -2138,10 +2206,22 @@ def inside(self, x: NDArray[float], y: NDArray[float], z: NDArray[float]) -> NDA dist_z = np.abs(z - z0) return (dist_x <= Lx / 2) * (dist_y <= Ly / 2) * (dist_z <= Lz / 2) - def intersections_with(self, other: Shapely) -> list[Shapely]: + def intersections_with( + self, other: Shapely, cleanup: bool = True, quad_segs: Optional[int] = None + ) -> list[Shapely]: """Returns list of shapely geometries representing the intersections of the geometry with this 2D box. + Parameters + ---------- + other : Shapely + Geometry to intersect with. + cleanup : bool = True + If True, removes extremely small features from each polygon's boundary. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. If ``None``, uses + high-quality visualization settings. + Returns ------- List[shapely.geometry.base.BaseGeometry] @@ -2165,7 +2245,7 @@ def intersections_with(self, other: Shapely) -> list[Shapely]: dim = "xyz"[normal_ind] pos = self.center[normal_ind] xyz_kwargs = {dim: pos} - shapes_plane = other.intersections_plane(**xyz_kwargs) + shapes_plane = other.intersections_plane(cleanup=cleanup, quad_segs=quad_segs, **xyz_kwargs) # intersect all shapes with the input self bs_min, bs_max = (self.pop_axis(bounds, axis=normal_ind)[1] for bounds in self.bounds) @@ -2762,7 +2842,12 @@ def bounds(self) -> Bound: return (tuple(vertices.min(axis=1)), tuple(vertices.max(axis=1))) def intersections_tilted_plane( - self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4 + self, + normal: Coordinate, + origin: Coordinate, + to_2D: MatrixReal4x4, + cleanup: bool = True, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Return a list of shapely geometries at the plane specified by normal and origin. @@ -2774,6 +2859,11 @@ def intersections_tilted_plane( Vector defining the plane origin. to_2D : MatrixReal4x4 Transformation matrix to apply to resulting shapes. + cleanup : bool = True + If True, removes extremely small features from each polygon's boundary. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. If ``None``, uses + high-quality visualization settings. Returns ------- @@ -2786,6 +2876,8 @@ def intersections_tilted_plane( tuple(np.dot((normal[0], normal[1], normal[2], 0.0), self.transform)[:3]), tuple(np.dot(self.inverse, (origin[0], origin[1], origin[2], 1.0))[:3]), np.dot(to_2D, self.transform), + cleanup=cleanup, + quad_segs=quad_segs, ) def inside(self, x: NDArray[float], y: NDArray[float], z: NDArray[float]) -> NDArray[bool]: @@ -3026,7 +3118,9 @@ def to_polygon_list(base_geometry: Shapely, cleanup: bool = False) -> list[Shape unfiltered_geoms = [] if base_geometry.geom_type == "GeometryCollection": unfiltered_geoms = [ - p for geom in base_geometry.geoms for p in ClipOperation.to_polygon_list(geom) + p + for geom in base_geometry.geoms + for p in ClipOperation.to_polygon_list(geom, cleanup) ] if base_geometry.geom_type == "MultiPolygon": unfiltered_geoms = [p for p in base_geometry.geoms if not p.is_empty] @@ -3069,7 +3163,12 @@ def _bit_operation(self) -> Callable[[Any, Any], Any]: return result def intersections_tilted_plane( - self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4 + self, + normal: Coordinate, + origin: Coordinate, + to_2D: MatrixReal4x4, + cleanup: bool = True, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Return a list of shapely geometries at the plane specified by normal and origin. @@ -3081,6 +3180,11 @@ def intersections_tilted_plane( Vector defining the plane origin. to_2D : MatrixReal4x4 Transformation matrix to apply to resulting shapes. + cleanup : bool = True + If True, removes extremely small features from each polygon's boundary. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. If ``None``, uses + high-quality visualization settings. Returns ------- @@ -3089,17 +3193,26 @@ def intersections_tilted_plane( For more details refer to `Shapely's Documentation `_. """ - a = self.geometry_a.intersections_tilted_plane(normal, origin, to_2D) - b = self.geometry_b.intersections_tilted_plane(normal, origin, to_2D) + a = self.geometry_a.intersections_tilted_plane( + normal, origin, to_2D, cleanup=cleanup, quad_segs=quad_segs + ) + b = self.geometry_b.intersections_tilted_plane( + normal, origin, to_2D, cleanup=cleanup, quad_segs=quad_segs + ) geom_a = shapely.unary_union([Geometry.evaluate_inf_shape(g) for g in a]) geom_b = shapely.unary_union([Geometry.evaluate_inf_shape(g) for g in b]) return ClipOperation.to_polygon_list( self._shapely_operation(geom_a, geom_b), - cleanup=True, + cleanup=cleanup, ) def intersections_plane( - self, x: Optional[float] = None, y: Optional[float] = None, z: Optional[float] = None + self, + x: Optional[float] = None, + y: Optional[float] = None, + z: Optional[float] = None, + cleanup: bool = True, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Returns list of shapely geometries at plane specified by one non-None value of x,y,z. @@ -3111,6 +3224,11 @@ def intersections_plane( Position of plane in y direction, only one of x,y,z can be specified to define plane. z : float = None Position of plane in z direction, only one of x,y,z can be specified to define plane. + cleanup : bool = True + If True, removes extremely small features from each polygon's boundary. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. If ``None``, uses + high-quality visualization settings. Returns ------- @@ -3119,13 +3237,13 @@ def intersections_plane( For more details refer to `Shapely's Documentaton `_. """ - a = self.geometry_a.intersections_plane(x, y, z) - b = self.geometry_b.intersections_plane(x, y, z) + a = self.geometry_a.intersections_plane(x, y, z, cleanup=cleanup, quad_segs=quad_segs) + b = self.geometry_b.intersections_plane(x, y, z, cleanup=cleanup, quad_segs=quad_segs) geom_a = shapely.unary_union([Geometry.evaluate_inf_shape(g) for g in a]) geom_b = shapely.unary_union([Geometry.evaluate_inf_shape(g) for g in b]) return ClipOperation.to_polygon_list( self._shapely_operation(geom_a, geom_b), - cleanup=True, + cleanup=cleanup, ) @cached_property @@ -3283,7 +3401,12 @@ def bounds(self) -> Bound: ) def intersections_tilted_plane( - self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4 + self, + normal: Coordinate, + origin: Coordinate, + to_2D: MatrixReal4x4, + cleanup: bool = True, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Return a list of shapely geometries at the plane specified by normal and origin. @@ -3295,6 +3418,11 @@ def intersections_tilted_plane( Vector defining the plane origin. to_2D : MatrixReal4x4 Transformation matrix to apply to resulting shapes. + cleanup : bool = True + If True, removes extremely small features from each polygon's boundary. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. If ``None``, uses + high-quality visualization settings. Returns ------- @@ -3306,11 +3434,18 @@ def intersections_tilted_plane( return [ intersection for geometry in self.geometries - for intersection in geometry.intersections_tilted_plane(normal, origin, to_2D) + for intersection in geometry.intersections_tilted_plane( + normal, origin, to_2D, cleanup=cleanup, quad_segs=quad_segs + ) ] def intersections_plane( - self, x: Optional[float] = None, y: Optional[float] = None, z: Optional[float] = None + self, + x: Optional[float] = None, + y: Optional[float] = None, + z: Optional[float] = None, + cleanup: bool = True, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Returns list of shapely geometries at plane specified by one non-None value of x,y,z. @@ -3322,6 +3457,11 @@ def intersections_plane( Position of plane in y direction, only one of x,y,z can be specified to define plane. z : float = None Position of plane in z direction, only one of x,y,z can be specified to define plane. + cleanup : bool = True + If True, removes extremely small features from each polygon's boundary. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. If ``None``, uses + high-quality visualization settings. Returns ------- @@ -3335,7 +3475,9 @@ def intersections_plane( return [ intersection for geometry in self.geometries - for intersection in geometry.intersections_plane(x=x, y=y, z=z) + for intersection in geometry.intersections_plane( + x=x, y=y, z=z, cleanup=cleanup, quad_segs=quad_segs + ) ] def intersects_axis_position(self, axis: float, position: float) -> bool: diff --git a/tidy3d/components/geometry/mesh.py b/tidy3d/components/geometry/mesh.py index 5cb4cba58c..db0b4560fa 100644 --- a/tidy3d/components/geometry/mesh.py +++ b/tidy3d/components/geometry/mesh.py @@ -530,7 +530,12 @@ def bounds(self) -> Bound: return self.trimesh.bounds def intersections_tilted_plane( - self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4 + self, + normal: Coordinate, + origin: Coordinate, + to_2D: MatrixReal4x4, + cleanup: bool = True, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Return a list of shapely geometries at the plane specified by normal and origin. @@ -542,6 +547,10 @@ def intersections_tilted_plane( Vector defining the plane origin. to_2D : MatrixReal4x4 Transformation matrix to apply to resulting shapes. + cleanup : bool = True + If True, removes extremely small features from each polygon's boundary. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. Not used for TriangleMesh. Returns ------- @@ -557,7 +566,12 @@ def intersections_tilted_plane( return path.polygons_full def intersections_plane( - self, x: Optional[float] = None, y: Optional[float] = None, z: Optional[float] = None + self, + x: Optional[float] = None, + y: Optional[float] = None, + z: Optional[float] = None, + cleanup: bool = True, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Returns list of shapely geometries at plane specified by one non-None value of x,y,z. @@ -569,6 +583,10 @@ def intersections_plane( Position of plane in y direction, only one of x,y,z can be specified to define plane. z : float = None Position of plane in z direction, only one of x,y,z can be specified to define plane. + cleanup : bool = True + If True, removes extremely small features from each polygon's boundary. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. Not used for TriangleMesh. Returns ------- @@ -623,7 +641,7 @@ def intersections_plane( "Using bounding box instead." ) log.warning(f"Error encountered: {e}") - return self.bounding_box.intersections_plane(x=x, y=y, z=z) + return self.bounding_box.intersections_plane(x=x, y=y, z=z, cleanup=cleanup) def inside( self, x: np.ndarray[float], y: np.ndarray[float], z: np.ndarray[float] diff --git a/tidy3d/components/geometry/polyslab.py b/tidy3d/components/geometry/polyslab.py index 33c30236ce..1edcd87805 100644 --- a/tidy3d/components/geometry/polyslab.py +++ b/tidy3d/components/geometry/polyslab.py @@ -587,7 +587,11 @@ def _move_axis_reverse(arr: NDArray) -> NDArray: @verify_packages_import(["trimesh"]) def _do_intersections_tilted_plane( - self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4 + self, + normal: Coordinate, + origin: Coordinate, + to_2D: MatrixReal4x4, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Return a list of shapely geometries at the plane specified by normal and origin. @@ -599,6 +603,8 @@ def _do_intersections_tilted_plane( Vector defining the plane origin. to_2D : MatrixReal4x4 Transformation matrix to apply to resulting shapes. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. Not used for PolySlab geometry. Returns ------- @@ -641,7 +647,7 @@ def _do_intersections_tilted_plane( path, _ = section.to_2D(to_2D=to_2D) return path.polygons_full - def _intersections_normal(self, z: float) -> list[Shapely]: + def _intersections_normal(self, z: float, quad_segs: Optional[int] = None) -> list[Shapely]: """Find shapely geometries intersecting planar geometry with axis normal to slab. Parameters @@ -2567,7 +2573,12 @@ def _dilation_value_at_reference_to_coord(self, dilation: float) -> float: return z_coord def intersections_tilted_plane( - self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4 + self, + normal: Coordinate, + origin: Coordinate, + to_2D: MatrixReal4x4, + cleanup: bool = True, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Return a list of shapely geometries at the plane specified by normal and origin. @@ -2579,6 +2590,10 @@ def intersections_tilted_plane( Vector defining the plane origin. to_2D : MatrixReal4x4 Transformation matrix to apply to resulting shapes. + cleanup : bool = True + If True, removes extremely small features from each polygon's boundary. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. Not used for PolySlab. Returns ------- @@ -2592,7 +2607,9 @@ def intersections_tilted_plane( [ base.Geometry.evaluate_inf_shape(shape) for polyslab in self.sub_polyslabs - for shape in polyslab.intersections_tilted_plane(normal, origin, to_2D) + for shape in polyslab.intersections_tilted_plane( + normal, origin, to_2D, cleanup=cleanup, quad_segs=quad_segs + ) ] ) ] diff --git a/tidy3d/components/geometry/primitives.py b/tidy3d/components/geometry/primitives.py index 667ef5cb1a..004d820b66 100644 --- a/tidy3d/components/geometry/primitives.py +++ b/tidy3d/components/geometry/primitives.py @@ -28,7 +28,7 @@ _N_SAMPLE_CURVE_SHAPELY = 40 # for shapely circular shapes discretization in visualization -_N_SHAPELY_QUAD_SEGS = 200 +_N_SHAPELY_QUAD_SEGS_VISUALIZATION = 200 # Default number of points to discretize polyslab in `Cylinder.to_polyslab()` _N_PTS_CYLINDER_POLYSLAB = 51 @@ -71,7 +71,12 @@ def inside( return (dist_x**2 + dist_y**2 + dist_z**2) <= (self.radius**2) def intersections_tilted_plane( - self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4 + self, + normal: Coordinate, + origin: Coordinate, + to_2D: MatrixReal4x4, + cleanup: bool = True, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Return a list of shapely geometries at the plane specified by normal and origin. @@ -83,6 +88,11 @@ def intersections_tilted_plane( Vector defining the plane origin. to_2D : MatrixReal4x4 Transformation matrix to apply to resulting shapes. + cleanup : bool = True + If True, removes extremely small features from each polygon's boundary. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. If ``None``, uses + ``_N_SHAPELY_QUAD_SEGS_VISUALIZATION`` for high-quality visualization. Returns ------- @@ -91,6 +101,9 @@ def intersections_tilted_plane( For more details refer to `Shapely's Documentation `_. """ + if quad_segs is None: + quad_segs = _N_SHAPELY_QUAD_SEGS_VISUALIZATION + normal = np.array(normal) unit_normal = normal / (np.sum(normal**2) ** 0.5) projection = np.dot(np.array(origin) - np.array(self.center), unit_normal) @@ -106,13 +119,18 @@ def intersections_tilted_plane( u /= np.sum(u**2) ** 0.5 v = np.cross(unit_normal, u) - angles = np.linspace(0, 2 * np.pi, _N_SHAPELY_QUAD_SEGS * 4 + 1)[:-1] + angles = np.linspace(0, 2 * np.pi, quad_segs * 4 + 1)[:-1] circ = center + np.outer(np.cos(angles), radius * u) + np.outer(np.sin(angles), radius * v) vertices = np.dot(np.hstack((circ, np.ones((angles.size, 1)))), to_2D.T) return [shapely.Polygon(vertices[:, :2])] def intersections_plane( - self, x: Optional[float] = None, y: Optional[float] = None, z: Optional[float] = None + self, + x: Optional[float] = None, + y: Optional[float] = None, + z: Optional[float] = None, + cleanup: bool = True, + quad_segs: Optional[int] = None, ) -> list[BaseGeometry]: """Returns shapely geometry at plane specified by one non None value of x,y,z. @@ -124,14 +142,22 @@ def intersections_plane( Position of plane in x direction, only one of x,y,z can be specified to define plane. z : float = None Position of plane in x direction, only one of x,y,z can be specified to define plane. + cleanup : bool = True + If True, removes extremely small features from each polygon's boundary. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. If ``None``, uses + ``_N_SHAPELY_QUAD_SEGS_VISUALIZATION`` for high-quality visualization. Returns ------- List[shapely.geometry.base.BaseGeometry] List of 2D shapes that intersect plane. For more details refer to - `Shapely's Documentation `_. + `Shapely's Documentation ``. """ + if quad_segs is None: + quad_segs = _N_SHAPELY_QUAD_SEGS_VISUALIZATION + axis, position = self.parse_xyz_kwargs(x=x, y=y, z=z) if not self.intersects_axis_position(axis, position): return [] @@ -139,7 +165,7 @@ def intersections_plane( intersect_dist = self._intersect_dist(position, z0) if not intersect_dist: return [] - return [shapely.Point(x0, y0).buffer(0.5 * intersect_dist, quad_segs=_N_SHAPELY_QUAD_SEGS)] + return [shapely.Point(x0, y0).buffer(0.5 * intersect_dist, quad_segs=quad_segs)] @cached_property def bounds(self) -> Bound: @@ -430,7 +456,11 @@ def _update_from_bounds(self, bounds: tuple[float, float], axis: Axis) -> Cylind @verify_packages_import(["trimesh"]) def _do_intersections_tilted_plane( - self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4 + self, + normal: Coordinate, + origin: Coordinate, + to_2D: MatrixReal4x4, + quad_segs: Optional[int] = None, ) -> list[Shapely]: """Return a list of shapely geometries at the plane specified by normal and origin. @@ -442,6 +472,9 @@ def _do_intersections_tilted_plane( Vector defining the plane origin. to_2D : MatrixReal4x4 Transformation matrix to apply to resulting shapes. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. If ``None``, uses + ``_N_SHAPELY_QUAD_SEGS_VISUALIZATION`` for high-quality visualization. Returns ------- @@ -452,6 +485,9 @@ def _do_intersections_tilted_plane( """ import trimesh + if quad_segs is None: + quad_segs = _N_SHAPELY_QUAD_SEGS_VISUALIZATION + z0, (x0, y0) = self.pop_axis(self.center, self.axis) half_length = self.finite_length_axis / 2 @@ -471,7 +507,7 @@ def _do_intersections_tilted_plane( r_bot = 0 z_bot = z0 + self._radius_z(z0) / self._tanq - angles = np.linspace(0, 2 * np.pi, _N_SHAPELY_QUAD_SEGS * 4 + 1) + angles = np.linspace(0, 2 * np.pi, quad_segs * 4 + 1) if r_bot > 0: x_bot = x0 + r_bot * np.cos(angles) @@ -525,13 +561,18 @@ def _do_intersections_tilted_plane( path, _ = section.to_2D(to_2D=to_2D) return path.polygons_full - def _intersections_normal(self, z: float) -> list[BaseGeometry]: + def _intersections_normal( + self, z: float, quad_segs: Optional[int] = None + ) -> list[BaseGeometry]: """Find shapely geometries intersecting cylindrical geometry with axis normal to slab. Parameters ---------- z : float Position along the axis normal to slab + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. If ``None``, uses + ``_N_SHAPELY_QUAD_SEGS_VISUALIZATION`` for high-quality visualization. Returns ------- @@ -540,6 +581,8 @@ def _intersections_normal(self, z: float) -> list[BaseGeometry]: For more details refer to `Shapely's Documentation `_. """ + if quad_segs is None: + quad_segs = _N_SHAPELY_QUAD_SEGS_VISUALIZATION static_self = self.to_static() @@ -550,7 +593,7 @@ def _intersections_normal(self, z: float) -> list[BaseGeometry]: return [] _, (x0, y0) = self.pop_axis(static_self.center, axis=self.axis) - return [shapely.Point(x0, y0).buffer(radius_offset, quad_segs=_N_SHAPELY_QUAD_SEGS)] + return [shapely.Point(x0, y0).buffer(radius_offset, quad_segs=quad_segs)] def _intersections_side(self, position: float, axis: int) -> list[BaseGeometry]: """Find shapely geometries intersecting cylindrical geometry with axis orthogonal to length. diff --git a/tidy3d/components/geometry/utils.py b/tidy3d/components/geometry/utils.py index 7fc847e41b..9a483278ec 100644 --- a/tidy3d/components/geometry/utils.py +++ b/tidy3d/components/geometry/utils.py @@ -93,6 +93,8 @@ def merging_geometries_on_plane( plane: Box, property_list: list[Any], interior_disjoint_geometries: bool = False, + cleanup: bool = True, + quad_segs: Optional[int] = None, ) -> list[tuple[Any, Shapely]]: """Compute list of shapes on plane. Overlaps are removed or merged depending on provided property_list. @@ -107,6 +109,11 @@ def merging_geometries_on_plane( Property value for each structure. interior_disjoint_geometries: bool = False If ``True``, geometries of different properties on the plane must not be overlapping. + cleanup : bool = True + If True, removes extremely small features from each polygon's boundary. + quad_segs : Optional[int] = None + Number of segments used to discretize circular shapes. If ``None``, uses + high-quality visualization settings. Returns ------- @@ -122,7 +129,7 @@ def merging_geometries_on_plane( shapes = [] for geo, prop in zip(geometries, property_list): # get list of Shapely shapes that intersect at the plane - shapes_plane = plane.intersections_with(geo) + shapes_plane = plane.intersections_with(geo, cleanup=cleanup, quad_segs=quad_segs) # Append each of them and their property information to the list of shapes for shape in shapes_plane: diff --git a/tidy3d/components/grid/corner_finder.py b/tidy3d/components/grid/corner_finder.py index 59386c5595..6cf8b34339 100644 --- a/tidy3d/components/grid/corner_finder.py +++ b/tidy3d/components/grid/corner_finder.py @@ -16,6 +16,10 @@ from tidy3d.constants import inf CORNER_ANGLE_THRESOLD = 0.25 * np.pi +# For shapely circular shapes discretization. +N_SHAPELY_QUAD_SEGS = 8 +# whether to clean tiny features that sometimes occurs in shapely operations +SHAPELY_CLEANUP = False class CornerFinderSpec(Tidy3dBaseModel): @@ -132,7 +136,12 @@ def _merged_pec_on_plane( medium_list = [PEC for _ in geometry_list] # merge geometries merged_geos = merging_geometries_on_plane( - geometry_list, plane, medium_list, interior_disjoint_geometries + geometry_list, + plane, + medium_list, + interior_disjoint_geometries, + cleanup=SHAPELY_CLEANUP, + quad_segs=N_SHAPELY_QUAD_SEGS, ) return merged_geos diff --git a/tidy3d/components/grid/grid_spec.py b/tidy3d/components/grid/grid_spec.py index bcffe6b2f3..a5e3f76404 100644 --- a/tidy3d/components/grid/grid_spec.py +++ b/tidy3d/components/grid/grid_spec.py @@ -1922,8 +1922,8 @@ def _find_vertical_intersections( cells_valid = [] # for each polygon vertex find the index of the first grid line on the right - grid_lines_on_right = np.argmax(grid_x_coords[:, None] >= poly_vertices[None, :, 0], axis=0) - grid_lines_on_right[poly_vertices[:, 0] >= grid_x_coords[-1]] = len(grid_x_coords) + # Use searchsorted for O(n log m) instead of O(n * m) with argmax broadcasting + grid_lines_on_right = np.searchsorted(grid_x_coords, poly_vertices[:, 0], side="left") # once we know these indices then we can find grid lines intersected by the i-th # segment of the polygon as # [grid_lines_on_right[i], grid_lines_on_right[i+1]) for grid_lines_on_right[i] > grid_lines_on_right[i+1] diff --git a/tidy3d/components/grid/mesher.py b/tidy3d/components/grid/mesher.py index f66be82606..589e4a82a8 100644 --- a/tidy3d/components/grid/mesher.py +++ b/tidy3d/components/grid/mesher.py @@ -2,6 +2,7 @@ from __future__ import annotations +import bisect import logging import warnings from abc import ABC, abstractmethod @@ -271,7 +272,13 @@ def parse_structures( # Rtree from the 2D part of the bounding boxes tree = self.bounds_2d_tree(struct_bbox) - intervals = {"coords": list(domain_bounds), "structs": [[]]} + intervals = { + "coords": list(domain_bounds), + "structs": [[]], + "min_steps": [ + np.inf + ], # Track min step per interval for O(1) lookup; inf until structures added + } """ Build the ``intervals`` dictionary. ``intervals["coords"]`` gets populated based on the bounding boxes of all structures in the list (some filtering is done to exclude points that will be too close together compared to the absolute lower desired step). At every point, the @@ -364,6 +371,7 @@ def parse_structures( in_domain = np.nonzero((coords >= domain_bounds[0]) * (coords <= domain_bounds[1]))[0] intervals["coords"] = [intervals["coords"][int(i)] for i in in_domain] intervals["structs"] = [intervals["structs"][int(i)] for i in in_domain if i < num_ints] + intervals["min_steps"] = [intervals["min_steps"][int(i)] for i in in_domain if i < num_ints] # Compute the maximum allowed step size in each interval max_steps = [] @@ -372,11 +380,9 @@ def parse_structures( if max(intervals["structs"][coord_ind]) >= num_unenforced: max_step = structure_steps[max(intervals["structs"][coord_ind])] max_steps.append(max_step) - # otherwise, define the max step as the minimum over all medium steps - # of media in this interval + # otherwise, use the cached minimum step for this interval - O(1) lookup else: - max_step = np.amin(structure_steps[intervals["structs"][coord_ind]]) - max_steps.append(max_step) + max_steps.append(intervals["min_steps"][coord_ind]) # Re-evaluate the absolute smallest min_step and remove intervals that are smaller than that intervals["coords"], max_steps = self.filter_min_step(intervals["coords"], max_steps) @@ -435,13 +441,14 @@ def insert_bbox( coords = intervals["coords"] structs = intervals["structs"] + cached_min_steps = intervals["min_steps"] min_step_check = MIN_STEP_SCALE * min_step # Left structure bound bound_coord = str_bbox[0, 2] - indsmin = np.nonzero(bound_coord <= coords)[0] - indmin = int(indsmin[0]) # coordinate is in interval index ``indmin - 1`` + # coordinate is in interval index ``indmin - 1`` + indmin = bisect.bisect_left(coords, bound_coord) is_close_l = self.is_close(bound_coord, coords, indmin - 1, min_step_check) is_close_r = self.is_close(bound_coord, coords, indmin, min_step_check) is_contained = self.is_contained(bound_coord, bbox_contained_2d) @@ -450,9 +457,8 @@ def insert_bbox( skip_unshadowed = False if unshadowed and indmin > 0: grid_size_str = structure_steps[str_ind] - min_grid_size = min( - (structure_steps[ind] for ind in structs[indmin - 1]), default=grid_size_str - ) + # Use cached min_steps for O(1) lookup instead of np.amin + min_grid_size = cached_min_steps[indmin - 1] if not (isclose(grid_size_str, min_grid_size) or grid_size_str < min_grid_size): skip_unshadowed = True @@ -468,11 +474,13 @@ def insert_bbox( # Copy the structure containment list to the newly created interval struct_list = structs[max(0, indmin - 1)] structs.insert(indmin, struct_list.copy()) + # Copy the cached min step to the newly created interval + cached_min_steps.insert(indmin, cached_min_steps[max(0, indmin - 1)]) # Right structure bound bound_coord = str_bbox[1, 2] - indsmax = np.nonzero(bound_coord >= coords)[0] - indmax = int(indsmax[-1]) # coordinate is in interval index ``indmax`` + # coordinate is in interval index ``indmax`` + indmax = bisect.bisect_right(coords, bound_coord) - 1 is_close_l = self.is_close(bound_coord, coords, indmax, min_step_check) is_close_r = self.is_close(bound_coord, coords, indmax + 1, min_step_check) is_contained = self.is_contained(bound_coord, bbox_contained_2d) @@ -481,9 +489,8 @@ def insert_bbox( skip_unshadowed = False if unshadowed and indmax < len(structs): grid_size_str = structure_steps[str_ind] - min_grid_size = min( - (structure_steps[ind] for ind in structs[indmax]), default=grid_size_str - ) + # Use cached min_steps for O(1) lookup instead of np.amin + min_grid_size = cached_min_steps[indmax] if not (isclose(grid_size_str, min_grid_size) or grid_size_str < min_grid_size): skip_unshadowed = True @@ -498,14 +505,35 @@ def insert_bbox( # Copy the structure containment list to the newly created interval struct_list = structs[min(indmax - 1, len(structs) - 1)] structs.insert(indmax, struct_list.copy()) + # Copy the cached min step to the newly created interval + cached_min_steps.insert( + indmax, cached_min_steps[min(indmax - 1, len(cached_min_steps) - 1)] + ) # Add the current structure index to all intervals that it spans, if it is not # contained in any of the latter structures - for interval_ind in range(indmin, indmax): - # Check at the midpoint to avoid numerical issues at the interval boundaries - mid_coord = (coords[interval_ind] + coords[interval_ind + 1]) / 2 - if not self.is_contained(mid_coord, bbox_contained_2d): - structs[interval_ind].append(str_ind) + current_step = structure_steps[str_ind] + num_intervals = indmax - indmin + if num_intervals > 0: + # Vectorized containment check for better performance with many intervals + if bbox_contained_2d: + # Compute all midpoints at once + mid_coords = np.array( + [(coords[i] + coords[i + 1]) / 2 for i in range(indmin, indmax)] + ) + # Vectorized containment: check all midpoints against all containing boxes + contained_mask = self._is_contained_vectorized(mid_coords, bbox_contained_2d) + else: + # No containing boxes - nothing is contained + contained_mask = np.zeros(num_intervals, dtype=bool) + + # Add structure to non-contained intervals + for i, interval_ind in enumerate(range(indmin, indmax)): + if not contained_mask[i]: + structs[interval_ind].append(str_ind) + # Update cached min step incrementally - O(1) instead of recomputing + if current_step < cached_min_steps[interval_ind]: + cached_min_steps[interval_ind] = current_step return indmin >= indmax @@ -760,6 +788,40 @@ def is_contained(normal_pos: float, contained_2d: list[ArrayFloat1D]) -> bool: contain_box[0, 2] <= normal_pos <= contain_box[1, 2] for contain_box in contained_2d ) + @staticmethod + def _is_contained_vectorized( + normal_positions: ArrayFloat1D, contained_2d: list[ArrayFloat1D] + ) -> ArrayFloat1D: + """Vectorized containment check for multiple positions at once. + + Parameters + ---------- + normal_positions : ArrayFloat1D + Array of positions along the meshing direction to check. + contained_2d : list[ArrayFloat1D] + List of bounding boxes that could contain the positions. + + Returns + ------- + ArrayFloat1D + Boolean array where True means the position is contained in at least one box. + """ + if not contained_2d: + return np.zeros(len(normal_positions), dtype=bool) + + # Stack all z-bounds: shape (num_boxes, 2) + z_bounds = np.array([[box[0, 2], box[1, 2]] for box in contained_2d]) + + # Check containment: positions[:, None] broadcasts against z_bounds[None, :, :] + # Result shape: (num_positions, num_boxes) + positions = np.asarray(normal_positions) + in_range = (positions[:, None] >= z_bounds[None, :, 0]) & ( + positions[:, None] <= z_bounds[None, :, 1] + ) + + # Any box containing the position means it's contained + return np.any(in_range, axis=1) + @staticmethod def filter_min_step( interval_coords: list[float], max_steps: list[float]